Three-dimensional GRMHD Simulations of Rapidly Rotating Stellar Core-Collapse

We present results from fully general relativistic (GR), three-dimensional (3D), neutrino-radiation magneto-hydrodynamic (MHD) simulations of stellar core collapse of a 20 M$_\odot$ star with spectral neutrino transport. Our focus is to study the gravitational-wave (GW) signatures from the magnetorotationally (MR)-driven models. By parametrically changing the initial angular velocity and the strength of the magnetic fields in the core, we compute four models. Our results show that the MHD outflows are produced only for models (two out of four), to which magnetic field strengths of 10$^{12}$ G and rotation rates of 1 or 2 rad s$^{-1}$ are initially imposed in the core. Seen from the direction perpendicular to the rotational axis, a characteristic waveform is obtained exhibiting a monotonic time increase in the wave amplitude. As previously identified, this stems from the propagating MHD outflows along the axis. We show that the GW amplitude from anisotropic neutrino emission becomes more than one order-of-magnitude bigger than that from the matter contribution, whereas seen from the rotational axis, both of the two components are in the same order-of-magnitudes. Due to the memory effect, the frequency of the neutrino GW from our full-fledged 3D-MHD models is in the range less than $\sim$10 Hz. Toward the future GW detection for a Galactic core-collapse supernova, if driven by the MR mechanism, the planned next-generation detector as DECIGO is urgently needed to catch the low-frequency signals.


