A New Kilohertz Gravitational-wave Feature from Rapidly Rotating Core-collapse Supernovae

We present self-consistent three-dimensional core-collapse supernova simulations of a rotating 20M ⊙ progenitor model with various initial angular velocities from 0.0 to 4.0 rad s−1 using the smoothed particle hydrodynamics code SPHYNX and the grid-based hydrodynamics code FLASH. We identify two strong gravitational-wave features with peak frequencies of ∼300 Hz and ∼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 only present 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 antineutrinos. 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 initial peak frequency is around 800 Hz, and it 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 for estimating the initial angular velocity of a collapsing core.

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 generalrelativistic (GR) hydrodynamical simulations that include a wide range of angular velocities of a 12M ⊙ progenitor.However, due to limits on computational resources, Richers et al. (2017) calculated only up to the first ∼ 50 ms after the core bounce, focused on the bounce signals, and conducted the calculations exclusively with the s12 progenitor.Pajkos et al. (2019Pajkos et al. ( , 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 gridbased 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 ⊙ progen-itor, 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.

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.
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 codecomparison 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 FLASH2 that includes the IDSA for neutrino transport (Pan et al. 2016(Pan et al. , 2018(Pan et al. , 2019(Pan et al. , 2021)).We also use the same effective GR potential and LS220 equation of state.

Initial Setup
In this work, we employ the same 20M ⊙ 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 ext = 1.8 × 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 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 electronflavor neutrinos.The heavier µ and τ neutrinos and antineutrinos 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, where Ω is the angular velocity as a function of the spherical radius r, A = 1.03376 × 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 Ω 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 200k, 675k, and 1600k 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 1600k particles.Furthermore, we used the simulations with 675k particles to investigate a larger parameter space of Ω 0 .
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 ∼ 0 • .2− 0 • .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 µ and τ neutrinos.We use the same s20 progenitor from Woosley & Heger (2007) and take the same rotational profile as in Equation (1) with Ω 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, Ω 0 .Finally, L and H stand for low and high resolution, respectively, for SPHYNX models.The column f peak shows the resulting peak GW frequency around kHz in the time interval t 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 ρ ≥ 10 6 g cm −3 at the initial time.The last column ∆x 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.

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: where G is gravitational constant, c is the speed of light, and D is the distance to the source.The reduced quadrupole (I lm ) is defined in Cartesian coordinates as, where l and m are components of the position vector x, δ lm is a δ-Kronecker, and ρ is the local density.The components of the reduced quadrupole in spherical coordinates are related to the Cartesian by (Oohara et al. 1997), Therefore, to obtain the GW waveforms we compute the six non-zero components of Ï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 Ïlm .Nevertheless, time derivatives cause numerical difficulties due to two main reasons: the numerical noise introduced by the discretization and the magnification of the highfrequency 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 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 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(Pan et al. , 2021)).

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 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 /ρ i , where m i and ρ i are the mass and density carried by the SPH particles, respectively.The orthonormal harmonic basis functions, Y m l , are expressed as and P m l (cos θ) 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.

Dynamics Overview
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 Ω 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 (∆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 Ω 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 pb ∼ 40 and 80 ms, respectively, transitioning from infall to outflow.The Ω 0 = 3 cases (dashed lines) show behavior similar to that of the Ω 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 (∆x ∼ 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 Ω 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 Ω 0 = 3 rad s −1 , the deviation in the PNS mass M 11 increases further after t 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 Ω 0 = 3 rad s −1 , when angular momentum keeps propagat-  ing 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 pns ∼ 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 sphericallyaveraged 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.

GW Features and Origin in Rapidly Rotating CCSNe
In this section, we discuss the GW features in the models with Ω 0 = 2 and 3 rad s −1 , using SPHYNX with 1600k 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 Py-CWT3 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.
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 ringdown signals, we can observe that the GW emissions exhibit quasi-sinusoidal time oscillations starting at t 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. 2005Ott et al. , 2007;;Scheidegger et al. 2008Scheidegger et al. , 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 angu- 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 Ω 0 = 2 rad s −1 (F20) and Ω 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 ρ < 10 11 , 10 11 ≤ ρ < 10 13 , and ρ ≥ 10 13 g cm −3 for models S30H and F30, seen along the pole to eliminate bounce and ringdown 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 ≤ ρ < 10 13 g cm −3 , while the kHz signal is mainly from the PNS inner core where ρ ≥ 10 13 g cm −3 .
In Figure 5, we also evaluate the spherical harmonic components a 11 and a 22 in the regions of 10 11 ≤ ρ < 10 13 and ρ ≥ 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 mode,m .The mode frequency is related to the pattern frequency by f pat,m = f mode,m /m (Watts et al. 2005).Note that the mode frequency of a 11 is doubled in Figure 5 2020).On the other hand, the kHz signal resembles the ∼ 930 Hz signal found by Ott et al. (2007) in CCSN simulations of a 20M ⊙ 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 log 10 (h + ) Figure 5. GW spectrograms of the plus mode emitted from regions with density ρ < 10 11 g cm −3 (left), 10 11 ≤ ρ < 10 13 g cm −3 (middle), and ρ ≥ 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 a11 and a22, respectively, where the frequency of a11 is doubled to ease comparison.the 300 Hz and kHz signals, but the interactions among these GW features are more complex than that in model S30H.
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) where ϖ 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 ϖ = 15 km.We note that the relative difference in the mean density, C 0 , is below 30% before t pb ∼ 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 pb ≲ 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 pb = 30 ms.We evaluate the Lamb frequency via and the Brunt-Väisälä frequency using the Ledoux criterion (Buras et al. 2006;Ott et al. 2013) where l is taken to be 1, c s is the local speed of sound, dΦ is the local gravitational potential where the approximation dΦ/dr ∼ −GM (r)/r 2 is adopted.The Ledoux criterion reads (Ledoux 1947) (16) where we approximate the lepton fraction by the electron fraction, Y l ∼ 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 ∼ 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 O(v/c) terms of the neutrino pressure are included (see Equation 24in 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): where ⟨E νe ⟩ = Z νe /Y νe is the mean energy of trapped electron anti-neutrinos.We use the mean energy ⟨E νe ⟩, instead of Y νe or Z νe , to avoid the sharp decline in electron antineutrino fractions, which could introduce strong numerical noise.We evaluate the spherical harmonic components in the shell R 13 ± 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  Figure 8. Spectrograms of the normalized spherical harmonic coefficients, a11 (white, doubled frequency) and a22 (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.
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.

GW Dependence on Initial Angular Velocity
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 675k 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 dyn ∼ 2 R 3 13 /GM 13 /2π, for the domain with ρ ≥ 10 13 g cm −3 for models with Ω 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 Ω 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 675k particles accurately capture the main GW features observed in our simulations with 1600k 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 675k 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.
As a result, from Figure 10 we find that the 300 Hz signal starts to appear when the initial angular velocity Ω 0 ≥ 1 rad s −1 and it persists consistently across different angular velocities.However, the kHz signal appears only in a range within 1.5 ≤ Ω 0 ≤ 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 ρ > 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, (ρ − ρ)/ρ, at t pb = 50 ms in Figure 11, where ρ 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 ≲ ρ ≲ 10 13 g cm −3 , all models with Ω 0 ≥ 1.0 rad s −1 show the spiral structure in the density fluctuation plot, but not in models with Ω 0 ≤ 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 ρ > 10 13 g cm −3 , where the kHz signal emanated, the m = 1 deformation develops only in the models with 1.5 ≤ Ω 0 ≤ 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 ρ > 10 13 g cm −3 .
Figure 12 shows the mean energy of trapped electron antineutrinos 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 ∼ 20 km in the models with 1.5 ≤ Ω 0 ≤ 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 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 (∼ 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.In this section, we discuss the detectability of the GW features from rapid-rotating CCSNe using the current groundbased GW detectors.Figure 13 shows the amplitude spectral density (ASD) of the plus mode of GW emissions between t pb = −10 and 100 ms in the SPHYNX simulations with 675k 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.
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 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 ≤ Ω 0 ≤ 3.5 rad s −1 .The cases with Ω 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 ≤ Ω 0 ≤ 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 Ω 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 pb = −10 and 100 ms.We can see that the GW emissions have a qualitatively similar spectral density distribution between the SPHYNX simulations with 200k (S30L), 675k (S30), and 1600k 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.

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 ⊙ progenitor with different initial angular velocities (Ω 0 ), which are parametrically added to the progenitor model.In this work, we computed models with 0 ≤ Ω 0 ≤ 4 rad s −1 using the smoothed particle hydrodynamics code SPHYNX, and models with Ω 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 Ω 0 ≥ 1.0 rad s −1 , which mainly emanated from the region of 10 11 ≤ ρ < 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 ≤ Ω 0 ≤ 3.5 rad s −1 , originates mainly in the region of ρ ≥ 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 ≤ ρ < 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 ≤ Ω 0 ≤ 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.

Figure 1 .
Figure1.Panels from left to right and from top to bottom describe the time evolution of the central density (ρc), the averaged shock radius (R sh ), the mass accretion rate measured at r = 200 km and 500 km ( Ṁ200 and Ṁ500), the enclosed mass within a density contour of 10 11 g cm −3 and 10 13 g cm −3 (M11 and M13), and the averaged isodensity radii corresponding to M11 and M13 (R11 and R13) in the models S20H, S30H, F20, and F30.

Figure 2 .
Figure 2. Entropy distribution on the XZ-plane at t pb = 50 ms in models S30H (top) and F30 (bottom).The arrows represent the velocity field.

Figure 3 .
Figure3.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.lar 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 Figure4, 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 Figure4, we present the GW strains obtained using the FLASH code with Ω 0 = 2 rad s −1 (F20) and Ω 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

Figure 4 .
Figure4.GW strains and spectrograms of the plus mode for the models with Ω0 = 2 and 3 rad s −1 , using SPHYNX with 1600k 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).
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 pat,1 ≃ f 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 ⊙ 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. (

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

Figure 7 .
Figure 7. Radial profiles of rotational (blue), Brunt-Väisälä (orange), and Lamb (green) frequencies at t 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 9 .
Figure 9.Time evolution of the PNS mass with density ρ ≥ 10 13 g cm −3 (top), the corresponding averaged isodensity radius (middle), and the doubled dynamical frequency (bottom) in the SPHYNX simulations with 675k particles (solid lines).Different colors represent simulations with different initial angular velocities.Dashed lines show the counterpart models (S20H and S30H) with 1600k particles for comparison.

Figure 10 .
Figure10.GW strains and spectrograms of the plus mode seen along the pole at a distance of 10 kpc.Different panels represent the SPHYNX 675k particle simulations with different initial angular velocity.The white lines represent the doubled dynamical frequency for the domain with ρ ≥ 10 13 g cm −3 in each model.

Figure 11 .
Figure11.Normalized density fluctuation in the equatorial plane at t pb = 50 ms for models using SPHYNX with 675k particles and different initial angular velocities.The average density, ρ, is taken over the azimuthal direction in the plane.The black dashed curves denote the isodensity contours at ρ = 10 11 , 10 12 , 10 13 , and 10 14 g cm −3 .

Figure 12 .
Figure12.Mean energy of trapped electron anti-neutrino distribution on the XZ-plane at t pb = 50 ms for models using SPHYNX with 675k 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 ρ = 10 12 , 10 13 , and 10 14 g cm −3 .

3. 4 .Figure 13 .
Figure13.Amplitude spectral density of the plus mode of GW emissions between t pb = −10 and 100 ms, for models using SPHYNX with 675k 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 .
Figure14.Similar to Figure13, but for the Ω0 = 3.0 models using different codes and resolutions.For the comparison, only the GW emissions between t pb = −10 and 100 ms are used for analysis.

Table 1 .
Summary of the models.Columns from left to right show the model name, numerical code, initial angular velocity (Ω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 (∆xmin,20ms).