Characterizing the Directionality of Gravitational Wave Emission from Matter Motions within Core-collapse Supernovae

We analyze the directional dependence of the gravitational wave (GW) emission from 15 3D neutrino radiation hydrodynamic simulations of core-collapse supernovae (CCSNe). Using spin weighted spherical harmonics, we develop a new analytic technique to quantify the evolution of the distribution of GW emission over all angles. We construct a physics-informed toy model that can be used to approximate GW distributions for general ellipsoid-like systems, and use it to provide closed form expressions for the distribution of GWs for different CCSN phases. Using these toy models, we approximate the protoneutron star (PNS) dynamics during multiple CCSN stages and obtain similar GW distributions to simulation outputs. When considering all viewing angles, we apply this new technique to quantify the evolution of preferred directions of GW emission. For nonrotating cases, this dominant viewing angle drifts isotropically throughout the supernova, set by the dynamical timescale of the PNS. For rotating cases, during core bounce and the following tens of milliseconds, the strongest GW signal is observed along the equator. During the accretion phase, comparable—if not stronger—GW amplitudes are generated along the axis of rotation, which can be enhanced by the low T/∣W∣ instability. We show two dominant factors influencing the directionality of GW emission are the degree of initial rotation and explosion morphology. Lastly, looking forward, we note the sensitive interplay between GW detector site and supernova orientation, along with its effect on detecting individual polarization modes.


INTRODUCTION
The end of massive stellar evolution is marked by a corecollapse supernova (CCSN).These stellar explosions, and in some cases implosions, are dynamic events influenced by a variety of physical processes, from neutrino processes, to hydrodynamic turbulence, to magnetic field interactions.Moments after the supernova is launched, the birth of a compact object called the protoneutron star (PNS) occurs, which will ultimately cool to a neutron star or-with sufficient mass accretion-will collapse to a black hole.The PNS contains the neutron rich remnants of the once-stellar iron core and is characterized by nuclear matter at densities ≳ 2 × 10 14 g cm −3 near the center.As turbulent downflows interact with the PNS surface and convection within the PNS develops, the PNS can oscillate, generating gravitational waves (GWs) (Sotani & Takiwaki 2016;Kotake & Kuroda 2017).The principal source driving PNS oscillations-PNS convection (Andresen et al. 2019;Mezzacappa et al. 2023) or accretion downflows (Vartanyan et al. 2023)-remains an ongoing discussion.In the case of accreting matter, the luminosity of the GWs is proportional to the square of the rate at which turbulent kinetic energy is accreted by the PNS (Thorne 1996) (Also see Müller (2017) for order of magnitude estimates of GW features).
These oscillations can exhibit different 'modes' which correspond to unique frequencies of emission, based on the restoring force driving the PNS.Some examples include gmodes driven by gravity, p-modes driven by pressure, rmodes driven by rotation, and w-modes driven by oscilla-tions in spacetime.As each mode has a unique restoring force, GW observables can encode different characteristics describing the supernova system, such as the average density of the PNS (g-modes) or degree of rotation (r-modes) (Unno et al. 1989).
Beyond PNS oscillations, other sources of GWs within CCSNe include the 'bounce signal' from a deformed rotating PNS (Dimmelmeier et al. 2008), hydrodynamic turbulence (Pajkos et al. 2019), asymmetric emission of neutrinos (Vartanyan & Burrows 2020), compact object ejection (Burrows & Hayes 1996), and GWs from fluid instabilities like the standing accretion shock instability (SASI) (Kuroda et al. 2016;Andresen et al. 2019), or the low T/|W| instability (Shibagaki et al. 2020).While generating predicted waveforms for CCSNe can be valuable for improving detectability of the next Galactic event, connecting characteristics of the signals to internal source physics is vital.Each of these sources also encodes physical characteristics of the supernova center, a region previously unobservable with electromagnetic (EM) radiation.GWs from fluid instabilities provide information regarding the strength of convection and the timescale on which it occurs (Andresen et al. 2017;Takiwaki et al. 2016;Radice et al. 2019).Neutrino sources provide the degree of asymmetric neutrino production (and to some extent mass accretion) (Vartanyan & Burrows 2020).The bounce signal is directly related to the rotational content of the PNS just after core bounce (Dimmelmeier et al. 2008).Furthermore, GWs from CCSNe depend on the hot nuclear equation of state (EOS) (Pan et al. 2018).Several works have also proposed a joint analysis of GWs and neutrinos, to better constrain physical conditions of the PNS (Nagakura & Vartanyan 2022), SASI activity (Kuroda et al. 2017), EOS insights (Vartanyan et al. 2019), and stellar compactness along with explosion properties (Warren et al. 2020).
With the recent Supernova 2023ixf discovery from Itagaki (Teja et al. 2023) and the fourth observing run for LIGO (LIGO Scientific Collaboration et al. 2015) underwayalong with observations by Virgo (Acernese et al. 2015) and KAGRA (Kagra Collaboration et al. 2019)-supporting GW observers with theoretical predictions is paramount.With various physical insights from CCSN gravitational radiation, considering these signals in the context of observability is important to distinguish between signal features that are interesting and signal features that are detectable, therefore more valuable.One often considered factor is the GW amplitude, or the GW strain, commonly denoted h.Depending on the specifications of the GW observatory, detectors have certain sensitivities, or a minimum threshold below which GW amplitudes fall below the noise floor.Supernova models often produce time domain waveforms (TDWFs) to observe the predicted GW amplitudes over time for a given event.
The GW frequency content is equally valuable.GW detectors have limiting factors that constrain the frequency range of observable signals.For example, for the Laser Interferometer Gravitational-wave Observatory (LIGO), shot noise imparted by the individual laser photons on the test mass limit the upper limit of the frequency range ∼ 8 kHz.By contrast, seismic noise establishes the lower limit of the frequency range ∼ 10 Hz (Aasi et al. 2013).One common tool for determining the frequency content of the CCSN GW signal includes spectrograms, which display the frequency content of the signal as time evolves.Characteristic strain plots are another method that displays the frequency distribution of the signal summed over a given time interval.Characteristic strain plots are often used to compare cumulative GW frequency content to GW detector sensitivity curves (Moore et al. 2014).
Directly measuring the polarization of the GW signal is an ongoing area of research.According to Einstein's theory of general relativity (GR), GWs exhibit two polarization states: plus (h + ) and cross (h × ) modes.These states can serve as a set of bases, a linear combination of which can describe any GW signal; as a parallel with EM radiation, GWs too can be circularly or elliptically polarized (e.g., Shibagaki et al. 2021).In principle, arrays of GW detectors at different orientations can be used to detect GW polarization.However, current detector sensitivities have difficulty discerning GW polarization even from sources with stronger amplitudes than CCSNe, such as black hole mergers (Gair et al. 2013).
One tool that has emerged to predict the polarization from numerical models includes 'polograms', which describe the relative strength between the h + and h × modes.Hayama et al. (2016) connect the circular polarization of GWs with rotation near the center of the supernova.Hayama et al. (2018) conclude the circular polarization from nonrotating CCSNe can encode information regarding the SASI and PNS g-mode oscillations.Chan & Hayama (2021) quantify the degree of circular polarization through quantities such as Stoke's parameters, within existing GW detection pipelines, such as Coherent WaveBurst (cWB) (Klimenko et al. 2021).
Before the advent of modern high performance computing platforms, simplified analytic models were used to estimate GW production from compact objects.One example is treating rotating compact objects as rotating, triaxial ellipsoids (Chin 1965;Chau 1967;Saenz & Shapiro 1978;Zimmermann & Szedenits 1979;Gal'tsov & Tsvetkov 1984;Thorne 1996;de Freitas Pacheco 2010), with deformations from magnetic fields (Bonazzola & Gourgoulhon 1996;Konno et al. 2000;Palomba 2001) and even neutron star (NS) crust interaction (Ushomirsky et al. 2002).This simplified description has also been used to describe so called 'bar-mode' instabilities proposed to form in CCSN centers with sufficient rotation (New et al. 2000).Beyond simple physical expres-sions, with increased access to modern computing resources, more recent multidimensional simulations have noted the impact of viewing angle on expected GW amplitudes.
In the context of the low T/|W| instability, Ott et al. (2005) note a stronger GW strain along the axis of rotation for short term purely hydrodynamic models without radiation transport.Scheidegger et al. (2008) and Scheidegger et al. (2010) find similar results for rotating, magnetized models.Similarly, Kuroda et al. (2014); Takiwaki & Kotake (2018); Powell & Müller (2020); Shibagaki et al. (2021) corroborate these results when comparing polar and equatorial GW strains.Kotake et al. (2011) consider the role of rotation influencing neutrino contributions to GW signals.While providing valuable first steps towards understanding the directionality of GWs from CCSNe, these works only consider viewing angles along the equator and axis of rotation.Furthermore, certain observationally-focused works consider different CCSN source orientations (Abbott et al. 2020;Szczepańczyk et al. 2021Szczepańczyk et al. , 2023)).Even with these considerations, there stands open questions: are there viewing angles-not necessarily aligned with the cardinal axes-along which GW amplitudes are strongest, and, if so, do they evolve in time?Furthermore, what are the physical processes that dictate these preferred viewing angles?
Due to cheaper computational cost, the use of axisymmetric simulations has the advantage of pushing to later times for a wider variety of initial conditions, but can only compute h + , due to the geometry of the domain.Axisymmetric models have also been shown to develop large asymmetries along the symmetry axis, leading to GW amplitudes larger by an order of magnitude, compared to 3D works (O'Connor & Couch 2018).Nevertheless, extensive work has been completed leveraging these 2D works in the context of exploding CCSN scenarios (Murphy et al. 2009;Yakunin et al. 2010;Müller et al. 2013;Yakunin et al. 2015;Morozova et al. 2018).Marek et al. (2009), Pan et al. (2018), andEggenberger Andersen et al. (2021) explore EOS-dependent features of CCSN evolution and the corresponding imprint on the GW signal.Pajkos et al. (2021) use multiple components of the GW signal to constrain the progenitor compactness.
Recently, Mezzacappa et al. (2020) conducted a 3D study outlining the different sources of GW emission from matter motion within a 15 M ⊙ star.Raynaud et al. (2022) use GWs to constrain convective motion within magnetized PNSs.Mezzacappa et al. (2023) conduct a numerical experiment with three 3D models, constraining the physical engine generating GWs at different frequencies.Afle et al. (2023) use multidimensional models to analyze the detectability of f-mode GW frequencies with third generation GW detectors.More exotically, Urrutia et al. (2023) explore the impact of viewing angle of GW emission from long GRB jets.
At later times in the supernova evolution, Mueller & Janka (1997) note the large imprint of neutrino asymmetries on the expected CCSN GW signal; Vartanyan & Burrows (2020) and Vartanyan et al. (2023) investigate the directional dependence of GW emission generated directly from asymmetric neutrino emission.Richardson et al. (2022) investigate the directional dependence of GW amplitudes from neutrinos and matter, for a single model.Mukhopadhyay et al. (2022) investigate the observability of these low frequency signals for the DECi-hertz Interferometer Gravitational Wave Observatory (DECIGO) (Kawamura et al. 2011).Likewise, work continues to detect CCSN GWs with next generation GW detectors (Srivastava et al. 2019), with some focusing on decihertz frequencies (Arca Sedda et al. 2020).Various works have also investigated the expected contributions to the stochastic GW background from CCSNe (Buonanno et al. 2005;Finkel et al. 2022).
The main goal of this work is to characterize the time evolution of the directional dependence of GW emission for both rotating and nonrotating nonmagnetized CCSNe.To fill the need to understand how the preferred direction of GW observations evolves in time for a variety of initial conditions, we present a new method to understanding the directionality of GW emission from CCSNe for arbitrary viewing angles.We determine an evolving preferred direction of GW emission during CCSNe, for both rotating and nonrotating cases.We observe possible correlations with the presence of the low T/|W| instability.We use simplified analytic estimates to show the distribution of GW amplitudes over viewing angle, depending on the supernova phase.Identifying preferred viewing angles of dominant GW emission would help inform future population studies of GW sources to the GW background.Likewise, understanding the emergence of a preferred direction of GW emission-if any-during specific phases of the supernova evolution (e.g., bounce, prompt convection, accretion phase, explosion, etc.) would connect CCSN source physics to observability results seen in GW CCSN observation papers (e.g., Szczepańczyk et al. 2023).For example, two relevant questions we answer: what are the CCSN physical conditions that leave detectability relatively agnostic to source orientation?For each phase of CC-SNe, what are progenitor characteristics that allow for a preferred direction of GW emission, which may in turn affect detectability?
This paper is organized as follows: in Section 2 we review the methods used to set up our simulations and perform our analysis.In Section 3.1 we review current techniques used to analyze GW signals.Section 3.2 introduces a new visualization method to investigate the directional dependence of the GW emission at a given instance in time.Section 4 uses simplified analytic arguments to obtain closed form solutions for the GW direction for different phases of the supernova.
Section 5 visually uses this method to track the time evolution of the directionality of GWs.Section 6 quantifies these observations and connects their source physics.Section 7 quantifies the evolution of the directionality when accounting for both polarizations.Section 8.1 discusses the physical implications generating the directional dependence.Section 8.2 discusses future benefits for observability.Section 8.3 provides a worked example discussing the influence of directionality on a GW bounce signal detection.Section 8.4 discusses directionality considerations in the context of signal injection into data analysis pipelines.Finally, in Section 8.5, we summarize.

Numerical Models
This work examines 15 3D neutrino radiation hydrodynamic simulations of CCSNe.The first set of five models are referred to as the mesa set, which are nonrotating 20 M ⊙ models used in O' Connor & Couch (2018).They make use of the SFHo EOS (Hempel et al. 2012;Steiner et al. 2013), M1 neutrino transport (O'Connor & Couch 2018), and are evolved for ∼ 0.5 sec post-bounce (pb).Five additional simulations, one nonrotating and four rotating, are analyzed from 40 M ⊙ models from Pan et al. (2021) and 20 M ⊙ models.These use the LS220 EOS (Lattimer & Swesty 1991).To account for neutrino interactions, they use parameterized deleptonization (Liebendörfer 2005) through collapse, with the isotropic diffusion source approximation (IDSA) for neutrino transport (Liebendörfer et al. 2009) after bounce; we abbreviate this combination as Ye(ρ)+IDSA.These 10 models are completed with the FLASH multiphysics code (Dubey et al. 2009;Fryxell et al. 2010).Likewise, they use the Newtonian multipole solver from Couch et al. (2013), supplemented by the general relativistic effective potential (GREP) proposed by Marek et al. (2006).We also make use of five models from Burrows et al. (2019), Radice et al. (2019)  1 ,  and Vartanyan et al. (2023) 2 .These simulations make use of the Fornax code (Skinner et al. 2019) and have comparable microphysics to FLASH.For neutrino interactions, Fornax accounts for inelastic scattering off electrons and nucleons.The neutrino physics used in the version of FLASH from the mesa suite does not account for inelastic processes (O'Connor & Couch 2018;O'Connor & Couch 2018), which may contribute to why the models from the mesa suite fail to explode, whereas many models from Vartanyan et al. (2023) successfully explode.The specific initial conditions are outlined in the aforementioned works.One advantage of using models from FLASH and Fornax is they employ dras-tically different grids: FLASH uses the Paramesh (MacNeice et al. 2000) adaptive mesh refinement library on a Cartesian grid, and Fornax uses a static, spherical, dendritic grid.This diversity of the dataset also allows us to consider multiple zero-age main-sequence masses, rotation rates, EOSs, grid perturbations at collapse, neutrino treatments, and GWs from multiple supernova phases through explosion.For a concise review of the relevant input parameters, see Table 1.As a point of emphasis, the GW analyses in this work examine GW sources that arise from matter contributions only.For insightful works regarding GW emission from neutrino asymmetries, see Vartanyan & Burrows (2020), Richardson et al. (2022), andVartanyan et al. (2023).