INTRODUCTION
The final evolution of massive stars ends with a catastrophic collapse of the central core, which is followed by various explosive phenomena.One such phenomenon is a core-collapse supernova (CCSN), gigantic stellar explosion (see Janka et al. 2016;Müller 2016;Radice et al. 2018;Burrows et al. 2020, for recent reviews).Though not fully understood yet, it becomes almost certain that the two major mechanisms are among the best candidates.The first one is the best-studied neutrino mechanism, which is believed to account for core-collapse supernovae with a typical explosion energy of around ∼ 10 51 erg(≡ 1 Bethe = 1 B in short).The second one is known as magnetorotational (MR) explosion, which occurs when the progenitor core rotates rapidly and possesses strong magnetic fields (see Bisnovatyi-Kogan 1970;LeBlanc & Wilson 1970;Meier et al. 1976;Müller & Hillebrandt 1979, for the original idea).The MR mechanism is expected to account for a subclass of CCSNe called, namely, hypernova (HN), the observation of which presents high explosion energy of ∼ 10 B (Stritzinger et al. 2018;Nomoto et al. 2006).
The MR explosion is characterized by bipolar outflows and inherently an non-spherical phenomenon (Ardeljan et al. 2000;Kotake et al. 2004;Burrows et al. 2007;Takiwaki et al. 2009;Scheidegger et al. 2010;Winteler et al. 2012;Sawai & Yamada 2016;Bugli et al. 2020;Varma et al. 2021).The bipolar structure originates from the amplification of magnetic fields primarily due to rotational winding along the rotational axis.The amplified magnetic fields if the strength dominates over the ram pressure of the infalling matter, lead to the ejection of matter towards the axis.By efficiently converting the available differential rotational energy of the proto-neutron star (PNS), the magnetic field can be amplified up to the equipartition.This efficient conversion leads to an explosion that is typically one order of magnitude more energetic than the neutrino-driven explosion models.It is commonly reported that the recent three-dimensional (3D) neutrino-radiation magneto-hydrodynamics (MHD) models produce collimated and slightly weaker outflows than those previously obtained in 2D (Mösta et al. 2014;Obergaulinger & Aloy 2020, 2022;Kuroda et al. 2020), however, the results are acute sensitive to the initial conditions, such as progenitor mass, rotation and magnetic fields and the numerical methods (e.g., very energetic explosion is found in Aloy & Obergaulinger 2021;Obergaulinger & Aloy 2021).
One of the crucial factors for the MR explosion is the pre-collapse rotation and magnetic fields in the stellar core.In recent years, multidimensional (multi-D) hydrodynamic simulations of CCSN progenitors shortly before core collapse have been actively performed and discussed their impacts on the successful explosion of CCSNe, especially for non-rotating massive stars (Couch et al. 2015;Müller et al. 2017;Yoshida et al. 2019Yoshida et al. , 2021b;;Fields & Couch 2020, 2021;McNeill & Müller 2020;Yadav et al. 2020).For rotating massive stars, since the seminal study of Kuhlen et al. (2003) using simplified 3D simulations, the interplay between rotation and convective motion has been investigated in 2D (Arnett & Meakin 2010;Chatzopoulos et al. 2016) and in 3D (Yoshida et al. 2021a;Fields 2022;McNeill & Müller 2022).Yoshida et al. (2021a) and Fields (2022) found that the angular momentum distribution in the 3D fast rotating model changes from rigid rotation, which is often seen in 1D stellar evolution simulations of rotating stars, to roughly constant specific angular momentum primarily due to the turbulent motion in the Si/O layer.On the other hand, McNeill & Müller (2022) presented that convection steepens specific angular momentum gradients.Possible ingredients to account for the discrepancy include the differences in the initial conditions, simulation times, and the numerical schemes.Varma & Müller (2021) performed a 3D MHD simulation of convective oxygen and neon shell burning in a non-rotating star shortly before core collapse, and observed that the initial magnetic fields are strongly amplified in the oxygen shell due to convective and turbulent flows.Several 3D CCSN simulations based on 3D progenitor models showed that 3D progenitor structure before core collapse significantly affects the dynamical evolution after core collapse (Couch et al. 2015;Müller et al. 2017;Bollig et al. 2021;Vartanyan et al. 2022).These works indicate that CCSN simulations based on 3D progenitors models are essential to link progenitors to their outcomes.
Depending on the strength of progenitor rotation, rotational flattening of the unshocked core generates the burst GW signature near at bounce (Dimmelmeier et al. 2008;Abdikamalov et al. 2014;Richers et al. 2017).If the ratio of rotational kinetic energy to gravitational energy of the PNS is O(1)%, the PNS can be subject to the low-/| | instability (Shibata et al. 2003;Ott et al. 2005), which causes nonaxisymmetric PNS deformation and produces quasi-periodic variation in GWs and neutrino signals (Takiwaki & Kotake 2018;Takiwaki et al. 2021;Shibagaki et al. 2020Shibagaki et al. , 2021;;Pan et al. 2021;Bugli et al. 2023;Micchi et al. 2023).The characteristic frequency of the GW is determined by the rotational frequency since the quadrupole deformation is a significant contribution to the GW, whereas the one for the neutrino is determined by not only the rotational frequency but also the dominant non-axisymmetric deformation mode of the PNS along the azimuthal direction (but see also Takiwaki et al. (2021), where their most rapidly rotating model shows three characteristic frequencies due to the coupling between two different modes).
The pioneering studies about the GW emission from MR explosion were conducted based on multi-D simulations with simplified neutrino treatments (Yamada & Sawai (2004); Kotake et al. (2004); Obergaulinger et al. (2006a,b); Shibata et al. (2006); Cerdá-Durán et al. (2007); Scheidegger et al. (2008Scheidegger et al. ( , 2010)); Takiwaki & Kotake (2011), and see Kotake et al. (2006); Abdikamalov et al. (2022); Powell et al. (2023) for collective references therein).Because of such approximations, these studies are basically valid only up to an early post-bounce time.Recently, MHD CCSN simulations with neutrino transport have been performed.The systematic 2D MHD CCSN simulations (Jardine et al. 2022) have been performed and confirmed that the GW signals of their MR explosion model deviate from the well-known fitting formula for the  / mode oscillation of the PNS (Müller et al. 2013).This result is confirmed in the 3D MHD simulations (Powell et al. 2023) for the fitting formula (Müller et al. 2013) and for the universal relation (Torres-Forné et al. 2019, 2021;Sotani et al. 2021).Bugli et al. (2023) showed that the GW signals from their 3D MR explosion models have a broad-band spectral shape, all of which, at least, cannot be explained only by the low-/| | instability.Different from the methods used in the above simulations, Raynaud et al. (2022) have explored the GW signals from PNS convection and its dynamo using long-term 3D anelastic simulations, where neutrino transport effects were approxinately treated by means of transport coefficients of the fluid.They found that the low-frequency excess in the GW spectrogram (≲ 100 Hz) appears at the transition to the strong field dynamo regime.The numerical study on MR explosion based on sophisticated neutrino transport has just started and further investigation along this line is mandatory to cover all potential GW signatures from CCSNe.
In the present paper, we extend the study of Kuroda (2021), who presented a detailed analysis of the GW signal from the MR-driven explosion models based on 3D-GR, MR CCSN simulations of a 20  ⊙ star with spectral neutrino transport, by employing further systematic variations of the initial angular velocity and the initial magnetic field strength.In addition to the rotating strongly magnetized model (R10B12) and the rotating, ultra-strongly magnetized model (R10B13) obtained in Kuroda (2021), we newly compute two more models, one with more initial rapid rotation with strong magnetic fields (R20B12) and the second one with slower initial rotation and strong magnetic fields (R05B12).Among the four models, MHD outflows are produced for R10B12 and R20B12, for which we initially impose most rapid rotation and strong magnetic fields in this work.We present a detailed analysis of the GW properties from our full-fledged 3DGR-MHD models for the first time, in this kind.
Our paper is structured as follows.In Section 2, we briefly introduce our general relativistic neutrino-radiation magnetohydrodynamics code, including the input physics and the initial conditions implemented.In Section 3, we present the results of our simulations.We compare the dynamics, GWs, and neutrino-GWs of our models.We also explore the detectability of their GW signals.We summarize our results in Section 4.

NUMERICAL METHOD AND INITIAL MODELS
Our numerical models are obtained using the full 3D-GR neutrinoradiation ideal MHD code of Kuroda (2021).It solves the metric equations based on the BSSN formalism (c.f.Shibata & Nakamura 1995;Baumgarte & Shapiro 1999;Marronetti et al. 2008, and references therein) on a static Cartesian spatial mesh, for which we employ here the same resolution as in Kuroda (2021).The computational domain extends to 1.5 × 10 4 km from the center, in which 2:1 ratio nested boxes with 10 refinement levels are embedded.Each nested box contains 64 3 cells so that the finest resolution at the center achieves 458 m.The induction equation for the magnetic field is solved by a constrained transport (CT) method (details can be found in Evans & Hawley 1988).The GR spectral neutrino transport is based on the two-moment scheme with M1 analytical closure (c.f.Shibata et al. 2011, and references therein), with the same neutrino energy resolution of 12 bins in the range of 1-300 MeV, as was used in Kuroda (2021).The weak rates used for the collision integral of the neutrino transport equation are given in Kotake et al. (2018).For the present study, the same SFHo nuclear relativistic mean-field equation of state (EOS) of Steiner et al. (2013) is employed as in the previous study of Kuroda (2021), taking into account the electrons/positrons, and photons contributions.Further details of our model can be found in Kuroda et al. (2016), Kuroda et al. (2020), andKuroda (2021).
We employ the solar-metallicity model of the 20  ⊙ star "s20a28n" from Woosley & Heger (2007), which is frequently used in CCSN studies.Assuming a cylindrical rotation profile, we impose the initial angular momentum of the core as where   is the time component of the contravariant four-velocity, , and  0 is set as 10 8 cm.For the initial magnetic fields, we use the following purely toroidal vector  1. summary of our models.From left to right, the columns represent the model name, the initial rotation rate parameter (equation ( 1)), the initial magnetic field strength parameter (equation ( 2)), the postbounce time at the end of the simulation, the diagnostic explosion energy, the PNS mass; the average rotation rate, the poloidal and toroidal components of the average magnetic field at the surface of the PNS.The last five quantities are measured at the end of the simulation. potential: which gives nearly uniform magnetic field parallel to the z-axis for  <  0 and dipolar magnetic field for  >  0 .In this study,  0 is set as 10 8 cm.For the nuclear EOS, we use SFHo of Steiner et al. (2013).
The 3D computational domain is a cubic box with 3 × 10 4 km width, in which nested boxes with 10 refinement levels are embedded in the Cartesian coordinates.Each box contains 64 3 cells and the minimum grid size near the origin is Δ = 458 m.The neutrino energy space logarithmically covers from 1 to 300 MeV with 12 energy bins.

RESULT
In this section, we start to briefly overview the dynamical evolution of the computed models in Section 3.1.Then in Section 3.2, we focus on the evolution of the PNS.In Sections 3.3 and 3.4, we investigate the characteristics of the GWs originating from matter and neutrino and discuss their detectability.