Gravitational Wave Analysis
As the gravitational treatments used in these models do not formally evolve a spacetime metric, the GW generation must be calculated during a post processing step.We make use of the generic quadrupole formulae presented in Oohara et al. (1997) and where Qij represents the second time derivative of the reduced mass quadrupole moment in an orthonormal basis, D is the distance to the source, G is Newton's gravitational constant, and c is the speed of light.For all calculations of h + and h × , we assume a fiducial distance of D = 10 kpc.Expanding the angular Q values in terms of Cartesian components yields In practice, simulation data tracks each of the Qij components through numerical integration (Finn & Evans 1990), and a finite difference in time is performed to construct Qij .These values are then applied to all altitudinal angles of 0 < θ ≤ π and azimuthal angles 0 < ϕ ≤ 2π to construct the surface plots outlined in Section 3.2.As a note, θ is measured with respect to the positive z axis, and ϕ is measured with respect to the positive x axis.The rotating models in this work align the rotation axis with the positive z axis.

Spin Weighted Spherical Harmonic Decomposition
A key part of our analysis lies in quantifying the distribution of GW strain at each point in time.To perform this decomposition, similar to Gossan et al. (2016) and Ajith et al. (2007), we calculate coefficients for spin -2 weighted spherical harmonics characterized by coefficients where −2 Y * lm (θ, ϕ) is the complex conjugate of the spin -2 weighted spherical harmonics, D is the source distance (to cancel out the factor of 1/D in Equation (1) and Equation (2)), and i is the imaginary number.In the plots below, we select the real component of −2 a lm .The use of spin weighted spherical harmonics is advantageous because it provides an orthogonal set of basis functions to describe multidimensional surfaces and has been used to great success describing GW distributions in the numerical relativity community (Boyle et al. 2014).Likewise, for varying values of m, they conveniently describe differing altitudinal angles; for example, m = 0 peaks at θ = π/2 (the equator of a sphere), whereas m = ±2 peaks at θ = 0, π (the pole of the sphere), a characteristic that will prove useful for this work.The expressions for −2 Y lm (θ, ϕ) can be found in Table 2 (Ajith et al. 2007).

VISUALIZING GRAVITATIONAL WAVE EMISSION
In this section, we begin with existing methods quantifying GW data, and sequentially generalize these techniques to eventually consider GW directionality.All visualizations in this section make use of models s40o[0, 0.5, 1].

Time Domain Waveforms
We begin with the TDWF from s40o0, s40o0.5, and s40o1.Figure 1 displays h + throughout the simulation duration-note the different y scales between each case.The left column corresponds to the expected GW signal for an observer along the line of site of the CCSN equator (θ = π/2) and along the positive x axis (ϕ = 0).The right column corresponds to an observer along the supernova north pole (θ = 0), which represents the axis of rotation for rotating models.The top row represents the nonrotating model s40o0, the central row represents the slow rotating model s40o0.5, and the bottom row represents the fast rotating model s40o1.As rotation increases, coherent motion of the fluid generates stronger oscillations of the PNS.Likewise, in the fast rotating model, beginning ∼ 150 ms pb, the presence of the low T/|W| instability excites the PNS to generate GWs with larger amplitudes as well (Pan et al. 2018).
When comparing signals of a given rotation rate, but differing viewing angle, there are only moderate differences for s40o0.This result is expected because there is no centrifugal support to deform the PNS along a particular plane (e.g., a more oblate PNS in the xy plane for an axis of rotation in the z direction).However, as rotation increases, the GW signal detected by observers at different viewing angles becomes more pronounced.For the rotating cases, the GW signal just after bounce, t pb = 0 sec, commonly called the bounce signal, has a large GW amplitude when viewed along the equator, but has an amplitude of 0 when viewed along the axis of rotation.This effect is explained through Equation (1) and Equation (2).When viewed along the equator, the oblate PNS will have a cross section similar to an ellipse.As the infalling material deforms the PNS at the time of bounce, the cross section changes to a more spherical shape, and Qxz , Qyz , and Qzz will take on correspondingly larger valuesa higher amplitude bounce signal is generated.By contrast, when observed along the axis of rotation, the PNS cross section is nearly circular.As the infalling material does not have an azimuthal dependence just after bounce, the cross section of the PNS retains its circular shape.The dominant Qxx , Qxy , and Qyy retain smaller values, thereby preventing an observed GW signal.Hundreds of milliseconds after bounce, the GW amplitude related to the PNS oscillations also displays variations in amplitude, depending on viewing angle.
Examining TDWFs can be valuable in investigations such as this one.However, they are inherently limited to observing one viewing angle at a time.In general, there is no guarantee that the maximum GW signal must be along the equator or pole.Furthermore, as the mechanism of GW generation can change throughout the supernova evolution, in principle, a viewing angle of preferred GW emission may evolve through time.

Visualizing GWs in Multiple Dimensions
To address the need to understand GW emission along every viewing angle, we introduce the strain surface plot.This figure is inspired from those in Vartanyan & Burrows (2020) that display the GW emission by color and concavity on a 3D surface.It can be thought of as an extension to 3D of Figure 3.7 of Maggiore (2007) and applied to strain.In practice, it can be difficult to compare color or concavity by eye, so we offer another form of visualization, with examples provided in Figure 2a, Figure 3a, Figure 4a, and Figure 6a.
These surfaces are taken from different points in time for s40o1.To interpret these plots, consider Figure 2a at t pb ∼ 2 ms.The distance from the origin to any point on the green surface corresponds to h + along that viewing angle.The purple star indicates the direction along which there is a maximum h + .The x and y axes form the plane of the supernova equator, and the z axis is the axis of rotation.Figure 2a displays similar h + emission when viewed at any angle in line with the equator.By contrast, when viewed along the axis of rotation, the surface does not extend far from the origin, indicating a weak GW amplitude.Physically, this behavior is justified.As explained in Section 3.1 and shown in the top row of TDWFs in Figure 1, the bounce signal is detected when observed along the equator due to the geometry of the deformed PNS; a more detailed explanation of the PNS dynamics is reserved for Section 4. Figure 2a offers a compact way to conclude this directional dependence at a given point in time without relying on multiple TDWFs.
In Figure 3a, we notice a deviation from azimuthal symmetry, around 13 ms post-bounce.At this point, a combination of the ringdown from the bounce and emergence of prompt convection begins to skew the symmetry of the emission.The direction of maximum h + amplitude is roughly aligned with the x axis at this time.Furthermore, a moderate GW amplitude can be observed along the z axis, in contrast to the bounce phase.
Figure 4a displays GW amplitudes nearly 270 ms later.The configuration of the h + surface preferentially lies along the axis of rotation, with a distinct four lobed structure.There are also non-negligible deformations of the PNS along all three spatial degrees of freedom, likely due to the turbulent nature of the accreting matter.During the accretion phase, for models s40o.5 and s40o1, this four-lobed structure varies in amplitude, but retains similar morphology.Figure 6a displays the strain surface plot for s40o0 roughly 230 ms pb.We see the emergence of nonaxisymmetric behavior, paired with GW maxima misaligned from all coordinate axes.This is due to accretion downflows striking the PNS along a direction misaligned from the principal axes.
For convenience, we offer methodology to translate from strain surface plots to TDWFs: 1. Construct a strain surface plot at a time of interest t, 2. identify a desired line of site (θ, ϕ).In this example we pick a viewing angle in the equator: the x axis (π/2, 0), 3. the distance from the origin to the green surface along this line of site represents the injected GW strain h(θ, ϕ, t) for a detector, 4. the values (t, h) correspond to the ordered pair on the TDWF.
As a point of emphasis, these surface plots represent GW strain detected by an observer far from the GW source.Equation (1) and Equation (2) are applicable only for distant observers from slow moving (v ≪ c) sources.These do not represent the GW amplitudes generated just inside the supernova.To calculate these quantities, a more robust treatment of relativity is needed, with numerical models that track quantities of the spacetime metric.

ANALYTIC ESTIMATES WITH TOY MODELS
Armed with established GW strain surfaces from Section 3.2, we now move to describing their origin by means of a simplified example.In this section, we step through the different phases of GWs from CCSNe: bounce, ringdown, and accretion phase.We explain the directionality of GWs through the use of an intuitive example of a wobbling, oscillating ellipsoid toy model to reproduce h + along a given direction.Considering both polarizations will be addressed in Section 7.

Introducing the Wobbling, Oscillating Ellipsoid
To quantify GW amplitudes with closed form solutions, we approximate the PNS as a uniform density, solid-body ellipsoid with three principal axes a, b, and c for the x, y, and z axes, respectively.We acknowledge that PNSs retain complex convective structures and undergo differential rotation that can cause them to deviate from solid-body rotation; for more detail, see Appendix B. To more realistically account for the complex dynamics-while retaining closed form solutions-we allow the ellipsoid to evolve in three ways: rotation, tilt, and size changes.To describe how each physical component of the ellipsoid contributes to GW generation, consider a simplified mathematical example.As h ∝ Q, a simplified expression of the reduced mass quadrupole moment of an ellipsoid is Q ∼ M R 2 , for a given mass M , and radius R. Applying two time derivatives yields a useful reference estimate Here the time derivatives represent different physical mechanisms that contribute to the GW generation: Ṁ the mass ac- /By] (c) Entropy slice through the pole (z axis) from model s40o1 2 ms after bounce.Notice the PNS surface (solid black line: 10 11 g cm −3 ) slightly prolate, after being deformed from the infalling matter.The PNS will ringdown, finally achieving an equilibrium as an oblate spheroid, due to centrifugal effects.For scale, every cm on the page corresponds to 5.6 × 10 9 cm s −1 .cretion rate, M the rate of change of the mass accretion rate, Ṙ the velocity of the PNS surface, and R the acceleration of the PNS surface.When considering these order of magnitude estimates, it is important to remember contributions of these terms-for example, velocities Ṙ-must contribute to nonspherical mass motions.As this work is primarily concerned with the dynamics of the PNS, we pay special attention to the fixed mass case, when M = Ṁ = 0, and constant speed case R = 0.The remaining quantities-M , R, and Ṙ-are free parameters.In the upcoming sections, we will assign values to these parameters to mimic conditions in our numerical simulations.The values of M are taken from our simulation outputs of the PNS baryonic mass at a given time.The values of R along each direction correspond to a, b, and c values; they correspond to the simulation values of the PNS radius along the x, y, and z axes, respectively.Values of velocities near the PNS surface, Ṙ, along each direction correspond to ȧ, ḃ, and ċ.For intuitive explanations, we also refer to a mass quadrupole moment matrix of the form In Equation A2we explicitly define this quadrupole moment matrix for our toy model case, which will pick up offdiagonal components when considering a tilted, rotating ellipsoid.To connect the reduced mass quadrupole moment Q used in Equation (1) and Equation ( 2), with the mass quadrupole moment I used in Appendix A, one can perform This toy model will be used to draw connections between PNS dynamics during each phase of the supernova and the associated distribution of GWs.Furthermore, it can be used to build an intuition to identify under which conditions different dynamics may dominate the GW signal, for example, surface oscillations versus nonaxisymmetric rotation.For future astrophysics studies that explore GW sources approximated as rotating, dynamic, tilted, ellipsoids, this model can be used to estimate the directionality of GW emission, beyond simple first order estimates that rely on tilt angles θ tilt ≪ 1.

The Bounce Phase
Occurring at bounce, before significant hydrodynamic instabilities manifest, the centrifugally deformed PNS is nearly axisymmetric.With a nearly spherical infall of matter, the ram pressure imparted on the PNS deforms it, making it less oblate-more spherical.As the PNS rings down, the PNS will oscillate between a prolate and oblate configuration over the timescale of milliseconds, losing energy by exerting pressure on the overlying accreting fluid (or PdV work), and eventually settle in an oblate configuration, depending on its degree of rotation.When viewed at a viewing angle perpendicular to the axis of rotation (along the equatorial plane), the PNS cross section appears as an ellipse with varying principal axes.At this viewing angle, the mass quadrupole moment is accelerating.When viewed perpendicular to the equator (along the axis of rotation), the PNS cross section resembles a circle with varying radius.However, the ratio between the principal axes remains nearly constant; no GWs are generated along the axis of rotation.Although the PNS remains nearly axisymmetric during this bounce and subsequent ringdown, the dynamics of the ringdown are what generate the GW burst mostly along the equator-the Ṙ and R terms in Equation ( 7).Hence, a clearer picture emerges: with axisymmetric bounce dynamics, a resulting axisymmetric GW distribution is generated for this set of coordinate axes (Kotake & Kuroda 2017).For clarity, refer to Figure 2c.This image is a slice along the axis of rotation of model s40o1 2 ms after bounce.Color corresponds to entropy, and arrows denote the velocity of the fluid.The solid black line is the PNS surface chosen where the density ρ = 10 11 g cm −3 .This visual shows a prolate spheroid with equatorial radius ∼ 60 km and polar radius ∼ 62 km. Figure 2d is taken at the same time but shows a slice through the equator, displaying the axisymmetry of the matter profile.It is this axisymmetric matter configuration, with evolving principal axes, that motivates our toy ellipsoid model of M ellipsoid = 1M ⊙ , with polar radius c = 62 km and equatorial radii a = b = 60 km.We select a velocity at the poles of ċ = 600 km s −1 , taken from the simulation output.Figure 2b displays the GW output with these toy parameters.The correspondence between the distribution of GW generation agrees quite well with Figure 2a, displaying an axisymmetric GW bounce signal, with maximum amplitude along the equator and no GW amplitude along the pole.Referring to Equation ( 7) it is the 2M Ṙ2 term generating this GW distribution.However, we note that magnitude of h + differs by nearly 2 orders of magnitude between the simulation data and toy model.While the 2M RR term may partially contribute and additional model tuning could achieve closer numbers, our goal is not to exactly reproduce GW amplitudes with these toy models.Our goal of this section is to give instructive examples connecting the directions of the GW strain to the internal matter dynamics.

The Ringdown Phase
We define the ringdown as the tens of ms following bounce, in which the rotating PNS finds a new equilibrium after being perturbed from redirecting initially infalling material to an outgoing shock.For tens of ms pb, as the PNS achieves a new equilibrium, nonaxisymmetric structure may arise as the deleptonization due to the neutrino burst leaves a negative lepton gradient.In the case when rotation is not  dominant enough to stabilize the PNS system, according to the Solberg-Hoiland stability criterion (Endal & Sofia 1978), prompt convection can occur (Mazurek 1982).Once again, we refer to simulation outputs from s40o1 to understand the physical picture.Figure 3c displays a nearly circular cross section along the axis of rotation.Compared to the bounce phase, the PNS is transitioning towards a more stable configuration that is less prolate compared to Figure 2c, as expected.Small perturbations due to prompt convection allow for slight perturbations in the y-z plane as well.In Figure 3d, the nonaxisymmetry becomes more clear.In particular, prompt convection creates 4 maxima, causing deviations from the circular surface during the bounce phase.While not a true ellipsoid, the degree of nonaxisymmetry still motivates our numerical choices for our toy model with M = 1.3M ⊙ , with equatorial radii a = 83 km and b = 81 km.The polar radius is c = 77 km.By asserting the majority of the motion will occur along the equator, we choose ȧ = ḃ = 5000 km s −1 and ċ = 0.For the angular velocity, we observe a value of Ω = 100 rad s −1 from the simulation output.As a caveat, we acknowledge 5000 km s −1 is a large value for an equatorial velocity tens of ms after bounce.In our simulation output we only note typical values ∼ 100 km s −1 .However, because this toy model is solid-body, the Q term will be more strongly influenced by the rotational motion about its axis of rotation, compared to a less solid-body, hydrodynamically evolving PNS.Hence, a higher ȧ and ḃ are required to more strongly contribute to the Ṙ component of Q ∼ 2M Ṙ2 of the toy model.The exact angular velocity profile during the ringdown phase is shown in Appendix B.
Figure 3b shows the GW distribution from the toy model.The combination of rotation, paired with a nonaxisymmetric object, is what drives the second characteristic distribution of GWs during the ringdown phase.Figure 3a displays this distribution characterized by a maximum along the equator, similar to the bounce phase.However, within the equator, a key difference arises π/2 radians from either GW maximum, showing h + amplitudes that are smaller by a factor of 2. Additionally, nonnegligible h + along the pole arises.Thus there is a superposition of two dynamics.
The first form of dynamics is the remaining ringdown that will complete over the timescale of tens of ms, causing the prominent GW amplitudes along the CCSN equator, as described in the previous section.The second is the effect of a rotating, nonaxisymmetric body.Consider the work of Creighton & Anderson (2011) (see their Section 3.4), who outline the GW emission for a solid-body, fixed-size, rotating ellipsoid, and for an angular velocity Ω, distance D, altitudinal angle θ, and azimuthal angle ϕ.In this example, the terms generating the GWs are due to the rotation of the PNS.At a given point in time, note h + varies with altitudinal angle θ as 1 + cos 2 θ.This implies expected maxima closest to the axis of rotation, for a solid-body rotator.Secondly, varying with azimuthal angle ϕ, the cos 2ϕ term indicates 4 GW maxima equally spaced by π/2 radians.In Figure 3b, the 1 + cos 2 θ term contributes to the nonzero GW amplitudes along the axis of rotation, arising from the nonzero difference I xx − I yy , mathematically quantifying the deviation from axisymmetry.
The modulated GW amplitude when viewed along the CCSN equator is characterized by the cos 2ϕ term in Equation ( 10).
The reason h + is no longer axisymmetric, is due to the mismatch between CCSN equatorial radii a ̸ = b generating a nonzero I xx − I yy .For the closed form accounting for both ringdown and nonaxisymmetric rotation, we appeal to our closed form solution from Appendix A, in particular Equation (A33).When selecting the noted toy model parameters, the strain surface plot for this rotating, oscillating ellipsoid-Figure 3b-shows good agreement with the distribution of h + compared to the simulation output: Figure 3a.