Postbounce Dynamics
Figure 1 shows the time evolution of the maximum shock radius,  sh,max , and the diagnostic explosion energy,  exp , of our models.Following Müller et al. (2012b), we define the GR(M)HD version of the diagnostic explosion energy as where is the lapse function,  is the relativistic energy density including the magnetic energy density but excluding the rest mass contribution (see Kuroda et al. 2020, for the definition),  is the rest mass density,  is the Lorentz factor, and  is the determinant of the threedimensional spatial metric.The most slowly rotating model, R05B12 (red solid line), does not produce an MHD jet explosion nor a successful explosion in the simulation time.After the shock wave stalls, the  average shock radius keeps ∼ 200 km until the end of the simulation (the post-bounce time,  pb = 200 ms).The strongest magnetic-field model, R10B13 (green dashed line), presents prompt explosion right after bounce, but is largely different from the usual bipolar jet explosion.At core bounce, the strength of the magnetic fields at the PNS surface is already so strong to keep the shock wave expanding.This results in the prompt explosion before the development of the collimated bipolar jets.The shock surface of this model is, therefore, more roundish than jet explosion models and the high-entropy regions show peculiar small-scale fragmented structure (see Kuroda 2021, for more details).Finally, the maximum shock radius of this model reaches ∼4400 km at the final simulation time ( pb = 316 ms) and it is still growing while the explosion energy is already satu- The left two panels of Figure 2 show the snapshots of entropy on the - plane at the time when both the north and south polar jets reach ∼300 km for the MHD jet models, R20B12 (top) and R10B12 (bottom).The faster initial rotation model, i.e., R20B12, strengthens the magnetic fields faster than R10B12 does by quickly winding up the magnetic field lines along the rotational axis, which results in the earlier jet launch of R20B12 at  pb ∼60 ms than the one of R10B12 at  pb ∼100 ms.The equatorial shock surface reaches 300 km a few tens milliseconds after this time.The right two panels of Figure 2 show the snapshots of entropy on the - plane at the final simulation time for the MHD jet models, R20B12 (top) and R10B12 (bottom).As we see in Figure 1, the MHD jets reach ∼ 1 × 10 4 km in model R20B12 and ∼ 4×10 3 km in model R10B12 before the end of their simulation times.
Our jet explosion models show a slight deviation from the axisymmetric bipolar jet flow.Mösta et al. ( 2014) found a more pronounced non-axisymmetric structure, which was claimed due to the kink instability, and that instability disrupts the axial jet at  pb = 186 ms (see their Figure 1).Meanwhile, the jet in our model R20B12 is not disrupted even at  pb = 544 ms (see top right panel of Figure 2).To clarify the difference, systematic studies on dimension and perturbation are necessary as Mösta et al. ( 2014) and Bugli et al. (2021).We leave this issue to the future investigation.

PNS Evolution
Now we move on to the PNS evolution.In this paper, we define the PNS surface as the isosurface with  = 10 11 g cm −3 .The top panel of Figure 3 displays the evolution of the PNS masses for our four models.The evolution of the PNS masses depends on the initial rotation rate.The model R20B12 reaches ∼ 1.5 M ⊙ at  pb ∼ 300 ms while the models R10B12 and R10B13 do ∼ 1.6 M ⊙ at that time.After this time, their mass does not grow because the explosion takes place not only toward the polar direction, i.e., the direction of the jet, but also toward the equatorial direction and suppresses the mass accretion (see Figure 2).In the model R05B12 the growth of the PNS mass does not stop until the end of the simulation due to the short simulation time and the failure of the shock revival, but the mass is already exceeding 1.6 M ⊙ at  pb ∼ 200 ms.The faster the progenitor core rotates, the less the PNS mass growth rate is due to the stronger centrifugal forces and the higher conversion efficiency of the poloidal to toroidal magnetic fields, which could result in the early MHD jet explosion and the efficient angular momentum transport to the external layers surrounding the PNS.
The evolution of the average surface magnetic fields and average rotation rates of the PNS is shown in the middle and bottom panels of Figure 3, respectively.Here we define the average rotation rate as the moment-of-inertia-weighted average rotation rate.In the early phase, the strength of the toroidal magnetic fields (thin lines) is larger than the poloidal ones (thick lines) except for the model R10B13.But at a later time, the strength of the poloidal magnetic fields exceeds the toroidal ones.This is because, in the early phase, the field wrapping efficiently converts the poloidal magnetic fields to the toroidal ones, which results in a larger strength of the toroidal magnetic fields than the poloidal ones.However, once the Maxwell stress decelerates the rotation of the PNS to make the rotation rate quite small, the mechanism of the field wrapping does not work anymore and the toroidal component of the magnetic fields becomes dominant.Mass accretion also plays a role of angular momentum supply to the PNS, but, because the shock expansion toward all directions substantially weakens the mass accretion (see the PNS mass evolution shown in the top panel of Figure 3), the mass accretion does not significantly affect the PNS rotation in our models after their explosion.Bugli et al. (2021) does not show such inversion of the magnetic field strength.This difference is likely because the difference in the mass accretion history changes the evolution of the PNS rotation.The model L1-0 in Bugli et al. (2021) shows a long-term mass accretion until  pb ∼ 400 ms, which results in the ∼ 1.9 M ⊙ PNS.This indicates a more continuous supply of the angular momentum by the mass accretion to the PNS in their model than ours.In fact, the angular momentum of their model is not as small as our models at the end of simulations.
To see the angular momentum loss of the PNS, we plot the time evolution of angular momentum flux due to the Maxwell stress at PNS surface for our models in Figure 4.Here we separately measure the angular momentum loss in the equatorial region, defined as /4 <  < 3/4, and in the polar region, defined as 0 <  < 3/4 or 0 <  < .Except for the model R10B13, the angular momentum loss commonly occurs mainly in the equatorial region (thin lines) and becomes maximum at  pb ∼ 60 ms.In the model R10B13, the large amount of the angular momentum of the PNS is instantaneously transferred right after core bounce toward all directions due to the large magnetic fields.Therefore, the PNS of this model slowly rotates after the bounce, and this model has no chance to generate large R20B12 R10B12 Figure 5. GW strains of plus (red solid lines) and cross (blue dashed lines) modes and spectrograms of their characteristic strains for models R20Bp12 (top panels) and R10Bp12 (bottom panels) seen along the equator (left panels) and the pole (right panels) at a source distance of 10 kpc.toroidal magnetic fields after bounce.To support that the spin down, e.g.observed in model R20B12 (see the bottom panel in Figure 3), is mainly caused by the angular momentum transfer via Maxwell stress, we compare the initial total angular momentum and the total angular momentum lost.From Figure 4, we can estimate the angular momentum loss during 20 ms≲  pb ≲ 150 ms for model R20B12, as ∼ 5 × 10 48 g cm 2 s −1 .The total angular momentum reaches a maximum just after bounce ( pb ∼ 30 ms) and its value is ∼ 7 × 10 48 g cm 2 s −1 , which is comparable to the estimated total angular momentum loss.These support that the spin down seen in Figure 3 is most likely caused by the angular momentum flux via the Maxwell stress.