The Accretion Phase with Strong Rotation
As ringdown subsides and continued accretion onto PNS ensues, the PNS will be deformed from matter downflows onto its surface and-depending on the degree of rotationpossibly from PNS convection.Figure 4c shows the entropy profile ∼ 285 ms pb for model s40o1.With continued angular momentum accretion, the PNS has become extremely oblate with polar radius ∼ 40 km and equatorial radii ∼ 60 km.The velocity profile of the overlying material also shows the deformation of the PNS in the upper right quadrant due to downflows along the axis of rotation.The equatorial slice in Figure 4d shows slight deviations from axisymmetry as the PNS attempts to attain equilibrium after deformations from turbulent downflows.To gain a clearer understanding regarding the multidimensionality of the flow, we refer to Figure 5.This view into the inner 100s of km of the CCSN denotes higher values of the magnitude of the total spatial velocity in brighter hues.The most inner contour appears near the PNS surface, illustrating how oblate the compact object has become from ample angular momentum accretion.Just above the PNS surface, even higher velocity material is present.As the PNS surface is nearly concentric with the brighter, higher velocity material just above, one is provided with another depiction of the angular momentum accretion generating deviations of the PNS from spherical symmetry.This oblateness paired rotating nonaxisymmetry motivates our choice of toy Figure 4: Strain surface plots and entropy slices for model s40o1 during the accretion phase.The x and y axes form the equator of the supernova; the z axis indicates the axis of rotation.For the top two panels, the distance from the origin to a point on the green surface represents the detected h + along that direction.The bottom two panels are entropy slices, with overlaid velocity profiles (arrows).A clearer 3D picture of the dynamics can be seen in Figure 5.  4c and Figure 4d motivates our approximation of ȧ = ḃ = ċ = 0.In reality, the PNS will be dynamic; however, this assumption implies that rotation is dominant, or that the rotation of the nonaxisymmetric features (rather than their radial motion) contribute more importantly to the 2M Ṙ2 term in Equation ( 7).
Figure 4b shows the GW distribution for the solid rotating ellipsoid, displaying a roughly four lobed 'clover-like' structure, with relative GW maxima along the axis of rotation.To explain the behavior in the strain surface plots, we appeal to Equation (10).Entering as 1 + cos 2 θ, for θ = 0, π this term exhibits a maximum, consistent with the axis of rotation.To explain the four lobed structure, observe the cos 2ϕ term.Similar to the explanation during the ringdown phase, as an observer circumscribes the CCSN in the azimuthal di-rection, the cos 2ϕ term will exhibit two sets of maxima and minima: one occurring at ϕ = 0, π and π/2, 3π/2.These two sets are reflected by the orthogonal lobes depicted in Fig- ure 4b.However, the key difference between the ringdown and accretion phase is that rotation dominates the dynamics.Thus, the GW amplitude at ϕ = 0, π/2, π, and 3π/2 are of similar amplitude.The behavior of h × displays similar maxima along the axis of rotation and 4 GW maxima in the azimuthal direction.However, h × is out of phase by π/4 radians, compared to h + .When comparing the toy model of Figure 4b with the simulation output of Figure 4a, we notice good agreement in the overall GW distributions.Once again, the amplitudes differ by two orders of magnitude, indicating non-solid-body effects becoming more dominant.Similar to before, the solid-body rotation in the toy model contributes more matter rotating, leaving a larger Q ∼ 2M Ṙ2 , compared to a more differentially rotating PNS.As a point of empha- Figure 6: Strain surface plots and entropy slices for model s40o0 during the accretion phase.The x and y axes form the equator of the supernova; the z axis indicates the pole.For the top two panels, the distance from the origin to a point on the green surface represents the detected h + along that direction.The bottom two panels are entropy slices, with overlaid velocity profiles (arrows).sis, we acknowledge rotating PNSs do not exactly mimic a solid, rotating ellipsoid.We only use this simplified model to describe nonaxisymmetries, which revolve around a single axis, in a closed expression.Detailed illustrations of the angular velocity profile can be found in Appendix B.

The Accretion Phase Without (or with Weak) Rotation
The last degree of freedom we explore is due to the tilt of an ellipsoid.In the nonrotating or slowly rotating case, during the accretion phase, the PNS will still be perturbed by turbulent downflows and possibly more influenced by PNS convection.However, without centrifugal support to maintain a mass quadrupole moment, the ȧ, ḃ, and ċ terms will directly contribute to the 2M Ṙ2 term from Equation ( 7) and there will be no (less) rotational contribution from Ω for the non (slowly) rotating case.Observing Figure 6c, in this cross section, we note a deviation of the shock from circular symmetry and a nearly circular PNS cross section.Similarly, when viewing perpendicular to the equatorial plane, Figure 6d shows deviations along the azimuthal direction as well.It is this nonspherical postshock region that motivates our choice of an impulse to the PNS that is not aligned with the coordinate axes of the grid.
The physical picture is modeled with a nonrotating, spherical PNS.We assume a strong downflow perturbs the PNS at an angle θ = π/4 from the z axis.The ellipsoid principal axis is now aligned with θ = π/4 and receives an impulse causing c to evolve with ċ = 50 km s −1 .Figure 6b shows a modified picture, compared to before.The influence of tilt allows GW maxima to occur along directions misaligned from the x, y, and z axes.This behavior is contained in the analysis in Appendix A. Originating from rotating the initial matrix I ellipsoid , there will be nondiagonal components introduced for an observer outside the corotating ellipsoid frame.By introducing these off diagonal components, the expressions for h + and h × become less concise, and introduce a quadrupole matrix of the form Equation (A7).For stochastic accretion during the accretion phase, we expect this preferred angle to vary in time, depending on the direction from which the accretion flow originates.In Figure 6a, we notice similar behavior, with a direction of maximum GW amplitude that does not occur along the x, y, or z axes.

A Note on Similar Dynamics with Different Orientations
Notice the dynamics of the bounce toy model from Section 4.2 and those from the accretion phase in Section 4.5 are similar: they both involve a velocity ċ along the principal axis.The main difference is the nonrotating accretion case has the principal axis along c tilted by π/4 radians with respect to the z axis, whereas the principal axis along c is aligned with the z axis in the bounce case.Intuitively, one would expect a similar strain surface plot, simply tilted by π/4 radians.However, it is clear the strain surface plots in Figure 2b and Figure 6b are different.This difference arises because the relative contributions to the quadrupole moments have changed.These differences in Ï′TT ij propagate through to modify the angular dependence of h + and h × in Equations (A33-A34).To illustrate this point, see Figure 7.The three panels show the distributions of GW amplitudes for the tilted ellipsoid case from Section 4.5.From left to right, the panels show h + , h × , and h 2 + + h 2 × .The right panel displays the expected behavior of the GW distribution.In particular, since the motion of the ellipsoid is occurring along the axis π/4 radians from the z axis, the GW axisymmetry appears about this axis, as expected for transverse waves.As this tilt amounts to a change in coordinates, the contributions to the distribution of each polarization change (observer dependent), yet the distribution of the total physical signal h 2 + + h 2 × must remain unchanged (physics dependent).Thus, even for well-templated bounce signals, which arise from nearly axisymmetric dynamics, depending on the orientation, GW detectors may pick up plus and cross polarization modes.
A simple yet important subtlety arises for the future of GW astronomy with multiple GW detectors at different orientations.Consider uniquely oriented sets of observatories, who could detect both polarization modes-we acknowledge with current age detectors this is observationally challenging.For a given event, two sets of detectors with two coordinates systems will detect differing amounts of signal in each polarization.In particular, the key feature influencing the difference in observed GW polarizations for two given observatories is their relative orientation, rather than their relative position on Earth.Thus, when discussing the polarization content of a given event, the relevant physical quantity each will report, to make an accurate comparison, should be the total content of the signal h 2 + + h 2 × .As a point of emphasis for future CCSN GW theory works, in an attempt to further support a GW era where both polarizations can be detected, we recommend reporting information regarding both polarization modes; only reporting a single polarization mode would not give an orientation invariant quantity.