Gravitational Wave from Matter
In this section, we investigate the GW generated by hydrodynamic motion for our jet explosion models (R20B12 and R10B12).We extract the GWs with a standard quadrupole formula (Shibata & Sekiguchi 2003;Kuroda et al. 2014).To investigate the spectral evolution of the GWs, we evaluate the viewing-angle-dependent charac-teristic strain for an optimally oriented source (Shibagaki et al. 2021).Figure 5 shows the GW waveforms, ℎ +/× and the spectrograms for the characteristic strain of each model emitted along the equatorial (left-hand column) and the polar directions (right-hand column) for a source distance of 10 kpc, respectively.
The GWs observed along equatorial direction show the burst signal at core bounce due to the rotational flattening of the core, and also show the oscillation in the ringdown phase shortly after bounce.In the first few tens milliseconds after bounce, the GW due to the prompt convection is dominated in any direction.
For the model R10B12, there is a glitch at  pb ∼ 220 ms (see the left panel).This is a numerical artifact because we switched the EOS table to the EOS table remade at this time and this switch produces a sudden change in the thermodynamic quantities that generates the glitch in the GW signal.This is because the previous EOS returned unphysical thermodynamic quantities due to the existence of multiple roots at this time, so we remade the EOS table so that the thermodynamic variables in the EOS table are uniquely determined.Although the amplitude of the glitch is large, we believe that this artifact does not drastically change the result of our simulation because the overall trends of the GW spectrogram before and after the glitch look similar.
The model R20B12 has a clear GW signature due to the jet explosion at low frequency (∼ 50 Hz) in the top left panel (equatorial observer) as discussed in the previous works (e.g., Obergaulinger et al. 2006a;Takiwaki & Kotake 2011;Jardine et al. 2022;Powell et al. 2023).This type of signal is produced by the strong prorate explosion (see also Murphy et al. 2009).The low-frequency GW signal originated from the rise of ℎ + arises at  pb ∼ 100 ms on the GW spectrogram and finally reaches ℎ + ∼ 2 × 10 −21 .The rise of the ℎ + ceases at  pb ∼ 400 ms, which roughly corresponds to the time when the PNS stops rotating and no longer be able to energize the surrounding material with its rotational and/or magnetic energy (see the bottom panel of Figure 3).On the other hand, the GW signal at such a low frequency in the model R10B12 is not as prominent as the model R20B12.This is because the less energetic jet is launched in this model (see the bottom panel of Figure 1).
To understand the origin of the GW, we investigate the contribution of different spherical shells to the GW spectrogram as in Shibagaki et al. (2020).We confirm that the long-lasting low-frequency GW signal shown in the top left panel of Figure 5 is emitted from the spherical shell with  > 200 km.The other GW features mainly come from the spherical shell with  < 200 km.Therefore, the GW with the extended spectral shape starting from  pb ∼ 400 in the model R20B12 and the GW starting at  pb ∼ 100 ms with frequency increasing from ∼100 Hz to ∼400 Hz in the model R10B12 are likely due to PNS oscillation as observed in previous works (Jardine et al. 2022;Powell et al. 2023;Bugli et al. 2023).

Gravitational Wave from Neutrino
In this section, we focus on the GW generated by anisotropic neutrino emission.Following Müller et al. (2012a), we reconstruct the angle distribution of neutrinos and compute the GW by anisotropic neutrino emission.
To see the angle dependence of neutrinos, we plot the neutrino luminosities of   , ν , and   in the polar (solid lines) and equatorial (dotted lines) observer directions for our models in the top panel of Figure 6.Their luminosities tightly correlate with their accretion history.The most slowly rotating model (R05B12) has the highest neutrino luminosities for the largest mass accretion, while the most rapidly rotating model (R20B12) has the lowest neutrino luminosities for the smallest mass accretion (see the evolution of the PNS masses in Figure 3).The bottom panel of Figure 6 shows the difference between neutrino luminosities observed along the pole ( ,p ) and the equator ( ,e ) for each model.Except for R10B13, the neutrino luminosities in the polar direction are larger than the ones in the equatorial direction in the early phase.This is most likely because the rotational flattening of the PNS makes the apparent size of the neutrino sphere for the polar observer larger.As we discuss in Section 3.2, the PNS rotation is substantially decelerated by the magnetic braking and the low mass accretion after the explosion and becomes quite small after  pb ≃ 300 ms.This change in the PNS rotation reduces the rotational flattening of the PNS, and the difference between the luminosities observed along the pole and the equator becomes quite small in the late phase.On the other hand, model R10B13 shows a small difference between neutrino luminosities observed along the pole and the equator because the very strong magnetic fields of the PNS quickly decelerate the PNS rotation right after bounce (see the bottom panel of Figure 3).This effect should be considered in the gravitational wave emission in PNS cooling phase (Fu & Yamada 2022).Next, we move on to the GW signal from anisotropic neutrino emission for our jet explosion models.This harbors the nature of memory effect in the waveform, which we would like to name as the Epstein-Chrisotodolou memory effect, by honoring the two researchers who first pointed out the GW emission from neutrinos (Epstein 1978) and later provided the mathematical evidence that the non-linearity of GR leads to the memory feature (in general, Christodoulou (1991), see collective references in Kotake (2013)).The top and middle panels of Figure 7 display the neutrino GW strains of models R20B12 and R10B12 at a source distance of 10 kpc.For both models, the neutrino R20B12 R10B12 GW signal along the equatorial direction (middle panels) is more than one order of magnitude larger than the one along the polar direction (top panels).The neutrino GW strains emitted in the equatorial direction for R20B12 finally reach ∼ 2 × 10 −20 for ℎ + and ∼ 3 × 10 −21 for ℎ × while the ones for R10B12 finally reach ∼ 6 × 10 −21 for ℎ + and ∼ 1 × 10 −21 for ℎ × .The growth of the GW amplitude emitted toward an equatorial observer stops at  pb ≃ 300 ms for both models.mThis corresponds to the time when the PNS rotation is quickly decelerated and neutrino emission becomes less anisotropic.We can confirm this interpretation by monitoring the anisotropic parameter of neutrinos,  (Müller et al. 2012a), shown in Figure 8.The degree of the neutrino anisotropy for the plus mode is strongly enhanced only in the range of 0 <  pb < 300 ms for the equatorial observer.
Compared with the GW signal of neutrinos from the non-rotating CCSN models in Vartanyan & Burrows (2020), the maximum values of ℎ + in our jet models are much higher than their values.
In order to see the evolution of the neutrino GW spectra, we plot the characteristic GW spectral amplitudes of neutrinos for the four different time windows, [0, 100], [100,300], [300,500] and [0,  end ] ms for model R20B12 (left), and [0, 110], [110,220], [220,330] and [0,  end ] ms for model R10B12 (right) in the bottom panel of Figure 7.The first three time windows correspond to the phases before, during, and after increasing the neutrino GW amplitude for the equatorial observer, respectively.In the case of model R20B12, the second time window is the dominant source of the neutrino GW spectral amplitude for the equatorial observer at low frequency.Model R10B12 also shows that the second time window is the primary source for the equatorial observer.This analysis indicates that the rapid rise of the ℎ + generates the low-frequency GW amplitude seen in the bottom panel of Figure 7.All of the three time windows fairly contribute the neutrino GW spectral amplitude for the polar observer.
We note that, since the waveform of the GW strain of neutrinos is similar to the DC component in signal processing, the neutrino component at high frequencies could be just an aliasing effect.To circumvent the problem, the choice of a well-defined time window is important, for which we have paid special attention (see Appendix A for details).Here we use equation (A4) to perform the Fourier transform.
We proceed to discuss the detectability of the total GW signals.
The top panel of Figure 9 shows the GW spectral amplitudes for models R20B12 (solid) and R10B12 (dotted) seen from the polar observer at a distance of 10 kpc relative to the sensitivity curves of the advanced LIGO (aLIGO), advanced VIRGO (AdV), and KA-GRA (Abbott et al. 2018) as well as the next-generation GW detectors of Einstein Telescope (Hild et al. 2011, ET), Cosmic Explorer (Abbott et al. 2017, CE), B-DECIGO (Yagi & Seto (2011)), and DECIGO (Nakamura et al. (2016); Isoyama et al. (2018)).In this plot, we separately draw the matter contribution (blue), the neutrino contribution (green), and the sum of them (red).The GW signal of anisotropic neutrino emission for both models mainly contributes to the total GW signal only at low frequencies (  < 10 Hz), and their matter contribution plays a dominant role at  > 10 Hz.Their GW spectral amplitude at  > 30 Hz is larger than the sensitivities of the current-generation GW detectors such as aLIGO, AdV, and KAGRA.However, the GW spectral amplitude at low frequencies (  < 30 Hz) is smaller than their sensitivities.Thus, the next-generation GW detectors such as ET, CE, B-DECIGP, and DECIGO are necessary to detect the neutrino component at low frequencies.
The bottom panel of Figure 9 is the same as the top panel but for the equatorial observer.In this case, the neutrino component of the GW spectral amplitudes becomes dominant not only at low frequencies t pb [ms] plus(pol) cross(pol) plus(eq) cross(eq) Figure 8. Anisotropic parameter of GW by anisotropic neutrino emission for models R20B12 (top panel) and R10B12 (bottom panel) seen along the pole (solid lines) and the equator (dashed lines).The plus and cross modes are depicted in red and blue, respectively.The curves for the equator are offset by -80 %.
but also at high frequencies.Therefore, the neutrino component could be detectable even by the current-generation GW detectors.