FINDING PREFERRED DIRECTIONS OVER TIME
Strain surface plots can be powerful diagnostic tools for analyzing the directionality of GW emission at a specific instance in time.However, CCSNe are dynamic systems whose GW generation evolves over seconds.To track the evolution of the preferred direction emitting gravitational radiation, refer to Figure 8.Each point in this 3D 'strain surface scatter plot'-henceforth referred to as an S 3 plot-refers to the direction along which a maximum h + would be measured; that + + h 2 × are shown from left to right.While similar to the dynamics of the bounce phase, we note different contributions to each polarization.Whereas a coordinate system aligned with the principal ellipsoid axis along c will only produce h + emission (Figure 2b), a differently oriented GW source will produce GW amplitudes with differing contributions in each polarization.When considering h 2 + + h 2 × on the right, the axisymmetric GW amplitudes are recovered along the axis tilted by π/4 from the z axis.In short, two systems with similar dynamics will produce similar angular distributions of total amplitudes h 2 + + h 2 × , but differing amounts in each polarization.
is, the purple stars (like those seen in Figure 2a, Figure 3a, Figure 4a, and Figure 6a) from all time steps are recorded from the entire simulation duration and plotted.Points that have brighter colors occur later in the supernova evolution.
The nonrotating model, s40o0, is represented in Figure 8a.As dark points (early times) are close to the center of the h + distribution, this indicates a relative period of quiescence for the early time GW emission.As the supernova develops, GW amplitudes become more pronounced, shown through the brighter points (later times) near the outer parts of the distribution.We do not observe any particular structure or morphology from this distribution for the s40o0 model for either the plus or cross polarization.
The slowly rotating model, s40o0.5, is represented in Figure 8b.The slow rotating case exhibits noticeable h + amplitudes from the early evolution.These points originate from the bounce signal, generating a dark purple distribution along the equator.During the ringdown, one can imagine a strain surface, similar to Figure 3a rotating and tracing a path along the equator that modulates in distance from the axis of rotation.In the medium brightness, blue points display columnated behavior along the axis of rotation, which we attribute to the coherent matter motion circulating around the PNS due to rotation.The final few hundred ms of evolution displays a relatively stochastic preferred GW direction.In total, the entire distribution of preferred GW directions is slightly more prolate than the nonrotating case.
The plus polarization for the fast rotating model, s40o1, is represented in Figure 8c.Similar to the slow rotating case, the darkest points display preferred directions along the supernova equator, consistent with the dominant bounce signal.As time evolves, however, the dominant direction for h + becomes noticeably columnated along the axis of rotation.This conclusion is important because it breaks the narrative set by many axisymmetric CCSN GW works that investigate the detectability of GWs from CCSNe, assuming the 'optimal' source orientation for CCSNe lies along the equator.This is indeed true for detecting the bounce signal.However, as illustrated in Figure 8c, the dominant viewing angle changes ∼ 90 • to the axis of rotation.Furthermore, GW amplitudes along this direction are larger by the nearly a factor of two.These results draw two main conclusions: (1) GW amplitudes along the axis of rotation can be comparable, if not greater, than GW measurements viewed along the equator for rotating CCSNe.(2) There is no single optimal orientation to observe GWs over the entire duration of a CCSN.Point (1) is in accordance with Powell & Müller (2020) and Shibagaki et al. (2021), which is present in model s40o1.Similar to the slow rotating case, the dominant viewing angle also traces a coherent path around the axis of rotation.We notice this behavior for both h + and h × , Figure 8d.
As noted in Section 2, we also investigate the dominant viewing angle evolution for the 20 M ⊙ models from (O'Connor & Couch 2018).Similar to the nonrotating case   , 11, 12, 12.25, 23].(Top) Compared to the mesa models, this suite has similar GW amplitudes during the accretion phase, t pb ≲ 300 ms.After, the extreme matter asymmetries dominate the GW signal, by contributing to direct offset 'GW memory' effects.(Bottom) During the accretion phase, the preferred direction oscillates about ⟨cos 2 θ⟩ 10ms ∼ 0.5.As the matter asymmetries become more drastic, different morphologies and orientations affect directionality to settle at different values of ⟨cos 2 θ⟩ 10ms .We note the only nonexploding model of this set, v12.25, retains its stochastic, nonsecular behavior.
for the 40 M ⊙ model, they display an isotropic distribution of dominant viewing angles in h + as well.Models v[9BW,11,12,23] display linear structures protruding from isotropic distributions, which are caused by large explosion asymmetries.

QUANTIFYING THE DIRECTIONALITY OF MAX GW STRAIN
We quantify the directionality of the maximum GW amplitude by decomposing the S 3 plots into useful metrics, shown in Figure 9, Figure 10, and Figure 11.We begin by selecting a subset of the points in a 10 ms window.This window was chosen, the light travel time between both LIGO detectors, to provide an appropriate balance between capturing the stochasticity of CCSN GWs, while being able to follow the directional evolution of the GW signal.This window is then progressed forward in time, representing a moving average of the data from the S 3 plots.Figure 9 is for the s20 and s40 models.Figure 10 is for the nonrotating mesa models, and Figure 11 is for the v[ 9BW,11,12,12.25,23]models.

Quantifying GW Maxima
The top panel of all three figures refers to the mean value of the maximum GW amplitudes during this window, de-noted ⟨|h max + |⟩ 10 ms .Hence, it is a proxy for GW activity.Beginning with the simple case, we examine the s20 and s40 models.The top panel notes increasing average GW maxima for increased Ω 0 values assigned at collapse.While rotation can provide a stabilizing effect to convection and resulting GW signal in 2D (Endal & Sofia 1978;Pajkos et al. 2019), nonaxisymmetries may emerge in 3D, and can produce stronger GW amplitudes.As shown, this coherent matter motion circumscribing the PNS surface bolsters GW emission by roughly an order of magnitude between nonrotating and rapidly rotating cases.
In the top panel of Figure 10, the mesa models do not exhibit much GW activity, noted by the y axis scale smaller by nearly two orders of magnitude compared to Figure 9.Of the five models, m20p and m20pLR show early bursts of GW emission ∼ 210 ms post-bounce due to the accretion of velocity perturbations placed in their overlying shells.As noted in O' Connor & Couch (2018) the larger values of lateral (i.e., nonradial) components of kinetic energy excite stronger PNS oscillations.
The data set from Vartanyan et al. (2023) begins with relalow GW amplitudes.These models are not initialized with any artificial perturbations, relying on turbulent instabilities to seed PNS oscillations later in the CCSN evolution.After an initial quiescent phase, typical GW amplitudes of order 10 −23 are driven during the supernova accretion phase.For models v[9BW,11,12], beyond ∼ 350 ms pb, the GW amplitude becomes dominated by the 'memory effect', driven by large scale matter asymmetries in the explosion ejecta.This onset occurs around 1.75 sec pb for v23.Once again, we reiterate, the GW amplitudes calculated for the Vartanyan et al. (2023) models are only from matter contributions, not neutrinos; though, neutrino asymmetries can produce similar linear evolution in the GW amplitudes.

Quantifying GW Maxima Along the Poles
The bottom panel of Figure 9, Figure 10, and Figure 11 selects GW maxima points from S 3 plots and weights them by cos 2 θ; lastly, it is normalized by ⟨|h max + |⟩ 10 ms .This quantity tracks the tendency of the h + to cluster near the poles of the CCSN, or the axis of rotation.In the limit of infinite points isotropically distributed on the unit sphere, ⟨cos 2 θ⟩ = 1/2, marked by the black dashed line.For points tightly clustered around the poles, ⟨cos 2 θ⟩ ∼ 1.For GW amplitudes emitted near the plane of the supernova equator, ⟨cos 2 θ⟩ ∼ 0.
In Figure 9, model s40o0 displays small oscillations about the isotropic distribution, displaying no major directionality preference.For s40o0.5, ⟨cos 2 θ⟩ 10 ms begins near a value of zero, consistent with the bounce and ringdown producing maximum strains preferentially toward the supernova equator.Over the timescale of 100 ms, this GW emission transitions to preferential strains along the pole.From ∼ 100 ms pb to ∼ 400 ms pb, the preferred direction remains aligned with the axis of rotation.As noted in Pan et al. (2021), model s40o0.5 exhibits strong mass accretion ( Ṁ ) up until ∼ 400 ms pb.Paralleled with the S 3 plot in Figure 8b, this ∼ 300 ms interim corresponds to the dark blue, columnated collection of points.This is due to dynamics of PNS rotation, paired with oscillations from the high Ṁ , allowing PNS nonaxisymmetries to mimic rotating ellipsoid behavior-see Section 4.4.After this point, the mass accretion rate plateaus to a value of Ṁ ∼ 500 M ⊙ /sec.As mass accretion stagnates, the rate at which the PNS receives angular momentum plateaus as well.This reduction in the rate of angular momentum transfer allows for GW emission to deviate from preferential emission along the axis of rotation, allowing stochastic convective processes and downflows to excite GW emission isotropically, rather than from rotationally dominated dynamics of downflow.
The most rapidly rotating models: s40o1, s20o2, and s20o3, and they behave similarly.They begin with ⟨cos 2 θ⟩ 10 ms ∼ 0, due to the bounce and ringdown dynamics.Over the timescale of milliseconds, they transition directly to GW emission predominantly along the pole; the rapid amount of angular momentum accretion is responsible.As the PNS continuously spins up and receives perturbations from turbulent mass accretion, ⟨cos 2 θ⟩ 10 ms rapidly approaches a value of 1 for the duration of the simulation.This behavior can be attributed to continued matter motion around the PNS, induced by rotation.In particular, for model s40o1 the onset of the low T/|W| instability-noted in Pan et al. (2018)-provides this condition.While not formally forming a bar-mode, these perturbations, paired with rotation, provide the necessary conditions from Section 4.4, mimicking approximate rotating ellipsoid evolution.Figure 10 does not show noteworthy deviation from an isotropic distribution-that is, the evolution of the preferred GW viewing angle remains nearly isotropic for the nonrotating mesa models.
The CCSN set from Vartanyan et al. (2023) begins with relatively stochastic values, with two exceptions.During the prompt convective phase ∼ 50 ∼ 100 ms, models v11 and v23 show preferred GW emission along the equator.Continuing into the accretion phase, for all five models, PNS oscillations generate GW signals with ⟨cos 2 θ⟩ 10 ms oscillating about 0.5.Beyond ∼ 350 ms, the GW memory dominates the preferred direction of GW emission for models v [9BW,11,12].This direct offset signal settles around 1.75 sec pb for v23.As different models settle to different ⟨cos 2 θ⟩ 10 ms , the geometry of the matter motion dictates the final preferred viewing angle.The only exception is model v12.25, which retains stochastic, nonsecular ⟨cos 2 θ⟩ 10 ms behavior, as it continues accreting and does not successfully explode.As expected, highly asymmetric and energetic ejecta correspond to more extreme asymmetries in the GW emission.The variety in ⟨cos 2 θ⟩ 10 ms values in the bottom panel of Figure 11 shows yet another lens describing the highly nonspherical and almost chaotic differences that can manifest in CCSNe, due to differences in initial conditions.We note the secular, low frequency drifts seen in the ⟨cos 2 θ⟩ 10 ms values provide directionality considerations for future space-based GW observatories, sensitive to frequencies of order ≲ 10 Hz.When filtering out frequencies below 10 Hz, we see similar stochastic evolution in ⟨cos 2 θ⟩ 10 ms to the mesa suite of models in Figure 10.This behavior is expected, since the models from Vartanyan et al. (2023) do not exhibit rotationally dominated dynamics, similar to the mesa suite in Figure 10.

QUANTIFYING THE EVOLUTION OF STRAIN SURFACES
We now turn to the influence of polarization.Thus far, h + strain surface plots have allowed for an intuitive description of the physical mechanisms driving the GW generation and the evolution of the preferred direction.However, as predicted by general relativity, GW detectors will pick up a combination of h + and h × .Thus, we take inspiration from the numerical relativity community by decomposing the GW emission into specific modes; for example, see Boyle et al.
(2014).While the S 3 plots were useful for drawing connections between ellipsoid toy models, decomposing strain surface plots into spin weighted spherical harmonics grants an additional dimension of understanding.This is because the S 3 plots were constructed with GW maxima, neglecting the remainder of the surface, that indeed have nonzero GW strain.By observing s a lm values, we gain a better understanding by quantifying the GW strain of both polarizations over time.
We take the surfaces for both modes by constructing h = h + − ih × , where i is the imaginary number.With a real and imaginary component to each GW strain surface, we now convolve these complex numbers with a complex conjugate of the spin weighted spherical harmonics, as introduced in Section 2.3.A list of the spin -2 weighted spherical harmonics is found in the central column of Table 2.The right column of Table 2 lists the closed form expressions of −2 a 2m in terms of the reduced mass quadrupole moment time derivatives, Qij .The −2 a 20 mode is weighted by a factor of sin 2 θ, indicating preferential GW amplitudes along the equator, with no azimuthal dependence.For the m = ±1 modes, maxima occur for θ at π/3, 2π/3 and for ϕ at 0, π/2.Physically, these correspond to maxima occurring between the equator and axis of rotation.For −2 a 2±2 , these distri-   2 are the expressions for spin weighted spherical harmonics, −2 Y 2m .Physically, the more points there are closer to a value of 1, for a given m, the more similar the GW distribution resembles −2 Y 2m .
butions are weighted by (1 ± cos θ) 2 e ±2iϕ , which achieve maxima at θ set to integer multiples of π and ϕ set to integer multiples of π/2.Physically, this behavior will track strong GW amplitudes along the z axis (or axis of rotation for rotating models).
Figure 12 shows the spherical harmonic coefficients for model s40o1.Near bounce, the −2 a 20 mode dominates the GW strain.This is indicative of the strongest GW emission along the xy plane, when considering both polarizations, consistent with the visual provided in Figure 2a.After a quiescent period of ∼ 120 ms, moderate equatorial GW amplitudes arise.By contrast, −2 a 2±2 dominate the GW signal beyond ∼ 150 ms pb.This is due to rapid mass and angular momentum accretion, along with the presence of the low T/|W| instability, and is consistent with the columnated behavior in Figure 8c and Figure 8d.The contributions between the pole and equator, from the −2 a 2±1 modes, remain subdominant.In summary, the time scale that sets this transition is set by the transition of PNS dynamics from bounce and ringdown to excitations during the accretion phase: ∼ 150 ms.While comparing −2 a 2m s is illustrative for a single model, as different models emit differing amounts of GWs, a cross model comparison becomes difficult.As a remedy, we appeal to the power emitted in each mode.
We define the power in each mode m as the square of the real component of the spin weighted spherical harmonic coefficient, −2 P l|m| = ( −2 a l±m ) 2 .Note, in this context, we use the term 'power' not to represent the time rate of change of GW energy, but the relative contributions of each −2 a lm to the strain surface.In order to compare relative amplitudes between modes, we use a normalization factor P norm = i ( −2 a lmi ) 2 , which represents the power in all modes.Normalizing each power yields −2 Plm = ( −2 a lm ) 2 /P norm .As we are concerned with directionality, we group the power by the magnitude of m values |m|; explicitly,  The purpose of these figures is to gain a quick intuition to identify in which modes the power resides.For s40o0, Figure 13a shows relatively equal contributions for m = 0, |m| = 1, and |m| = 2 modes.By eye, there is a similar density of points near a value of 1 for all three rows.By contrast, for model s40o1 in Figure 13b, more interesting features emerge.In the bottom row, PEq shows a cluster of points near 1 for the first 20 ms.This is indicative of the bounce signal and ringdown.Moving attention to the middle row, PMid shows a distinct lack of power from ∼ 150 ms to ∼ 350 ms.Simultaneously in the top panel, there is a surplus of power in the PPo modes-rotation has facilitated a transition of power from the |m| = 1 modes to the |m| = 2 modes.
Figure 14a shows the same quantities, for v23.For the first ∼ 1.5 seconds, similar to the nonrotating s40o0, the accretion phase exhibits relatively equal emission in the m = 0, m = ±1, and m = ±2 modes.However, as the main source of GW emission transitions from PNS oscillations to explosion ejecta, the modes settle into a secular drift because the matter source of the GW amplitude is now dominated by expanding matter, rather than nearly chaotic PNS oscillations.
Noting the frequency sensitivity of current-age GW detectors to frequencies ≳ 10 Hz, in Figure 14b, we generate the same plot, after passing the Qij time series data through a 10 Hz high pass filter (specifically, a tenth order Butter filter using SciPy.signal(Jones et al. 2001-).We note a similar distribution of points to the nonrotating model s40o0 seen in Figure 13a.This behavior is expected as there is no coherent matter motion dominating the dynamics of the PNS, leaving a relatively stochastic signal.The difference in the directionality of the GW signal between the filtered and unfiltered v23 signal emphasizes that future space-based GW observatories, sensitive to signals of order ≲ 10 Hz, will be subject to explosion geometry-dependent directional factors.
Lastly, we address how these transitions compound over the entire supernova evolution.After calculating −2 P2|m| for all models over their entire duration, we integrate these quantities in time: dt −2 P2|m| .To obtain a normalization factor, we add all integrated powers together dt −2 Ptotal = i dt −2 P2|mi| .The quantity, dt −2 P2|m| / dt −2 Ptotal shows the relative contributions of each |m| mode over the entire CCSN evolution.For all models, this quantity is plotted in Figure 15.The datasets are grouped into thirds, separated by two vertical dashed lines.The leftmost section indicates the mesa models.These models act as a control case to establish behavior for nonrotating, failed supernovae.13 and Figure 14a for all models in this study, normalized by the total integrated power over all m values.Physically, this value corresponds to what fraction of time the GW strain surface resembles a given spin weighted spherical harmonic.Our datasets are split into thirds by the two vertical dashed lines.The left section indicates results from the mesa suite of models.These nonrotating models act as control cases, displaying equal contributions among m = 0, m = ±1, and m = ±2 modes.Between both dashed lines are the rotating models, listed in order of initial rotation rate.As rotation increases, there is a transition of GW amplitudes away from |m| = 1 modes (lower values of squares), towards |m| = 2 modes (increasing values of triangles).The right section contains the models that show varying degrees of asymmetry of the supernova ejecta.In particular, v9BW displays GW amplitudes resembling −2 Y 2±1 , whereas v11 GW distributions more strongly resemble −2 Y 2±2 , due to unique explosion geometries.While the remaining models show varying degrees of asymmetry, they show different contributions due to not only different morphologies, but also different relative orientations.In summary, two main factors that influence the directionality of GW emission from CCSNe are the degree of rotation, along with explosion morphology and explosion orientation.Among all models, the m = 0, |m| = 1, and |m| = 2 modes contain equal contributions ∼ 1/3.
Between both sets of dashed lines are models s40o[0,0.5,1]and s20o [2,3].The left-most model is nonrotating, and the rightmost model has the most rapid rotation.The nonrotating s40o0 behaves similar to the mesa models, as expected.With moderate rotation, s40o0.5 shows a slight decrease in the |m| = 1 modes, and slight increase in the |m| = 2 modes.Continuing rightward, the models with increasingly rapid initial rotation exhibit more power in the |m| = 2 mode and less in the |m| = 1 modes.We note relative strength of the m = 0 remains around ∼ 1/3, refining our physical picture.While rotation helps create stronger amplitudes along the axis of rotation, the contributions of GWs along the equator still remain important.
To the right of both dashed lines contains models v[ 9BW,11,12,12.25,23].Models v12, v12.25 show similar behavior to the mesa suite as well because they do not display drastic memory effects in the GW signal (see lower panel of Figure ( 11)).However, models v9BW, v11, v23 show a different picture.Model v9BW exhibits most fractional power in the |m| = 1 modes and roughly half as much in the m = 0 and |m| = 2 modes.This inversion is due to the effect of the direct offset GW signal.In particular, it exhibits higher values of PMid for the evolution, particularly after ∼ 400 ms pb.The physical cause of the imprint is the explosion ejecta geometry favoring GW emission in a direction between the equator and pole.By contrast, model v11 shows the large majority of its power in the |m| = 2 modes, with slightly more relative power in the |m| = 0 mode, compared to |m| = 1.Model v23 shows yet another unique ordering, with most power being deposited in the |m| = 1 modes, then m = 0, and least in |m| = 2. Due to the spread in dt −2 P2|m| / dt −2 Ptotal , this implies, in general, for each exploding supernova in Nature, these relative contributions will change, depending on the explosion outcome and geometry.
For this section, we summarize: two important factors when considering the directionality of GWs-for both polarizations-are the degree of rotation and the explosion geometry.Rotation has a tendency to transition stronger GW amplitudes along the pole.The directionality of GW amplitudes for matter contributions to direct offset signals is dependent on the explosion morphology and orientation.
As a final note, we emphasize the potential for spin weighted spherical harmonics to analyze GW directionality data for future studies.In the binary merger community, instead of multiple TDWFs at different viewing angles, −2 a lm coefficients are stored (Boyle et al. 2019).Then, when a waveform at a specific viewing angle is needed, the coefficients can be used as weights for the −2 Y lm expressions, for example given in Table 2; they are summed together, and a TDWF is generated.In practice, for CCSN models that use the quadrupole formula for generating GW data, sharing the time derivative(s) of the mass quadrupole moments would provide the most precise angular GW information because Equation (1) and Equation ( 2) can be applied without information loss during the integral over angles for the −2 a lm calculation.However, the value of tracking −2 a lm s-or calculating them from quadrupole data-is that they provide a metric to quantify the directionality of the GW distribution, something the quadrupole data cannot concisely provide.For those interested in performing their own decomposition of GW surfaces, we point them to the right column of Table 2, which contains the complex valued, closed form −2 a 2m expressions (Ajith et al. 2007).For this work, we are concerned with general angular distributions of GWs, that is, GW amplitudes along the equator, midlatitudes, or pole, for which l = 2 modes are sufficient.While an interesting question, investigating the contributions from l > 2 modes to the CCSN GW signal (in particular, for codes that directly evolve spacetime) is beyond the scope of this work and will be saved for future studies.
8. DISCUSSION AND SUMMARY

Relating GW Directionality to Specific Phases of Fluid Evolution
There are a variety of physical mechanisms responsible for driving the PNS to produce GWs.Vartanyan et al. (2023) mention the presence of a prompt convective GW signal present in ∼ 50 ms following bounce.The authors establish the strain during this phase is only present in h + .After examining v11, we note the distribution of h + is concentrated in the equator, or strongly favors the m = 0 spin weighted spherical harmonic.This behavior is indicative of nearly axisymmetric dynamics.While the negative entropy gradient near the stalled shock, just after bounce, is responsible for generating this prompt convection, the GW distribution appears similar to the strain surface plot at bounce in Figure 2a.In the bottom panel of Figure 11, between 0.05 and 0.1 sec, we note values for ⟨cos 2 θ⟩ 10 ms < 0.5 for most of the interim, corroborating these observations.Pan et al. (2021) also note several potential factors such as the SASI and the low T/|W| instability.In particular, it has been shown that the spiral modes of the SASI can redistribute angular momentum, causing the PNS to spin up in an initially nonrotating progenitor (Guilet & Fernández 2014;Pan et al. 2021).As the GW directionality for the mesa models with and without SASI were explored, with no clear preferred direction, we do not observe dominant GW viewing angle evolution associated with this hydrodynamic instability.This can be expected as O 'Connor & Couch (2018) note models exhibiting SASI,m[20,20LR,20p,20pLR], only yield slow rotating, ∼ 1 sec period PNSs from SASI spin-up.
For the initially nonrotating model s40o0, Pan et al. (2021) note a PNS forming with spin period ∼ 0.1 sec.We do not see preferred GW directionality with this model either (see bottom panel of Figure 9 and Figure 15).We conclude it is not the mere presence of a rotating PNS, but the accretion of significant amounts of angular momentum, which assist in directing stronger GW amplitudes along the axis of rotation.
In the fast rotating model, Pan et al. (2021) note the possible emergence of the low T/|W| instability, based on low frequency GW emission beginning ∼ 150 ms pb.As the fast rotating model exhibits a preferred direction of GW emission along the axis of rotation and coherent paths in the evolution of the preferred direction, it is possible this instability contributes.However, the presence of the low T/|W| instability is not requisite for these two behaviors as the slow rotating model-which shows no sign of low T/|W|-also remains slightly columnated and exhibits clear paths in Figure 8b.
In relation to the mesa models, we do not notice any coherent paths in Figure 10 for preferred GW directions during instabilities like the SASI or Lepton number Self-sustained Asymmetry (LESA).

Implications for Observability
These directional dependencies help prepare for the first direct detection of GWs from CCSNe.As detector sensitivities improve and signal search algorithms become more advanced, identifying which components of the signal can be reconstructed is vital.Importantly, the observability of a GW depends on not only the physical mechanism generating it but also the angle at which it is observed.For a rotating supernova, the bounce signal is most visible when viewed along the supernova equator.However, potentially stronger GW amplitudes along the axis of rotation may be missed during the accretion phase.By contrast, if viewed along the axis of rotation, the broadband bounce signal will not be detected, but high frequency oscillations from the PNS may have stronger amplitudes.For a randomly oriented su-0.0040.002 0.000 0.002 0.004 The light green, solid, curved line is the strain surface at t pb = −0.623ms for the 10 kpc source distance.The dark blue line is for the 80 kpc source distance.Shaded areas show the viewing angles where |h| > 10 −22 for the two source distances.As the CCSN source distance increases, the opening angle for detectable GW strains decreases; the darker blue 'detectable' region traces a smaller angle (∼ 53 • ) than the lighter green 'detectable' region (∼ 143 • ).pernova, the viewing angle will likely be off-axis, so GW detection algorithms should search for both components of the signal.
During the accretion phase, the GW emission remains roughly isotropic for nonrotating or slowly rotating casesi.e., the majority of expected supernovae.For rarer, rapid rotators, GW amplitudes along the pole can be stronger, by nearly a factor of two, compared to the equator.When considering rapid rotators, a balance must be considered.For the bounce phase, a well templated GW strain can be detected, with relatively weaker GW amplitudes during the accretion phase when viewed along the equator.Conversely, no well templated bounce signal would be detected when viewed along the pole but stronger accretion phase GW amplitudes would arise.This sensitive interplay as to which direction gives the best chance for detection depends on the orientation, source distance, supernova dynamics, and detection algorithm (Szczepańczyk & Zanolin 2022).
Lastly, for successful explosions, the explosion geometry will ultimately dictate along which direction a maximum GW amplitude is emitted.Likewise, it will produce certain direc-tions along which GW minima drop below detector sensitivity levels.Since the explosion ejecta morphology relies upon the interaction between initial progenitor convective perturbations, magnetized fluid evolution, and neutrino transport, we expect this effect to vary as widely as the progenitors themselves.While our work examined the influence of matter contributions to the GW signal, asymmetric neutrino production would also provide measurable direct offset signals, which also vary by viewing angle (Vartanyan et al. 2023).In this context, an exciting possibility exists for multiwavelength astronomy with the future space-based GW observatories, such as the Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al. 2017), to use high frequency GW observations with the LIGO-Virgo-KAGRA network (Abbott et al. 2018)-either during the bounce or accretion phaseas a warning sign for LISA to expect a low frequency counterpart in the following ∼second.
Contributions to the GW background are also affected by consequences of directionality.For detectors attempting to quantify the spectrum of the GW background (Christensen 2019), this work offers motivation to consider CCSN popu-lation studies.While supernova GW signals produce weaker amplitude GWs compared to binary mergers (Finkel et al. 2022), theoretical works could construct populations of randomly oriented and distributed CCSNe, with varying rotation rates and explosion fates, to better quantify which components of the signal may be detectable for GWs from a superposition of sources.

Applying Directionality Evolution to the Bounce Phase
To provide a concrete example considering the observability of GW amplitudes, we consider the bounce phase of model s40o1.The left panel of Figure 16 shows the GW strain |h| = h 2 + + h 2 × during the bounce.The light green line is for a CCSN 10 kpc from Earth, and the dark blue line is for a CCSN distance of 80 kpc.To produce a strain surface plot, we select a time when the light green line produces |h| ∼ 10 −21 , a GW amplitude above which most of the bounce signal is present.
In the right panel of Figure 16, we show a 2D strain surface plot at t pb ∼ −0.6 ms.The light green line is the strain surface for a source distance of 10 kpc, and the blue line is for 80 kpc.Suppose a GW detector needed a strain of |h| ∼ 10 −22 to robustly reconstruct a bounce signal.Then, the colored lines greater than 10 −22 away from the origin (outside the grey dashed circle) would represent viewing angles along which the GW amplitude is 'detectable'.For ease, we have marked the shaded regions as the viewing angles along which the signal would be detectable.The opening angle of detectability is ∆θ 10kpc detect ∼ 142 • for the 10 kpc (light green) case and ∆θ 80kpc detect ∼ 53 • for the 80 kpc (dark blue) case.An opening angle of ∆θ 10kpc detect ∼ 142 • corresponds to a solid angle of ∆Ω 10kpc detect = 11.9 steradians.An opening angle of ∆θ 80kpc detect ∼ 53 • corresponds to a solid angle of ∆Ω 80kpc detect ∼ 5.6 steradians.Armed with these numbers, consider a CCSN event at 10 kpc and 80 kpc.Assuming a random orientation for each, the likelihood the Earth falls within the GW detectability opening angle is 94% (∆Ω 10kpc detect /4π) for the 10 kpc case.When increasing the source distance out to 80 kpc, this chance drops to 42% (∆Ω 80kpc detect /4π).We emphasize these are not actual detection probabilities that consider signal reconstructions and frequency sensitivities.This is a simplified, worked example to show how the directionality work here can be used to better understand how viewing angle affects the likelihood of detecting CCSN GWs.We leave more robust population estimates that use multiple sources, consider different GW distributions beyond the bounce phase, and consider signal reconstruction algorithms for future work.

A Note on Signal Injection
Because of the aforementioned directional considerations, we encourage the observational community to continue us-ing gravitational waveforms from 3D studies injected at a variety of random angles (e.g., Szczepańczyk et al. 2023), rather than injecting an angle averaged amplitude.As we have shown, depending on the source characteristics and supernova phase, the expected GW amplitudes may not be isotropic over a given phase.An angle averaged GW amplitude may be appropriate for slow/non rotating CCSNe, during the accretion phase, as the emission is roughly isotropic over timescales of hundreds of milliseconds.The bounce amplitude, by contrast, varies with sin 2 θ pole , where θ pole is the angle between the pole, or axis of rotation, and the viewing angle.Consider the 80 kpc example from the previous Subsection 8.3.Sufficiently strong GW amplitudes ≳ 10 −22 can be detected for ∼ 42% of the viewing angles.However, if an angle averaged GW amplitude, h avg , were used, h avg ∼ π 0 1.25 × 10 −22 sin 2 θdθ/ π 0 dθ = 6.25 × 10 −23 , this source would fall below the 10 −22 threshold for all viewing angles, yielding a 0% probability of detection.Similar examples with closed form expressions for GW amplitudes, like rapidly rotating models during the accretion phase-Equation ( 10) and Equation ( 11)-or late time exploding models (for next generation, low frequency detectors)-see Richardson et al. (2022)-emphasize the importance of using a distribution of GW signal injections, rather than injecting an angle averaged quantity.In particular, they give more proper estimates of GW detectability for asymmetric explosions, like those seen in Vartanyan et al. (2023) and rapidly rotating models, like those seen in Powell & Müller (2020) and Powell et al. (2023).

Summary
This work has considered a new factor regarding the observability of GWs from CCSNe, particularly the impact of viewing angle.In particular, this work serves as a bridge between the theory and observation communities: explaining how the source physics influences viewing angle effects on CCSN GWs and quantifying the time evolution of preferred directions of GW emission.Here, we review the findings from our work: • strain surface plots can characterize the nature of the mechanism generating GWs at a given point in time, Figure 6, • viewing the collective distribution of preferred directions for GW emission shows CCSNe do not have a single 'optimal' viewing angle, Figure 8, • in particular, when considering all viewing angles, we quantitatively show nonrotating or slowly rotating cases exhibit roughly isotropic GW amplitudes throughout their accretion phase, whose prominent direction of emission transitions on the dynamical timescale of the PNS, Figure 10, • depending on the degree of rotation, when considering all viewing angles, we quantitatively show rotating CCSNe can exhibit strong GW amplitudes along the equator, then along the axis of rotation, with a transition timescale set by the time between ringdown and the accretion phase (order tens of ms), Figure 9, Figure 12, and Figure 13, • the dominant viewing angle of GW emission in rotating CCSNe may follow coherent paths that precess around the axis of rotation, Figure 8c, • while PNSs are not solid-body ellipsoids, simplified analytic models, motivated by PNS dynamics, can describe each phase of GW generation (see Richardson et al. (2022) for memory effects); if the geometry of a simple PNS deformation-or other ellipsoid-like astrophysical scenarios-can be predicted, a similar toy model can provide an estimate of the resulting GW directionality, Appendix A, • when considering spin weighted spherical harmonics (Table 2), rotation facilitates a transition of GW strain surfaces (e.g., Figure 3a) from m = ±1 to m = ±2 modes, or focusing stronger GW amplitudes closer to the axis of rotation, Figure 13, • two major factors influencing the directionality of h + and h × from matter sources are the degree of rotation and geometry of the explosion ejecta, Figure 15, • these two factors can inform future population studies that include CCSN contributions to the GW background, • we recommend future CCSN GW theory investigations to report both h + and h × , as the distribution h 2 + + h 2 × will remain invariant for different detectors, whereas the individual distributions of h + and h × will change for different sites.For example, a GW bounce signal with mostly h + at one detector may have significant h × components at a detector oriented differently, • we note that the use of angle averaged GW strains can affect detectability estimates for rapidly rotating and asymmetric exploding models; injecting TDWFs at a variety of source angles is recommended.
By using high fidelity 3D simulations of CCSNe from multiple code bases, we have reviewed a new method towards visualizing, and consequently analyzing, GW emission from CCSNe.While this study uses CCSNe, in general, any GW source can make use of strain surface plots.For those who have Qij from their own simulations, and wish to decompose the GW strain surfaces into −2 a 2m (to construct similar figures above), we recommend using the closed form expressions for spin weighted spherical harmonics specified in the right column of Table 2.
The strength of this work lies in the diversity of computational grids, initial progenitors, EOSs, rotation rates, neutrino treatments, initial perturbations, and unique explosion morphologies.Nevertheless, we acknowledge there is room for improvement.One major component is the influence of magnetic fields.These fields facilitate different fluid dynamics and provide another avenue for angular momentum transport away from the PNS (e.g., Powell et al. (2023)).Furthermore, including fully dynamical spacetimes would allow us to follow the CCSN evolution closer to, and beyond, black hole formation.While many of these models use sophisticated neutrino physics-M1-accounting for effects such as transitions between neutrino flavors-and associated flavor instabilities-remains an important factor not yet incorporated into 3D CCSN simulations.
Acknowledging these caveats, this study offers another observable factor to consider when generating gravitational waveforms from CCSNe.By considering both theoretical factors like CCSN source physics, alongside observational factors like the impact of the viewing angle and different detector sites, the GW community will be better prepared for the first direct detection of GWs from stellar explosions.
imperative.The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago.
Software: FLASH (see footnote 7) (Fryxell et al. 2000(Fryxell et al. , 2010)) In partial generality, we describe a toy system in which a tilted ellipsoid changes size and 'wobbles' around the axis of rotation-in our case, we select the z axis.That is, in general, the axis of rotation does not align with any of the three axes of symmetry, and the ellipsoid radii are not necessarily constant in time.We also neglect mass changes, ignoring the M and Ṁ terms from above.To account for a tilted ellipsoid, similar to Creighton & Anderson (2011), we apply a rotation matrix R to the quadrupole moment matrix I about the x axis, by a tilt angle τ From Equations (A33-A34) the strain surface plots for the wobbling, oscillating ellipsoid are constructed.
B. ANGULAR VELOCITIES OF NUMERICAL SIMULATIONS In Section 4.3 and Section 4.4, we note a mismatch between the amplitudes of the strain surface plots for the simulations and toy model.We attribute this difference to the solid-body assumption of the toy model.In Figure 17 we provide angular velocities for model s40o1 during the ringdown phase (∼ 13 ms pb) on the left and the accretion phase (∼ 285 ms pb) on the right.We denote the 10 11 g cm −3 density contour in white.The small white patches within each slice plot are due to negative angular velocities (clockwise rotation) not rendering on the log scale.
To justify our reasoning regarding differences in the GW amplitudes between the strain surfaces from the simulations and toy models, we provide slices of the angular velocities.In the left panel, we note a rapidly rotating inner 10 km, with nonuniformities in the angular velocity in the outer 90 km.For the accretion phase case on the right, we notice differences that vary radially outward.Neither case exhibits a uniform, solid-body angular velocity profile.Thus, we justify our difference in amplitudes for the strain surface plots in Figure 3   AngularVelocity (1/s)

Figure 1 :
Figure 1: Time domain waveforms (plus polarization) for s40o0 (top), s40o0.5 (middle), and s40o1 (bottom) (Pan et al. 2021).The left (right) column corresponds to an observed h + when viewed along the equator (pole).While TDWFs are useful for analyzing GW emission throughout the supernova simulation, they are inherently limited to GW signals at a single, fixed viewing angle.This motivates the need for additional analysis methods that can identify a dominant viewing angle beyond the equator or pole and determine if this direction evolves in time.Note the different scales between each y axis.Figure from Pajkos (2022).
+ , t pb = 2.02 ms (a) Strain surface plot for h+ just after bounce.The star denotes the direction along which the largest h+ would be received, what we denote as the 'preferred' viewing angle.Strain surface plot for a toy model of an axisymmetric prolate ellipsoid of mass 1 M⊙ with polar radius c = 62 km and equatorial radii a = b = 60 km, with velocity at the pole ċ = 600 km s −1 , ȧ = ḃ = 0 km s −1 , replicating the strain surface distribution for bounce seen in Figure 2a.We note the fluid-like nature of the PNS (deviating from a solid-body ellipsoid) makes assigning a single polar velocity difficult-thus the discrepancies in the magnitude of the estimate.Nevertheless, our goal is to describe the geometry of the GW distribution, which agrees quite well.
Entropy slice through the equator from model s40o1 2 ms after bounce display the axisymmetry of the PNS before convective instabilities create asymmetries.

Figure 2 :
Figure2: Strain surface plots and entropy slices for model s40o1 during the bounce phase.The x and y axes form the equator of the supernova; the z axis indicates the axis of rotation.For the top two panels, the distance from the origin to a point on the green surface represents the detected h + along that direction.The bottom two panels are entropy slices, with overlaid velocity profiles (arrows).
h + , t pb = 13.36 ms (a) Strain surface plot for h+ during ringdown.This rotating structure, which will modulate in scale for tens of ms during ringdown, will give the largest GW amplitudes along the equator of the supernova.(b) Strain surface plot for a toy model of an 'elastic ellipsoid'.Numbers taken from model s40o1 inform toy ellipsoid parameters of mass 1.3 M⊙, with principal axes c = 77 km, a = 83 km, b = 81, rotating with Ω = 100 rad s −1 , ȧ = ḃ = 5000 km s −1 , ċ = 0 km s −1 , which replicates the strain surface for rotating ringdown seen in Figure3a.In particular, the rotating nonaxisymmetry (i.e., a ̸ = b) of the ellipsoid contributes to the 'four lobed' structure.The changing principal axes ( ȧ ̸ = 0, ḃ ̸ = 0) contribute to the protruding GW surface along the equator.Entropy slice through the pole from model s40o1 ∼ 13 ms after bounce shows a slightly oblate PNS surface.For scale, every cm on the page corresponds to 5.6 × 10 9 cm s −1 .Entropy slice through the equator from model s40o1 ∼ 13 ms after bounce display the PNS surface deviating from axisymmetry (a circle in this cross section) as convective instabilities within the PNS create asymmetries.

Figure 3 :
Figure 3: Strain surface plots and entropy slices for model s40o1 during the ringdown phase.The x and y axes form the equator of the supernova; the z axis indicates the axis of rotation.For the top two panels, the distance from the origin to a point on the green surface represents the detected h + along that direction.The bottom two panels are entropy slices, with overlaid velocity profiles (arrows).
h + , t pb = 286.64ms (a) Strain surface plot for h+ during the accretion phase (∼ 285 ms pb).This rotating, four-lobed structure persists for the majority of the accretion phase, producing the largest GW amplitudes along the axis of rotation.(b) Strain surface plot for a 'rigid, solid-body ellipsoid'.Data taken from simulation outputs motivate toy model parameters of mass 2.1 M⊙, with principal axes c = 40 km, a = 62 km, b = 55 km, rotating with Ω =1000 rad s −1 , ȧ = ḃ = ċ = 0 km s −1 , replicating the strain surface during the accretion phase, seen in Figure4a.Once again, our choice of the 'solid-body' approximation captures the morphology of the GW strain, while overestimating the magnitude.Entropy slice through the CCSN pole for model s40o1 ∼ 285 ms pb.Arrows show the magnitude of the velocity in the yz plane, and the solid black line is the density contour for the PNS surface of ρ = 10 11 g cm −3 .After the accretion of sufficient angular momentum from overlying material, the oblateness of the PNS increases.For scale, every cm on the page corresponds to 5.6 × 10 9 cm s −1 .Entropy slice through the equator from model s40o1 ∼ 285 ms after bounce.Accretion downflows striking the PNS surface cause slight deviations from axisymmetry (a circle in this slice).The oblateness of the PNS, deviation of axisymmetry in the azimuthal direction, and rotation motivate our choice for modeling the source of the GW emission as a rotating ellipsoid.

Figure 5 :
Figure 5: Magnitude of the spatial velocity-color values centered ∼ 3.9 × 10 9 cm s −1 -near the PNS for model s40o1, taken ∼ 286 ms pb.The panel spans 200 km.The axis of rotation is aligned vertically with the page.The 3D structure more clearly displays the coherent flows around the PNS, whose oblateness can be seen by the inner contour, complementing the entropy slices seen in Figure 4. model parameters M = 2.1M ⊙ , c = 40 km, a = 62 km, b = 55 km, ȧ = ḃ = ċ = 0, and Ω = 1000 rad s −1 .Observing the velocity field in Figure4cand Figure4dmotivates our approximation of ȧ = ḃ = ċ = 0.In reality, the PNS will be dynamic; however, this assumption implies that rotation is dominant, or that the rotation of the nonaxisymmetric features (rather than their radial motion) contribute more importantly to the 2M Ṙ2 term in Equation (7).Figure4bshows the GW distribution for the solid rotating ellipsoid, displaying a roughly four lobed 'clover-like' structure, with relative GW maxima along the axis of rotation.To explain the behavior in the strain surface plots, we appeal to Equation (10).Entering as 1 + cos 2 θ, for θ = 0, π this term exhibits a maximum, consistent with the axis of rotation.To explain the four lobed structure, observe the cos 2ϕ term.Similar to the explanation during the ringdown phase, as an observer circumscribes the CCSN in the azimuthal di- h + , t pb = 229.76ms (a) Strain surface plot of h+ for nonrotating model s40o0 during the accretion phase (∼ 230 ms pb).We attribute the skewed maximum GW amplitude (not aligned with the coordinate axes) with downflows perturbing the PNS at a misaligned angle.Strain surface plot for nonrotating, tilted ellipsoid.Data taken from simulation outputs motivate toy model parameters of mass 2.1 M⊙, with principal axes a = b = c = 45 km, nonrotating, ȧ = ḃ = 0 km s −1 , and ċ = 50 km s −1 , tilted by π/4 radians, replicating the strain surface during the accretion phase, seen in Figure 6a.Once again, our choice of the toy model approximation captures the morphology of the GW strain, while differing in the magnitude.Entropy slice through the CCSN pole for model s40o0 ∼ 230 ms pb.Arrows show the magnitude of the velocity in the yz plane, and the solid black line is the density contour for the PNS surface of ρ = 10 11 g cm −3 .The skewed angle of the post shock region is what motivates our choice for the tilted ellipsoid.For scale, every cm on the page corresponds to 5.6 × 10 9 cm s −1 .Entropy slice through the equator from model s40o0 ∼ 230 ms after bounce.Similar asymmetries to Figure 6c of the shock are visible when viewing the supernova equator.

Figure 7 :
Figure7: Strain surface plots of the tilted toy model ellipsoid from Section 4.5.h + , h × , and h 2 + + h 2 × are shown from left to right.While similar to the dynamics of the bounce phase, we note different contributions to each polarization.Whereas a coordinate system aligned with the principal ellipsoid axis along c will only produce h + emission (Figure2b), a differently Dominant viewing angles of h+ for s40o0.The distribution of dominant GW viewing angles remains nearly isotropic for the duration of the supernova evolution.Dominant viewing angles of h+ for s40o0.5.The dark (early time) points distributed along the equatorial plane are generated from the bounce signal and ringdown.The dominant signal direction then aligns with axis of rotation and finally relaxes to an isotropic distribution.Dominant viewing angles of h+ for the fast rotating model s40o1.The preferred GW direction is noticeably columnated and follows coherent paths that precess around the axis of rotation.Same fast rotating model as in Figure8b, however, for the cross polarization h×.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :
Figure8: These strain surface scatter plots, or S 3 plots, compile the direction of maximum GW emission (e.g., purple stars in Figure2a, Figure3a, Figure4a, and Figure6a) throughout the supernova evolution from simulation outputs.Brighter points correspond to later times.The xy plane forms the equator of the supernova.The z axis identifies the axis of rotation for rotating cases.

Figure 12 :
Figure 12: The real value of the decomposition of h + − ih × , scaled by the distance 10 kpc, into spin weighted spherical harmonics −2 a 2m .Different rows correspond to different m numbers.Listed in Table 2 are the expressions for spin weighted spherical harmonics.At early times, −2 a 20 shows dominant contributions to the GW surfaces.Because of the sin 2 θ behavior, as noted in Table 2, this behavior indicates dominant GW emission along the equator of the supernova-characteristic of the axisymmetric bounce signal.Beyond 150 ms pb, the m = ±2 mode shows dominant contributions along the poles.m By eye, model s40o0 shows relatively equal contributions (i.e., similar density of points) from PPo , PMid , and PEq; physically, contributions from all three modes correspond to more isotropic emission.Compared to the top panel of Figure 13a, PPo for model s40o1 shows a larger contribution to the overall GW signal emitted.These increased m = ±2 modes indicate increased GW amplitudes closer to the axis of rotation.

Figure 13 :
Figure 13: The real value of the decomposition of h + − ih × into spin weighted spherical harmonics −2 a 2m .These coefficients are then squared and normalized by their sum, or −2 P2|m| = ( −2 a 2±m ) 2 / i ( −2 a 2±mi ) 2 .For convenience, we introduce, PPo = −2 P2|2| , favoring emission along the poles, PMid = −2 P2|1| , favoring emission along the midlatitudes, and PEq = −2 P20 , favoring emission along the equator.Different rows correspond to different magnitudes of m.Listed in Table2are the expressions for spin weighted spherical harmonics, −2 Y 2m .Physically, the more points there are closer to a value of 1, for a given m, the more similar the GW distribution resembles −2 Y 2m .
−2 Pl|m| = −2 Plm + −2 Pl(−m) .For convenience, we introduce shorthand symbols, PPo = −2 P2|2| which describes what fraction of the strain surface resembles −2 Y 2±2 , favoring emission along the poles.PMid = −2 P2|1| describes what fraction of the strain surface resembles −2 Y 2±1 , emission between the pole While the PNS oscillations provide the stochasticity in the signal, after t pb ∼ 1.5 seconds the asymmetry of the explosion ejecta dominates the GW signal, this secular drift of material provides slow changes in PEq and PMid , compared to the models in Figure 13.Similar to the left panel, however, after passing the signal through a 10 Hz high-pass filter.As expected, the power is evenly distributed among PEq, PMid , and PPo modes, similar to model s40o0 in Figure 13a.

Figure 14 :
Figure 14: Same as Figure 13, however for model v23.and equator-the midlatitudes.PEq = P20 describes what fraction of the strain surface resembles −2 Y 20 , emission along the equator.Armed with these quantities, we turn to Figure13.Each row corresponds to −2 Pl|m| for a given m.The purpose of these figures is to gain a quick intuition to identify in which modes the power resides.For s40o0, Figure13ashows relatively equal contributions for m = 0, |m| = 1, and |m| = 2 modes.By eye, there is a similar density of points near a value of 1 for all three rows.By contrast, for model s40o1 in Figure13b, more interesting features emerge.In the bottom row, PEq shows a cluster of points near 1 for the first 20 ms.This is indicative of the bounce signal and ringdown.Moving attention to the middle row, PMid shows a distinct lack of power from ∼ 150 ms to ∼ 350 ms.Simultaneously in the top panel, there is a surplus of power in the PPo modes-rotation has facilitated a transition of power from the |m| = 1 modes to the |m| = 2 modes.Figure14ashows the same quantities, for v23.For the first ∼ 1.5 seconds, similar to the nonrotating s40o0, the accretion phase exhibits relatively equal emission in the m = 0, m = ±1, and m = ±2 modes.However, as the main source of GW emission transitions from PNS oscillations to explosion ejecta, the modes settle into a secular drift because the matter source of the GW amplitude is now dominated by expanding matter, rather than nearly chaotic PNS oscillations.

Figure 15 :
Figure15: Integrated values of −2 P2|m| -the time series data seen in Figure13and Figure14afor all models in this study, normalized by the total integrated power over all m values.Physically, this value corresponds to what fraction of time the GW strain surface resembles a given spin weighted spherical harmonic.Our datasets are split into thirds by the two vertical dashed lines.The left section indicates results from the mesa suite of models.These nonrotating models act as control cases, displaying equal contributions among m = 0, m = ±1, and m = ±2 modes.Between both dashed lines are the rotating models, listed in order of initial rotation rate.As rotation increases, there is a transition of GW amplitudes away from |m| = 1 modes (lower values of squares), towards |m| = 2 modes (increasing values of triangles).The right section contains the models that show varying degrees of asymmetry of the supernova ejecta.In particular, v9BW displays GW amplitudes resembling −2 Y 2±1 , whereas v11 GW distributions more strongly resemble −2 Y 2±2 , due to unique explosion geometries.While the remaining models show varying degrees of asymmetry, they show different contributions due to not only different morphologies, but also different relative orientations.In summary, two main factors that influence the directionality of GW emission from CCSNe are the degree of rotation, along with explosion morphology and explosion orientation.

Figure 16 :
Figure16: (Left) TDWF for model s40o1 when viewed along the equator of the CCSN.Light green color is for a source distance of 10 kpc and dark blue is for a source distance of 80 kpc.The vertical line marks t pb = −0.623ms, when the strain for the 10 kpc source reaches 10 −21 .(Right) Since the bounce creates an axisymmetric GW distribution, we show a slice through a strain surface plot-the x-z plane.The distance from the origin to any point in this plane is the GW strain along that viewing angle.The grey dashed circle shows a contour for |h| = h 2 + + h 2 × = 10 −22 .The light green, solid, curved line is the strain surface at t pb = −0.623ms for the 10 kpc source distance.The dark blue line is for the 80 kpc source distance.Shaded areas show the viewing angles where |h| > 10 −22 for the two source distances.As the CCSN source distance increases, the opening angle for detectable GW strains decreases; the darker blue 'detectable' region traces a smaller angle (∼ 53 • ) than the lighter green 'detectable' region (∼ 143 • ).

Figure 17 :
Figure 17: Nonuniformities in the angular velocity profiles for model s40o1 during the ringdown phase (left) and accretion phase (right).Slice plots show the equatorial plane.White lines indicate the rest mass density contour for 10 11 g cm −3 .

Table 1 :
Setup information for all 15 models in this study.Column labels represent the following: M-zero age main sequence mass, Ω 0 -central rotation rate at collapse, A differential rotation parameter, EOS, neutrino (ν) treatment, and t pb end -simulation end time (post-bounce), and corresponding reference.