SUMMARY AND DISCUSSIONS
Having a focus on the GW signals, we have presented results from 3D GR-MHD core-collapse simulations of a 20 M ⊙ star with spectral neutrino transport.We have computed four models by parametrically changing the initial angular velocity and the strength of the magnetic fields in the core.Our results showed that the MHD outflows are produced only for models (two out of four), to which magnetic field strengths of 10 12 G and rotation rates of 1 or 2 rad s −1 are initially imposed in the core.Seen from the direction perpendicular to the rotational axis, a characteristic waveform was obtained, which exhibits a monotonic time increase in the wave amplitude.As previously identified, this stems from the propagating MHD outflows along the axis.We showed that the GW amplitude from anisotropic neutrino emission becomes more than one order-of-magnitude bigger than that from the matter contribution, whereas seen from the rotational axis, both of the two components are in the same order-of-magnitudes.Due to the Christodoulou-Epstein memory effect, the frequency of the neutrino GW from our full-fledged 3D-MHD models is in the range less than ∼ 10Hz.Toward the future GW detection for a Galactic core-collapse supernova, if driven by the MR mechanism, we point out that the planned next-generation detector as DECIGO is needed not to miss the low-frequency signals.
Finally, we shall discuss several major limitations of our work In this study, we present a sample of simulations, due to the computational expense, though extending the number of the models of Kuroda (2021).Apparently one of the most urgent tasks in the future is to perform more systematic simulations to obtain a clear understanding how the MHD dynamics is affected by (the subtle change in) the initial strength of rotation and magnetic fields in the core.As the explosion mechanism significantly depends on the strength of rotation and magnetic fields (classified as , -Ω and MR in Obergaulinger & Aloy 2020), so do the explosion energy (e.g., ∼ 10 52 erg in Obergaulinger & Aloy 2021 and ∼ 10 51 erg in other simulations).
Our 3D-GR-MHD simulations cover ∼300 ms post bounce.Recently, it is pointed out that galactic supernova neutrinos 1 , depending on the neutron star mass, can be observed over 10 s after bounce (Suwa et al. (2019), see also Wu et al. (2015); Nakazato & Suzuki (2020); Weishi Li et al. (2023) for recent work in the relevant context).Possible phase transition from hadronic to quark matter in the PNS (Fischer et al. 2018) can be imprinted in the GW signals (Zha et al. 2020;Kuroda et al. 2022).Apparently, long-term simulations are needed to follow the dynamics of 3D(-GR-MHD), blackhole/magnetar forming stellar collapse, which is expected to have a link to long-duration gamma-ray bursts in the context of collapsars (e.g.MacFadyen & Woosley 1999) and magnetars (e.g.Metzger et al. 2011).The quantitative GW- signals from such sources are yet to be clarified based on the first principle simulations, toward which, this study we believe makes a steady step forward.

Figure 1 .
Figure 1.Evolution of maximum shock radius (top) and evolution of explosion energy (bottom) for models R05B12 (red solid line), R10B12 (blue short dashed line), R10B13 (green long dashed line) and R20B12 (magenta dash-dotted line).

R20B12R10B12Figure 2 .
Figure 2. Snapshots of entropy on the  = 0 plane for models R20B12 (top) and R10B12 (bottom) around the time when both the north and south polar jets reach 300 km (left) and when the simulation ends (right).The black solid line represents an interface between the domains with different refinement levels.

Figure 3 .
Figure3.Time evolution of PNS masses (top), average surface magnetic fields (middle), and average rotation rates (bottom) for models R05B12 (red solid line), R10B12 (blue short dashed line), R10B13 (green long dashed line) and R20B12 (magenta dash-dotted line).The poloidal and toroidal magnetic fields in the middle panel are distinguished by line thickness.

Figure 6 .
Figure 6.Neutrino luminosity of   , ν and   (top, middle, and bottom small panels in top panel) seen along the pole (solid lines) and the equator (dotted lines), and difference between neutrino luminosities along the pole and the equator (bottom panel) for models R20B12 (red lines), R10B12 (blue lines), R10B13 (green lines) and R05B12 (magenta lines), sampled at a radius of 400 km.

Figure 7 .
Figure 7. Top and middle panels: neutrino GW strains of plus (red solid lines) and cross (blue dashed lines) modes observed along the pole (top panels) and along the equator (middle panels) for models R20B12 (left panel) and R10B12 (right panel) as a source distance of 10 kpc.Bottom panels: characteristic neutrino GW spectral amplitudes of models R20B12 (left panel) and R10B12 (right panel) for the four different time windows indicated in the legends of each panel.The GW sources are observed along the pole (solid lines) and along the equator (dotted lines) as a source distance of 10 kpc.