TY - JOUR
AU1 - Torres-Forné,, Alejandro
AU2 - Cerdá-Durán,, Pablo
AU3 - Passamonti,, Andrea
AU4 - Obergaulinger,, Martin
AU5 - Font, José, A
AB - Abstract Improvements in ground-based advanced gravitational wave (GW) detectors may soon allow us to observe the GW signal of a nearby core-collapse supernova. For most progenitors, likely with slowly rotating cores, the dominant GW emission mechanisms are the post-bounce oscillations of the proto-neutron star (PNS) before the explosion. We present a new procedure to compute the eigenmodes of the system formed by the PNS and the stalled accretion shock in general relativity including space–time perturbations. We apply our analysis to two core-collapse simulations and show that our improved method is able to obtain eigenfrequencies that accurately match the features observed in the GW signal and to predict the qualitative behaviour of quasi-radial oscillations. Our analysis is possible thanks to a newly developed algorithm to classify the eigenmodes in different classes (f, p, and g modes), improving our previous results. We find that most of the GW energy is stored in the lowest-order eigenmodes, in particular in the 2g1 mode and in the 2f mode. Our results also suggest that a low-frequency component of the GW signal attributed in previous works to the characteristic frequency of the standing accretion shock instability should be identified as the fundamental quadrupolar f mode. We also develop a formalism to estimate the contribution of quasi-radial (l = 0) modes to the GW quadrupolar component in a deformed background, with application to rapidly rotating cores. This work provides further support for asteroseismology of core-collapse supernovae and the inference of PNS properties based on GW observations. asteroseismology, gravitational waves, methods: numerical, stars: neutron, stars: oscillations, supernovae: general 1 INTRODUCTION With the LIGO/Virgo discovery of gravitational waves (GWs) from a binary neutron star (BNS) merger and the subsequent follow-up observations across the electromagnetic spectrum by dozens of astronomical facilities and several neutrino telescopes, the field of multimessenger astronomy has started (Abbott et al. 2017b). Significant advances towards explaining a number of open issues in relativistic astrophysics have been made thanks to the single observation of GW170817, from the mechanism behind short gamma-ray bursts (Abbott et al. 2017c) to the r-process-mediated nucleosynthesis of heavy elements in kilonovae (Abbott et al. 2017d), along with independent measures of cosmological parameters (Abbott et al. 2017a). Moreover, GWs from BNS mergers offer the opportunity to probe the properties of high-density matter, yielding constraints on the neutron star radii and on the equation of state (EOS) (Abbott et al. 2018; Annala et al. 2018; De et al. 2018; Fattoyev, Piekarewicz & Horowitz 2018). A type of transient GW signal that remains to be detected is the one produced in core-collapse supernovae (CCSNe), associated with the catastrophic collapse of the unstable iron cores that massive stars (from ∼8 M⊙ to ∼100 M⊙) develop at the end of their cycles of thermonuclear burning reactions. Collapsing stars produce rich and complex gravitational waveforms. When detected, these may provide important information about the phenomenology of the scenario, specially when combined with observations of their electromagnetic emission and neutrino emission. Their typical range of frequencies falls within the most sensitive region of the advanced LIGO and Virgo detectors. CCSNe are rare events, happening at rates of about one per century (see e.g. Gossan et al. 2016) observable within the Milky Way. This may explain why GWs from CCSNe have not yet been detected. In terms of signal amplitudes, it would be possible to detect magneto-rotational explosions up to ∼5 Mpc (Gossan et al. 2016), but the event rate for such a class of explosions is even lower, |${\sim } 10^{-4}\, \rm {yr}^{-1}$| (less than |$1{{\ \rm per\ cent}}$| of all CCSNe; see discussion in Woosley & Bloom 2006). Despite this constraint, there have been attempts to detect the GW emission of supernova signals with the current network of detectors, either through techniques that search generic short-duration burst signals without targeting specific supernova events or with targeted searches that use the known sky location and time window of electromagnetic or neutrino events (Thrane & Coughlin 2015; Klimenko et al. 2016). Unlike the case of binary black holes (BBH) mergers, it is currently not possible to relate uniquely and unambiguously the properties of the progenitor stars (such as mass, rotation rate, metallicity, or magnetic fields) with the resulting waveforms. The reasons are diverse: (i) the complex non-linear dynamics associated with the evolution of a fluid interacting with neutrino radiation, (ii) the stochastic and chaotic behaviour of instabilities (both during and prior to the collapse of the star), (iii) the uncertainties in the stellar evolution of massive stars (specially regarding the treatment of convection, magnetic fields, and angular momentum transport), and (iv) the uncertainties in the nuclear and weak interactions necessary for the high-density EOS and neutrino radiation, respectively. In CCSNe, GWs are produced during the hydrodynamical bounce and the evolution of the instabilities that occur in the cavity formed by the newborn neutron star and the accretion shock. The violent dynamics excite different modes of oscillation of the proto-neutron star (PNS) and of its surroundings, including the shock wave. There is a wide literature on the topic of neutron star oscillations of different kinds (see e.g. Kokkotas & Schmidt 1999; Friedman & Stergioulas 2013, and references therein). In this work we focus on oscillation modes of hot and stratified PNSs. McDermott, van Horn & Scholl (1983) were the first to consider the oscillation modes of a hot PNS and pointed out the presence of g modes associated with its surface. Soon it was realized (Reisenegger & Goldreich 1992) that the presence of composition stratification would also allow for a different class of g modes associated with the PNS core. Having in mind the possibility of detecting the GWs associated with these oscillations, Ferrari, Miniutti & Pons (2003) computed the oscillation spectrum of sequences of equilibrium models describing the thermal evolution during the first minute of the life of a neutron star. Subsequent investigations dealt with the effects of rotation (Ferrari et al. 2004), general relativity (GR; Passamonti et al. 2005), non-linearities (Dimmelmeier, Stergioulas & Font 2006), the presence of a phase transition at the crust formation (Krüger, Ho & Andersson 2015), and a realistic EOS (Camelio et al. 2017). What these works have in common is that they use simplified background neutron star models (equilibrium configurations) and that they focus on the cooling phase spanning from the supernova explosion to the formation of the crust. Sotani & Takiwaki (2016) were the first to consider the oscillation modes of a PNS before the explosion. Their background models were based on simplified fits to numerical simulation results and they focused on the high-frequency modes (p and w modes), not considering g modes. Moreover, they considered the PNS as an isolated object, neglecting the presence of a stalled shock. One of the modes that attracted more interest is the f mode (see e.g. Lau, Leung & Lin 2010; Doneva & Kokkotas 2015), especially because of the possibility of inferring information about the rotation rate in fast-rotating progenitors (Abdikamalov et al. 2014; Fuller et al. 2015; Richers et al. 2017). However, we focus in this work on developing methods valid for the non-rotating case, which should represent the vast majority of CCSN events (Woosley & Bloom 2006). In a previous paper (Torres-Forné et al. 2018) we studied for the first time the complete eigenmode spectrum of newly-born PNSs in the Cowling approximation (including g modes) and in general relativity, using as a background the results from CCSN numerical simulations and considering the PNS–shock wave system as a whole. Comparing the eigenmode spectrum with the GW spectra from the same simulations, we showed that both were closely related and that the different p modes and g modes could be identified directly in the time–frequency distribution of the GW signal. This result showed that it is possible to perform CCSN asteroseismology and serves as a starting point to carry out inference of astrophysical parameters of PNSs using the information contained in the gravitational waveforms to be detected in future observations of current or third-generation GW detectors. We note that these results have been recently extended by Morozova et al. (2018), who included space–time perturbations in their analysis, albeit only of the lapse function, showing that space–time perturbations indeed have an impact on the calculation. In this paper, we build on the approach presented in Torres-Forné et al. (2018) and incorporate to the methodology an augmented set of space–time perturbation equations with respect to Morozova et al. (2018). The aim of this work is to understand the features observed in the GW signal of CCSN simulations, which have been interpreted as g modes of the PNS (Murphy, Ott & Burrows 2009; Cerdá-Durán et al. 2013; Müller, Janka & Marek 2013; Yakunin et al. 2015; Kuroda, Kotake & Takiwaki 2016; Andresen et al. 2017). It also attempts to shed some light on the imprint of the standing accreting shock instability (SASI; Blondin, Mezzacappa & DeMarino 2003; Foglizzo et al. 2007) in the GW signal observed in numerical simulations (Cerdá-Durán et al. 2013; Kuroda et al. 2016; Andresen et al. 2017). The paper is organized as follows: In Section 2 we develop our linear analysis formalism and estimate the GW emission for the case of a background deformed by rotation. In Section 3 we present the numerical methods used to solve the eigenvalue problem and the algorithms used in the classification of the eigenmodes. Section 4 presents our results and the comparison with the GW signal computed in the numerical simulations. Finally, in Section 5 we summarize our results and discuss their implications. The paper contains three appendices with technical details regarding numerical aspects. Throughout this paper we use a space-like metric signature (−, +, +, +) and c = G = 1 (geometrized) units, where c stands for the speed of light and G is Newton’s gravitational constant. As customary, Greek indices run from 0 to 3 and Latin indices from 1 to 3, and we use Einstein’s summation convention for repeated indices. 2 LINEAR PERTURBATIONS OF A SPHERICALLY SYMMETRIC BACKGROUND We start our analysis with the description of the perturbations of a spherically symmetric, self-gravitating equilibrium configuration. The interested reader is addressed to Kokkotas & Schmidt (1999) and Friedman & Stergioulas (2013) for detailed information on linear perturbations of compact stars and asteroseismology. Classically, this analysis was performed in Schwarzshild coordinates by Thorne & Campolattaro (1967). In our work we use isotropic coordinates instead, which are closer to the gauge condition used in the relativistic CCSN numerical simulation we employ for the comparison (Cerdá-Durán et al. 2013), which is based on the conformally flat approximation (Wilson, Mathews & Marronetti 1996; Isenberg 2008). Moreover, the derivation of the equations in these coordinates also bears resemblance to that of the equations in the Newtonian case (see Reisenegger & Goldreich 1992), which makes it easier to identify the role of the different terms in the equations and to interpret the solutions. This choice of gauge also makes it straightforward to perform the mode analysis of the Newtonian simulations we make use of (Obergaulinger, Just & Aloy 2018). The metric in a 3 + 1 foliation of the space–time in coordinates (t, xi) is \begin{eqnarray*} \mathrm{ d}s^2 = g_{\mu \nu }\mathrm{ d}x^\mu \mathrm{ d}x^\nu = (\boldsymbol{ \beta ^i}\beta _i-\alpha ^2) \mathrm{ d}t^2 + 2 \beta _i \mathrm{ d}t \mathrm{ d}x^i + \gamma _{ij} \mathrm{ d}x^i \mathrm{ d}x^j, \nonumber\\ \end{eqnarray*} (1) where βi is the shift 3-vector, α is the lapse function, and γij is the spatial 3-metric. In this work we consider the conformally flat condition approximation (CFC; Wilson et al. 1996; Isenberg 2008) of general relativity. This approximation has some advantages for simulations of CCSNe: (i) In spherical symmetry the CFC metric is exact and it is equivalent to choosing maximal slicing and isotropic coordinates. (ii) In the post-Newtonian regime (M/R < 1) it deviates from GR as the first post-Newtonian correction times the ellipticity squared (Kley & Schäfer 1999; Cordero-Carrión et al. 2009). (iii) Direct comparisons with full GR simulations of the collapse of rapidly rotating stellar cores have shown excellent agreement in both the dynamics and the GW signal (Shibata & Sekiguchi 2004; Ott et al. 2007a, b). Disagreement with GR has been estimated to be below |$1{{\ \rm per\ cent}}$| in those cases (Cerdá-Durán et al. 2005). (iv) The resulting equations can be applied directly to the analysis of numerical simulations in the CFC approximation and, with minor modifications, to Newtonian simulations. (v) In the Cowling approximation the coordinates coincide with those used in Torres-Forné et al. (2018), making possible a direct comparison with our previous results. In CFC, the spatial 3-metric is conformally flat, γij = ψfij, where ψ is the conformal factor and fij is the spatial 3-metric. The Einstein equations for the CFC metric form a purely elliptic system given by \begin{eqnarray*} \nabla ^2 \psi = -2\pi \psi ^5 \left(E + \frac{K_{ij} K^{ij}}{16 \pi } \right), \end{eqnarray*} (2) \begin{eqnarray*} \nabla ^2 Q = 2\pi Q \psi ^4 \left(E + 2 S + \frac{7 K_{ij} K^{ij}}{16 \pi } \right), \end{eqnarray*} (3) \begin{eqnarray*} \nabla ^2 \beta ^i + \frac{1}{3} \nabla ^i \nabla _j \beta ^j = 16 \pi \alpha \psi ^4 S^i + 2 \psi ^{10} K^{ij} \nabla _j \left(\frac{\alpha }{\psi ^{6}}\right) , \end{eqnarray*} (4) where ∇2 and ∇i are the Laplacian and nabla operators with respect to the flat 3-metric, respectively, Kij is the extrinsic curvature, and Q ≡ αψ. The energy–momentum content couples to the space–time geometry through the projections of the energy–momentum tensor, Tμν, on to the 3 + 1 foliation: \begin{eqnarray*} E \equiv \alpha ^2 T^{00}, \quad S_i \equiv - \alpha ^{-1}(\boldsymbol{ T}_{\mathbf{ 0}\boldsymbol{ i}} - \boldsymbol{ T}_{\boldsymbol{ ij}} \beta ^j), \end{eqnarray*} (5) \begin{eqnarray*} S_{ij} \equiv \boldsymbol{ T}_{\boldsymbol{ ij}}, \quad S \equiv S_{ij} \gamma ^{ij}. \end{eqnarray*} (6) We consider a perfect fluid, for which the energy–momentum tensor is given by \begin{eqnarray*} \boldsymbol{ T}^{\boldsymbol{ \mu \nu} } = \rho h \boldsymbol{ u}^{\boldsymbol{ \mu} } u^{\nu} + P g^{\mu \nu } \, , \end{eqnarray*} (7) where ρ is the rest-mass density, P is the pressure, uμ is the 4-velocity, h ≡ 1 + ε + P/ρ is the specific enthalpy, and ε is the specific internal energy. It is useful to define the energy density as e ≡ ρ(1 + ε). The mass conservation (continuity) and momentum conservation equations in GR read (Banyuls et al. 1997) \begin{eqnarray*} \frac{1}{\sqrt{\gamma }} \partial _t \left[ \sqrt{\gamma }D \right] + \frac{1}{\sqrt{\gamma }} \partial _i \left[ \sqrt{\gamma } D v^{*i} \right] = 0, \end{eqnarray*} (8) \begin{eqnarray*} \frac{1}{\sqrt{\gamma }} \partial _t \left[ \sqrt{\gamma } S_j \right] + \frac{1}{\sqrt{\gamma }} \partial _i \left[ \sqrt{\gamma } S_j \boldsymbol{ v^{*i}} \right] + \alpha \partial _j P = \frac{\alpha \rho h}{2} \boldsymbol{ u}^{\boldsymbol{ \mu }}u^{\nu } \partial _j g_{\mu \nu }, \nonumber\\ \end{eqnarray*} (9) where γ is the determinant of the three-metric. The conserved quantities are D = ρW and Sj = ρhW2|$v$|j, where |$W=1/\sqrt{1-v_iv^i}$| is the Lorentz factor. The Eulerian and ‘advective’ velocities are, respectively, \begin{eqnarray*} \boldsymbol{ v^i} = \frac{u^i}{W} +\frac{\boldsymbol{ \beta ^i}}{\alpha },\quad \quad \boldsymbol{ v^{*i}} = \frac{u^i}{u^0} =\alpha v^i - \boldsymbol{ \beta ^i}. \end{eqnarray*} (10) Let us consider a solution of the hydrodynamics equations that is in equilibrium (⁠|$\partial _t=0$|⁠) and is static (⁠|$v$|i = 0). In this case the metric equations (equations 2–4) and the hydrodynamics equations (equations 8–9) reduce to \begin{eqnarray*} \nabla ^2 \psi = -2\pi \psi ^5 E, \end{eqnarray*} (11) \begin{eqnarray*} \nabla ^2 Q = 2\pi Q \psi ^4 (E+2S) , \end{eqnarray*} (12) \begin{eqnarray*} \frac{1}{\rho h} \partial _i P = -\partial _i \ln \alpha \equiv \mathcal {G}_i, \end{eqnarray*} (13) and βi = 0. In the Newtonian limit |$\mathcal {G}_i$| is the gravitational acceleration, whose only non-zero component is |$\mathcal {G}_r\equiv \mathcal {G}$|⁠. The solution of equations (11)–(13) corresponds to the unperturbed state or background solution. Following Torres-Forné et al. (2018), we consider linear adiabatic perturbations of the hydrodynamics equations with respect to the background equilibrium configuration. We use the same notation and denote the Eulerian perturbations of the different quantities with δ (e.g. δρ) and without it for the background quantities (e.g. ρ). Our previous work was based on the Cowling approximation. Here, we consider the general case including also the metric perturbations, namely δα, δψ, and δβj. For a weak gravitational field (M/R ≪ 1), the leading term in the post-Newtonian expansion of the metric contributes only to α, while the first correction to ψ and βiappears at the first post-Newtonian level (Blanchet, Damour & Schaefer 1990). Therefore, for the mildly relativistic gravitational field of a PNS, we expect δα to be dominant in front of δψ and δβi. In this work we consider δβi = 0, which simplifies the equations significantly. Although this approach does not fully include the metric perturbations, it allows us to assess the effect of 1PN corrections with respect to a ‘Newtonian-like’ approach in which only δα is considered, as in Morozova et al. (2018). We denote as ξi the Lagrangian displacement of a fluid element with respect to its position at rest. Its value is related to the advective velocity as \begin{eqnarray*} \partial _t \xi ^i = \delta \boldsymbol{ v^{*i}}. \end{eqnarray*} (14) The Lagrangian perturbation of any quantity, e.g. ρ, is related to the Eulerian perturbations as \begin{eqnarray*} \Delta \rho = \delta \rho + \xi ^i\partial _i \rho . \end{eqnarray*} (15) The linearized version of equations (8) and (9) are \begin{eqnarray*} \frac{ \Delta \rho }{\rho } = -\! \left(\partial _i \xi ^i + \xi ^i \partial _i \ln \sqrt{\gamma } \right) - 6 \frac{\delta \psi }{\psi } \, , \end{eqnarray*} (16) \begin{eqnarray*} \rho h \, \partial _t \delta v_i + \alpha \partial _i \delta P = - \delta \left(\rho h \right) \partial _i \alpha - \rho h (\partial _i \delta \alpha - \delta \alpha \partial _i \ln \alpha) \, . \end{eqnarray*} (17) We use spherical coordinates, {r, θ, φ}, in which |$\sqrt{\gamma } = \psi ^6 r^2 \sin \theta$|⁠. The condition of adiabaticity of the perturbations implies that \begin{eqnarray*} \frac{\Delta P}{\Delta \rho } = \left. \frac{\partial P}{ \partial \rho } \right|_{\rm adiabatic} = h c^2_{{\rm s}} = \frac{P}{\rho }\Gamma _1, \end{eqnarray*} (18) where cs is the relativistic speed of sound and Γ1 is the adiabatic index. This allows us to write \begin{eqnarray*} \delta \left(\rho h \right) = \left(1 + \frac{1}{c_\mathrm{ s}^2} \right) \delta P - \rho h \, \xi ^i \mathcal {B}_{i} \, , \end{eqnarray*} (19) where \begin{eqnarray*} \mathcal {B}_i \equiv \frac{\partial _i e}{\rho h} - \frac{1}{\Gamma _1}\frac{\partial _i P}{P} \end{eqnarray*} (20) is the relativistic version of the Schwarzschild discriminant. Since the background is spherically symmetric, the only non-zero component is |${\mathcal {B}}_r\equiv {\mathcal {B}}$|⁠. The radial and angular parts of equation (17) are given by \begin{eqnarray*} \rho h \, \psi ^4 \alpha ^{-2} \frac{ \partial ^2 \xi ^r }{\partial t^2 } + \partial _r \delta P = \delta \left(\rho h \right) \mathcal {G}- \rho h \alpha ^{-1} (\partial _r \delta \alpha + \delta \alpha \mathcal {G}), \end{eqnarray*} (21) \begin{eqnarray*} \rho h \, \psi ^4 \alpha ^{-2} r^2 \frac{ \partial ^2 \xi ^{\theta } }{\partial t^2 } + \partial _{\theta } \delta P = - \rho h \alpha ^{-1} \partial _\theta \delta \alpha , \end{eqnarray*} (22) \begin{eqnarray*} \rho h \, \psi ^4 \alpha ^{-2} r^2 \sin ^2\theta \frac{ \partial ^2 \xi ^{\varphi } }{\partial t^2 } + \partial _{\varphi } \delta P = - \rho h \alpha ^{-1} \partial _\varphi \delta \alpha , \end{eqnarray*} (23) where we have used the fact that, in the coordinate basis, the covariant components of the velocity are given by δ|$v$|r = ψ4δ|$v$|r, δ|$v$|θ = r2ψ4δ|$v$|θ and δ|$v$|φ = r2sin2θψ4δ|$v$|φ. Additionally we need equations for the metric perturbations |$\delta \hat{Q}$| and |$\delta \hat{\psi }$|⁠, that can be obtained by linearising equations (2) and (3) and subtracting the background solution (equations 11 and 12), \begin{eqnarray*} \nabla ^2 \delta \psi = -10\pi e \psi ^4 \delta \psi - 2 \pi \psi ^5 \left[ \frac{\delta P}{c_s^2} - \rho h \xi ^r \mathcal {B} \right], \end{eqnarray*} (24) \begin{eqnarray*} \nabla ^2 \delta Q &=& 2 \pi (\rho h + 5 P) \psi ^4 (\delta Q + 4 \alpha \delta \psi) \nonumber \\ &&+ \,2\pi \alpha \psi ^5 \left[ \left(6 + \frac{2}{c_s^2} \right) \delta P - \rho h \xi ^r \mathcal {B} \right], \end{eqnarray*} (25) where we have used that δKij = 0 for δβi = 0. We perform an expansion of the perturbations with a harmonic time dependence of frequency σ and a spherical-harmonic expansion for the angular dependence \begin{eqnarray*} \delta P = \delta \hat{P} \, \, Y_{lm} e^{- i \sigma t}, \end{eqnarray*} (26) \begin{eqnarray*} {\bf \xi }^r = \eta _r \, \, Y_{lm} e^{- i \sigma t}, \end{eqnarray*} (27) \begin{eqnarray*} {\bf \xi }^\theta = \eta _\perp \, \, \frac{1}{r^{2}} \partial _\theta Y_{lm} e^{- i \sigma t}, \end{eqnarray*} (28) \begin{eqnarray*} {\bf \xi }^\varphi = \eta _\perp \, \, \frac{1}{r^2 \sin ^2\theta } \partial _\varphi Y_{lm} e^{- i \sigma t}, \end{eqnarray*} (29) and all scalar quantities (e.g. δQ and δψ) in a equivalent way to δP. The quantities ηr and η⊥ and the scalar perturbations with the hat (e.g. |$\delta \hat{P}$|⁠, |$\delta \hat{Q}$| and |$\delta \hat{\psi }$|⁠), depend on the radial coordinate only. 2.1 Non-radial perturbations (l ≠ 0) For l ≠ 0, by inserting the spherical-harmonic expansion into equations (21) and (22) we obtain \begin{eqnarray*} -\! \sigma ^2 q\, \eta _r + \partial _r \delta \hat{P} = \delta \widehat{\rho h} \mathcal {G}- \rho h \alpha ^{-1} (\partial _r \delta \hat{\alpha }+ \mathcal {G}\delta \hat{\alpha }), \end{eqnarray*} (30) \begin{eqnarray*} -\! \sigma ^2 q \, \eta _{\perp } + \delta \hat{P} = -\rho h \alpha ^{-1}\delta \hat{\alpha }, \end{eqnarray*} (31) where for convenience we have defined |$q \equiv \rho h \, \alpha ^{-2} \psi ^4$|⁠. From equation (31) it follows that \begin{eqnarray*} \delta \hat{P} = q \sigma ^2 \eta _\perp \, -\frac{\rho h}{\alpha } \delta \hat{\alpha }= q \sigma ^2 \eta _\perp \, +\rho h \left(\frac{\delta \hat{\psi }}{\psi } -\frac{\delta \hat{Q}}{Q} \right). \end{eqnarray*} (32) Using equations (32) and (19) to simplify equations (16) and (30) we obtain \begin{eqnarray*} &&{\partial _r \eta _{r} + \left[ \frac{2}{r} + \frac{1}{\Gamma _1}\frac{ \partial _r P }{P } + 6 \frac{ \partial _r \psi }{\psi } \right] \eta _r + \frac{\psi ^4 }{\alpha ^2 c^2_\mathrm{ s}} \left(\sigma ^2 - \mathcal {L}^2 \right) \eta _\perp} \nonumber \\ &&{\quad = \frac{1}{c_\mathrm{ s}^2} \frac{\delta \hat{Q}}{Q} - \left(6 + \frac{1}{c_\mathrm{ s}^2} \right)\frac{\delta \hat{\psi }}{\psi } \, ,} \end{eqnarray*} (33) \begin{eqnarray*} &&{\partial _r \eta _\perp - \left(1 - \frac{ \mathcal {N}\,^2 }{\sigma ^2} \right) \eta _r + \left[ \partial _r \ln q - \mathcal {G}\left(1 + \frac{1}{c^2_\mathrm{ s}} \right) \right] \eta _\perp} \nonumber \\ &&{\quad = \frac{\alpha ^2}{\psi ^4\sigma ^2} \left[ \partial _r (\ln \rho h) - \left(1 + \frac{1}{c_\mathrm{ s}^2} \right) \mathcal {G}\right] \left(\frac{\delta \hat{Q}}{Q} - \frac{\delta \hat{\psi }}{\psi }\right) \, , } \end{eqnarray*} (34) where |$\mathcal {L}$| is the relativistic Lamb frequency defined as \begin{eqnarray*} \mathcal {L}^2 \equiv \frac{ \alpha ^2 }{\psi ^4} c^2_\mathrm{ s} \frac{l\left(l+1\right)}{r^2} \, \end{eqnarray*} (35) and |$\mathcal {N}$| is the relativistic Brunt–Väisälä frequency \begin{eqnarray*} \mathcal {N}^2 \equiv \frac{ \alpha ^2 }{\psi ^4 } \mathcal {G}^i \mathcal {B}_i = \frac{\alpha ^2}{\psi ^{4}} \mathcal {B} \mathcal {G}. \end{eqnarray*} (36) Equations (24) and (25) for the metric perturbations can also be simplified using equation (32). After the spherical-harmonic expansion they read \begin{eqnarray*} \hat{\nabla }^2 \delta \hat{\psi } &=& -2\pi \psi ^5 \left[ \left(5 e + \frac{\rho h}{c_\mathrm{ s}^2} \right) \frac{\delta \hat{\psi }}{\psi } - \frac{\rho h}{ c_\mathrm{ s}^2} \frac{\delta \hat{Q}}{Q} \right] \nonumber \\ &&- \,2 \pi \rho h \psi ^5 \left(\frac{\psi ^4 \sigma ^2}{\alpha ^2 c_\mathrm{ s}^2}\eta _\perp - \mathcal {B} \eta _r \right), \end{eqnarray*} (37) \begin{eqnarray*} \hat{\nabla }^2 \delta \hat{Q} &=& 2 \pi (\rho h + 5 P) \alpha \psi ^5 \left(\frac{\delta \hat{Q}}{Q} + 4 \frac{\delta \hat{\psi }}{\psi }\right) \nonumber \\ &&+ \,2\pi \rho h \alpha \psi ^5 \left[ \left(6 + \frac{1}{c_\mathrm{ s}^2} \right) \left(\frac{\psi ^4\sigma ^2}{\alpha ^2 } \eta _\perp - \frac{\delta \hat{Q}}{Q} + \frac{\delta \hat{\psi }}{\psi } \right) - \eta _r \mathcal {B} \right], \nonumber\\ \end{eqnarray*} (38) where \begin{eqnarray*} \hat{\nabla }^2 \equiv \partial _{rr} + \frac{2}{r}\partial _r - \frac{l(l+1)}{r^2}. \end{eqnarray*} (39) The system of equations (33), (34), (37), and (38) can be reduced to a system of six first-order ODEs by introducing the variables |$K\equiv \partial _r\delta \hat{Q}$| and |$\Psi \equiv \partial _r \delta \hat{\psi }$|⁠. The resulting system of equations can be written in compact form as \begin{eqnarray*} \partial _r \boldsymbol u = {\boldsymbol{\sf A}} \boldsymbol u, \end{eqnarray*} (40) where |$\boldsymbol u \equiv (\eta _r,\eta _\perp , K,\delta \hat{Q}, \Psi ,\delta \hat{\psi })^{T}$| and |${\boldsymbol{\sf A}}$| is a 6 × 6 coefficient matrix whose components are explicitly written in Appendix A. 2.2 Radial perturbations (l = 0) In non-rotating stars, the Lagrangian displacement of radial oscillations has ξθ = ξφ = 0, i.e. η⊥ = 0. This case cannot be treated as a particular case of the general derivation for l ≠ 0 (previous section). In this case, the continuity and momentum equations are \begin{eqnarray*} \frac{ \delta \hat{\rho } }{\rho } = -\partial _r \eta _r - \left(\frac{2}{r} +\frac{\partial _r \rho }{\rho } + 6 \frac{\partial _r \psi }{\psi } \right) \eta _r - 6 \frac{\delta \hat{\psi }}{\psi }, \end{eqnarray*} (41) \begin{eqnarray*} &&{\quad \quad \quad -\sigma ^2 q \eta _r + \partial _r \delta \hat{P} = \mathcal {G}\delta \widehat{\rho h} - \rho h \alpha ^{-1} (\partial _r \delta \hat{\alpha } + \delta \hat{\alpha } \, \mathcal {G}) . } \end{eqnarray*} (42) From equations (18) and (19) the various scalar perturbations are given by \begin{eqnarray*} \delta \hat{P} = P \Gamma _1 \left[ \frac{\delta \hat{\rho }}{\rho } +\eta _r \left(\frac{\partial _r \rho }{\rho } - \frac{\partial _r P}{\Gamma _1 P}\right)\right], \end{eqnarray*} (43) \begin{eqnarray*} \delta \widehat{\rho h} = \left(1 + \frac{1}{c_\mathrm{ s}^2} \right) \delta \hat{P} - \rho h \, \eta _r \, \mathcal {B} . \end{eqnarray*} (44) Combining equations (41)–(44) we obtain the following equations for ηr and |$\delta \hat{P}$| \begin{eqnarray*} \partial _r \eta _r + \left[ \frac{2}{r} +\frac{1}{\Gamma _1} \frac{\partial _r P}{ P} + 6 \frac{\partial _r \psi }{\psi } \right] \eta _r + \frac{1}{P\Gamma _1} \delta \hat{P} = - 6 \frac{\delta \hat{\psi }}{\psi }, \end{eqnarray*} (45) \begin{eqnarray*} \partial _r \delta \hat{P} + q (\mathcal {N}^2 - \sigma ^2) \eta _r - \mathcal {G}\left(1 + \frac{1}{c_\mathrm{ s}^2} \right) \delta \hat{P} = -\rho h \alpha ^{-1} (\partial _r \delta \hat{\alpha }+ \delta \hat{\alpha } \mathcal {G}). \nonumber\\ \end{eqnarray*} (46) Using equations (24) and (25) the metric perturbations are \begin{eqnarray*} \hat{\nabla }^2 \delta \hat{\psi } = -10\pi e \psi ^4 \delta \hat{\psi }- 2 \pi \psi ^5 \left[ \frac{\delta \hat{P}}{c_\mathrm{ s}^2} - \rho h \eta _r \mathcal {B} \right], \end{eqnarray*} (47) \begin{eqnarray*} \hat{\nabla }^2 \delta \hat{Q} &=& 2 \pi (\rho h + 5 P) \psi ^4 (\delta \hat{Q} + 4 \alpha \delta \hat{\psi }) \nonumber \\ &&+ \,2\pi \alpha \psi ^5 \left[ \left(6 + \frac{2}{c_\mathrm{ s}^2} \right) \delta \hat{P} - \rho h \eta _r \mathcal {B} \right]. \end{eqnarray*} (48) Similarly to the case with l ≠ 0, the system of equations (45)–(48) can be written as a system of six first-order ODEs by using K and Ψ and cast in the same compact form as equation (40), where in this case |$\boldsymbol u \equiv (\eta _r,\delta \hat{P}, K,\delta \hat{Q}, \Psi ,\delta \hat{\psi })^{T}$| and |${\boldsymbol{\sf A}}$| is a 6 × 6 coefficient matrix whose components are explicitly given in Appendix A. 2.3 Boundary conditions As in Torres-Forné et al. (2018), we impose zero-displacement boundary conditions at the shock location \begin{eqnarray*} \xi _r|_\textrm{shock} = 0 \quad \rightarrow \quad \eta _r|_\textrm{shock} = 0, \end{eqnarray*} (49) which is a consequence of the impossibility of perturbations propagating across the shock from the subsonic to the supersonic region. At the origin (r = 0) we impose regularity, which for l ≠ 0 (see Reisenegger & Goldreich 1992) results in1 \begin{eqnarray*} \eta _r|_{r=0} = \frac{l }{r} \eta _\perp |_{r=0} \propto r^{l-1} \end{eqnarray*} (50) and for l = 0 in \begin{eqnarray*} \eta _r|_{r=0} \propto r , \end{eqnarray*} (51) \begin{eqnarray*} \delta \hat{P}|_{r=0} = - P \Gamma _1 \left[ 6 \frac{\delta \hat{\psi }}{\psi } + \frac{\eta _r}{r} \right]_{r=0}. \end{eqnarray*} (52) The latter expression has been obtained by evaluating equation (45) at r = 0 and imposing the boundary condition for ηr given by equation (51). Regarding metric perturbations, regularity at r = 0 implies that \begin{eqnarray*} K|_{r=0} = \frac{l}{r} \delta \hat{Q}|_{r=0} \propto r^{l-1}; \quad \Psi |_{r=0} = \frac{l}{r} \delta \hat{\psi }|_{r=0} \propto r^{l-1}. \end{eqnarray*} (53) As a consequence all four metric functions are zero at r = 0 for l ≥ 2. Outside the sources (ρ = 0) the metric perturbations fulfil |$\hat{\nabla }^2 \delta \hat{\psi }=0$| and |$\hat{\nabla }^2 \delta \hat{Q}=0$| and hence the solution decays as r−(l + 1). Outside the shock, density is non-zero but, since ηr = η⊥ = 0 and |$\delta \hat{Q}$| and |$\delta \hat{\psi }$| decay radially, it can be shown that, sufficiently far away from the shock (but still inside the star), the solution also decays as r−(l + 1). For simplicity and given that outside the shock there is a considerable drop in density, we impose this behaviour at the shock location. As a consequence \begin{eqnarray*} {}[K+ (l+1) \delta \hat{Q} / r ]_{\rm shock} = 0, \end{eqnarray*} (54) \begin{eqnarray*} {}[ \Psi + (l+1) \delta \hat{\psi }/ r ]_{\rm shock} = 0. \end{eqnarray*} (55) We note that these boundary conditions differ with respect to those used by Morozova et al. (2018), which are imposed at the PNS surface, assuming it is a free surface (ΔP = 0). 2.4 Gravitational wave emission 2.4.1 Spherical background Let us consider a general linear perturbation of the spherically symmetric background considered in the previous sections as a combination of eigenfunctions, which we denote hereafter as δρlm (and so forth for other perturbed quantities): \begin{eqnarray*} \delta \rho = \sum _{l=2}^{\infty }\sum _{m=-l}^{+l} \delta \rho _{lm} = \sum _{l=2}^{\infty } \sum _{m=-l}^{+l} \delta \hat{\rho }_{lm} \, Y_{lm}, \end{eqnarray*} (56) where |$\delta \hat{\rho }_{lm}$| can be computed as \begin{eqnarray*} \delta {\hat{\rho }}_{lm} \approx \rho \left(\frac{\mathcal {N}^2}{ \mathcal {G}}\eta _r^{lm} + \frac{\sigma ^2}{c_\mathrm{ s}^2}\eta ^{lm}_\perp \right). \end{eqnarray*} (57) In the Newtonian limit, the energy stored in all the (l, m) modes with a certain amplitude can be approximated as \begin{eqnarray*} E_{lm}= \frac{\sigma ^2}{2}\int _0^{r_{\rm shock}}\mathcal {E}_{lm} r^2 \mathrm{ d}r, \end{eqnarray*} (58) where |$\mathcal {E}_{lm}$| is the eigenmode energy density defined as \begin{eqnarray*} \mathcal {E}^{lm}(r) = \rho \left[{\eta _r^{lm}(r)}^2+l(l+1)\frac{{\eta _\perp ^{lm}(r)}^2}{r^2} \right] . \end{eqnarray*} (59) The (l, m)-mass multipole moment (l ≥ 2) can be defined as \begin{eqnarray*} \delta D_{lm} \equiv \int \mathrm{ d}V r^l \delta \rho \, Y^*_{lm} = \int _0^{r_{\rm shock}} r^{l+2} \delta \hat{\rho }_{lm}\, \mathrm{ d}r. \end{eqnarray*} (60) As a consequence, only the (l, m) mode contributes to δDlm. As we address it in the next section, this is a direct consequence of considering a spherically symmetric background. Following Thorne (1969), it is possible to compute the radiated power in GWs as \begin{eqnarray*} P_{lm} = \frac{G}{8\pi c^{2l+1}} \frac{(l+1)(l+2)}{(l-1)l}\left[ \frac{ 4\pi \sigma ^{l+1}}{(2l+1)!!} \right]^2 \left(|\delta D_{lm}|^2 + \frac{1}{c^2}|\delta J_{lm}|^2 \right), \end{eqnarray*} (61) where δJlm are the current multipoles and we show the factors c and G explicitly for the sake of the discussion. The dominant GW emission channel is the mass quadrupole (l = 2). The contribution by higher-order mass multipoles and current multipoles is suppressed by at least a factor 1/c2 and is not relevant for GW detection. The implication is that if the background is spherically symmetric, only l = 2 oscillation modes are relevant for GW emission, and the radiated power of the relevant modes can be computed as \begin{eqnarray*} P_{2m} = \frac{G}{c^{5}} \frac{4\pi \sigma ^6}{75} |\delta D_{2m}|^2 . \end{eqnarray*} (62) The amplitude of the GW emitted for (l, m) = (2, 0) modes can be computed as \begin{eqnarray*} h_{+} = -\frac{1}{D}\sin ^2\Theta \frac{8\pi }{5} \sigma ^2 \delta D_{20}\, , \end{eqnarray*} (63) \begin{eqnarray*} h_{\times } = 0 \, , \end{eqnarray*} (64) where D is the distance to the source and Θ is the observation angle with respect to the symmetry axis of the mode. In order to compare the power emitted by different modes with l = 2, we define the GW emission efficiency as \begin{eqnarray*} \rm{(GW\,efficiency)}= \frac{P}{E f}, \end{eqnarray*} (65) where f is the frequency of the mode. This equation gives an idea of the fraction of the mode energy radiated in GWs per oscillation cycle. 2.4.2 Deformed background One of the limitations of our linear analysis is that it is only applicable if the background is spherically symmetric. However, in the collapse of rapidly rotating cores, deformations in the star induced by rotation are known to considerably enhance the GW emission, increasing the wave amplitude by a factor 10–1000 (see e.g. Fryer & New 2011). Another motivation to consider a deformed background is the possible presence of l = 0 (quasi-radial) modes in the GW signal. Cerdá-Durán et al. (2013) noted that some of the features observed in the GW signal of their simulation could be explained by an l = 0 mode. As we show next, the presence of a deformed background allows for modes with l ≠ 2 to contribute to the dominant GW channel (l = 2). Let us consider an axisymmetric background of the form \begin{eqnarray*} \rho (r,\theta) = \rho _0 (r)+ a \tilde{\rho }_2 (r) Y_{20} (\theta), \end{eqnarray*} (66) where ρ0 is the spherically symmetric contribution to the background, |$\tilde{\rho }_2$|⁠, whose volume integral is normalized to unity, gives the radial dependence of the quadrupolar deformation, and a is a parameter controlling the amount of deformation (a = 0 corresponds to the spherically symmetric case considered above). We consider only quadrupolar deformations of the background because those produce the strongest coupling between modes with l ≠ 2 and the quadrupolar component of the GW signal. This allows us to compute the leading contribution to the signal. For simplicity, we consider a ≪ 1. In particular we assume that the quadrupolar deformations are sufficiently small to be neglected when compared to the perturbations themselves, i.e. \begin{eqnarray*} a \tilde{\rho }_2 \lt \lt \delta \rho _{2m} \lt \lt \rho _0. \end{eqnarray*} (67) This allows us to compute the leading-order contribution to the GW signal of modes with l ≠ 2, without the added complexity of the linear analysis including rotation. Under these conditions, the equations for the linear perturbations are identical to the equations for the spherical background. This means that we can use the eigenvalues and eigenfunctions already computed to estimate the leading-order contribution to the GW signal. The only necessary step is to compute the contribution of modes with l ≠ 2 to the mass quadrupole δD2m. Since a ≪ 1, we can consider that the new background can be described as a continuous deformation of the spherically symmetric background given by ρ0. This deformation can be described by a displacement vector Xi defined such that, for surfaces of constant density, it corresponds to a Lagrangian perturbation : \begin{eqnarray*} a \tilde{\rho }_2 Y_{20} + \boldsymbol{ X^i} \partial _i \rho _0 = 0. \end{eqnarray*} (68) Since the deformations are proportional to Y20, Xi is purely radial, i.e. Xr = Nr(r)Y20 and Xθ = Xφ = 0. Note that this only happens for quadrupolar deformations of the background and simplifies the analysis significantly. Substituting in equation (68) one arrives to \begin{eqnarray*} N^r (r) = -\frac{a \tilde{\rho }_2}{\partial _r \rho _0}. \end{eqnarray*} (69) Let us consider a new set of coordinates (r′, θ′, φ′) defined by \begin{eqnarray*} r^{\prime } \equiv r + X^r = r + N^r Y_{20}; \quad \theta ^{\prime } \equiv \theta; \quad \varphi ^{\prime } \equiv \varphi . \end{eqnarray*} (70) By construction, in the new coordinates ρ(r, θ) = ρ(r′); i.e., they are adapted to the deformation. Using equation (69) we can write \begin{eqnarray*} \mathrm{ d}r^{\prime } = \left(1 + \partial _r N^r Y_{20} \right) \mathrm{ d}r; \quad \mathrm{ d}\theta ^{\prime } = \mathrm{ d}\theta; \quad \mathrm{ d}\varphi ^{\prime } = \mathrm{ d}\varphi , \end{eqnarray*} (71) which allows us to write the line element as \begin{eqnarray*} \mathrm{ d}l^2 \approx \left(1 - 2 \partial _r N^r Y_{20}\right) \mathrm{ d}r^{\prime 2} + r^{\prime 2} \left(1 - \frac{2 N^r Y_{20}}{r^{\prime }} \right) \mathrm{ d} \Omega ^{\prime 2}, \end{eqnarray*} (72) where we have used the fact that Xr/r ≪ 1 and we have neglected higher-order corrections. Here dΩ′ ≡ sin θ′dθ′dφ′. Finally, the volume element in the new coordinates reads \begin{eqnarray*} \mathrm{ d}V \approx \left[ 1 - 2 \left(\partial _r N^r + \frac{N^r}{r^{\prime }} \right) Y_{20} \right] r^{\prime 2} \mathrm{ d}r^{\prime } \mathrm{ d}\Omega ^{\prime } . \end{eqnarray*} (73) Therefore, using the fact that the eigenfunctions and eigenvalues are unchanged at this level of approximation, we can compute the mass quadrupole by considering the change in the volume element in equation (60). The resulting contribution to the mass quadrupole of a mode with (l′, m′) is \begin{eqnarray*} \delta D_{2m,l^{\prime }m^{\prime }} &=& \int \mathrm{ d}V r^2 \delta \hat{\rho }\, Y_{l^{\prime }m^{\prime }} \, Y^*_{lm} \nonumber \\ &\approx& \int \left[ 1 - 2 \left(\partial _r N^r + \frac{N^r}{r^{\prime }} \right) Y_{20} \right] r^{\prime 2} \mathrm{ d}r^{\prime } \mathrm{ d}\Omega ^{\prime } \, \delta \hat{\rho }\, Y_{l^{\prime }m^{\prime }} \, Y^*_{lm}. \end{eqnarray*} (74) For oscillation modes with l′ = 2 \begin{eqnarray*} \delta D_{2m,2m^{\prime }} = \delta _{mm^{\prime }} \int _0^{r_{\rm shock}} r^{\prime 4} \, \delta \hat{\rho }_{2m^{\prime }} \mathrm{ d}r^{\prime } + \mathcal {O}\left(a \right), \end{eqnarray*} (75) which is equivalent to the case with a spherical background given by equation (60). However, for l′ ≠ 2, which gives a negligible contribution to the GW signal for the spherical background, we now obtain a non-zero contribution given by \begin{eqnarray*} \delta D_{2m,l^{\prime }m^{\prime }} &=& \int r^{\prime 4} \mathrm{ d}r^{\prime } \mathrm{ d}\Omega ^{\prime } \left(2 \partial _r N^r + 4 \frac{N^r}{r^{\prime }} \right) \delta \hat{\rho }_{l^{\prime }m^{\prime }} Y_{20} Y_{l^{\prime }m^{\prime }} Y^*_{2m} \nonumber \\ &&+ \,\mathcal {O}\left(a^2 \right), \end{eqnarray*} (76) The angular part of the integral can be easily computed (Arfken & Weber 1995), and the only non-zero contributions are \begin{eqnarray*} \delta D_{20,00} = \delta \hat{D}_{00}, \end{eqnarray*} (77) \begin{eqnarray*} \delta D_{20,40} = \frac{6}{7}\delta \hat{D}_{40}, \end{eqnarray*} (78) \begin{eqnarray*} \delta D_{2\pm 1,4\pm 1} = \frac{\sqrt{30}}{7}\delta \hat{D}_{4\mp 1}, \end{eqnarray*} (79) \begin{eqnarray*} \delta D_{2\pm 2,4\pm 2} = \frac{\sqrt{15}}{7}\delta \hat{D}_{4\mp 2}, \end{eqnarray*} (80) where \begin{eqnarray*} \delta \hat{D}_{l^{\prime }m^{\prime }} \equiv \frac{1}{\sqrt{4\pi }}\int _0^{r_{\rm shock}} r^{\prime 4} \left(2 \partial _r N^r + 4 \frac{N^r}{r^{\prime }} \right) \delta \hat{\rho }_{l^{\prime }m^{\prime }} \mathrm{ d}r^{\prime }. \end{eqnarray*} (81) Note that these integrals are independent of m′, except for the mode amplitude, which can be different for each m′. We can extract some conclusions from these results: (i) The only oscillation modes contributing to the quadrupolar GW emission are those with l′ = 0, 2, 4, and (ii) modes with m′ only produce GW in the (2, m′) channel. We will focus next on the contribution by δD20,00 and δD20,40, since the remaining contributions are just proportional to the latter. In order to use equations (77)–(81) to compute the GW power of l = 2, 4 modes in a deformed background, one needs to compute |$\tilde{\rho }_2$| from the original multidimensional simulation. To simplify this process and develop our intuition about the meaning of the corrections due to the deformation we will make some assumptions. Let us consider the deformation at a fixed radius r. The polar and equatorial radius of the density isocontour corresponding to ρ0(r) are given by: \begin{eqnarray*} r_{\rm e} = r + b, \end{eqnarray*} (82) \begin{eqnarray*} r_{\rm p} = r - 2b, \end{eqnarray*} (83) \begin{eqnarray*} b \equiv - \sqrt{\frac{5}{16 \pi }}N^r. \end{eqnarray*} (84) The ellipticity of the system at this radius, considering an oblate form (rp < re), is \begin{eqnarray*} e \equiv \sqrt{1-\frac{r^2_{\rm p}}{r^2_{\rm e}}}= \frac{\sqrt{3 b (2 r - b)}}{r+b}. \end{eqnarray*} (85) To simplify further, we consider constant ellipticity at different radii. From equations (84) and (85) this implies that \begin{eqnarray*} N^r = - \sqrt{\frac{16 \pi }{5}} \frac{ (3 - e^2) - 3 \sqrt{1-e^2}}{3+e^2} r, \end{eqnarray*} (86) such that for e = 0 (spherical) it results in Nr = 0. Using this expression, the integral needed to compute the contribution to the mass quadrupole is \begin{eqnarray*} \delta \hat{D}_{l^{\prime }m^{\prime }} \equiv -\sqrt{\frac{4}{5}} \frac{ (3 - e^2) - 3 \sqrt{1-e^2}}{3+e^2} \int _0^{r_{\rm shock}} r^{\prime 4} \delta \hat{\rho }_{l^{\prime }m^{\prime }} \mathrm{ d}r^{\prime }. \end{eqnarray*} (87) Note that the integral is the same as for l = 2 modes, but using the corresponding multipole. 3 EIGENMODE CALCULATION 3.1 Background models To test the capabilities of linear perturbation analysis we compute the eigenmodes for two CCSN simulations, for which we have all the quantities necessary for the analysis and the GW signal, which is used for comparison. Model s20 is the result of a simulation of the core collapse of a star of |$20 \, \mathrm{M}_{\odot }$| and of solar metallicity (Woosley & Heger 2007) with a non-zero, but dynamically negligible, rotational velocity and a magnetic field that was studied by Obergaulinger et al. (2018). The simulation was performed in axisymmetry using the ALCAR/Aenus code (Just, Obergaulinger & Janka 2015) combining special relativistic magnetohydrodynamics, an approximately general relativistic gravitational potential (version ‘A’ of the TOV potential of Marek et al. 2006), the SFHo EOS (Steiner, Hempel & Fischer 2013), and an energy-dependent, two-moment neutrino transport scheme. The core does not launch a supernova explosion. The mass of the PNS increases continuously due to the ongoing mass accretion, but does not exceed the threshold for collapse to a black hole during the first second after core bounce. Model 35OC is a 2D core-collapse simulation performed by Cerdá-Durán et al. (2013) using the general relativistic code CoCoNuT (Dimmelmeier, Font & Müller 2002; Dimmelmeier et al. 2005). The progenitor is a low-metallicity 35M⊙ star at zero-age main sequence from Woosley & Heger (2006). This progenitor has a high rotation rate and is usually regarded as a progenitor of long-duration gamma-ray bursts (GRBs). The simulation used the LS220 EOS of Lattimer & Swesty (1991) to describe matter at high density along with a simplified leakage scheme to approximate neutrino transport. This model does not produce a supernova explosion. Instead, the PNS becomes unstable to radial perturbations and collapses to a black hole 1.6 s after bounce. Fig. 1 shows the time evolution of the location of the shock, the PNS surface, and the outer boundary of the inner core for both simulations. We define the shock as the location where the flow becomes subsonic, the PNS surface as the radius at which ρ = 1011 g cm−3, and the inner core as the region where Γ1 > 2, which is a good tracer of the region where the transition to nuclear matter has occurred. Fig. 2 shows the time evolution of the mass enclosed by different radii. Note that the mass between the PNS surface and the shock is a small fraction of the PNS mass. Figure 1. View largeDownload slide Time evolution of the location of the shock (blue), the surface of the PNS (orange), and the inner core (green), for models s20 (left-hand panel) and 35OC (right-hand panel). Figure 1. View largeDownload slide Time evolution of the location of the shock (blue), the surface of the PNS (orange), and the inner core (green), for models s20 (left-hand panel) and 35OC (right-hand panel). Figure 2. View largeDownload slide Time evolution of the mass within the shock (solid blue line, left scale), the surface of the PNS (solid orange line), and the inner core (solid green line), together with the central density (blue dashed line, right scale) for models s20 (left-hand panel) and 35OC (right-hand panel). Figure 2. View largeDownload slide Time evolution of the mass within the shock (solid blue line, left scale), the surface of the PNS (solid orange line), and the inner core (solid green line), together with the central density (blue dashed line, right scale) for models s20 (left-hand panel) and 35OC (right-hand panel). 3.2 Numerical integration For the numerical calculation of the eigenmodes (‘modes’ hereafter) we follow a procedure similar to that of Torres-Forné et al. (2018); i.e., we solve the perturbation equations for different values of σ and then search for the values of σ with vanishing radial displacement at the shock (ηr|shock = 0), by means of a bisection algorithm. The perturbation equations differ depending on whetherl = 0 or l ≠ 0, but in both cases the system can be cast as in equation (40). We integrate this system of six coupled ODEs outwards from the centre to the shock radius. The integration is performed using a second-order implicit method (trapezoidal rule). Numerical tests for this procedure can be found in Torres-Forné et al. (2018). Given that we are integrating a system of six ODEs, we have to provide six boundary conditions. Since we are using a staggered grid, our first integration point is not at r = 0. This implies that we have to provide non-zero values for all quantities at this point to perform the integration. Due to the regularity conditions seen in Section 2.3, it is sufficient to fix the value of three quantities, ηr, |$\delta \hat{Q}$|⁠, and |$\delta \hat{\psi }$|⁠, and the remaining three are automatically known. The value of ηr can be set arbitrarily and it fixes the amplitude of the eigenfunction. The two additional variables have to be fixed such that the boundary conditions at the shock location, equations (54) and (55), are fulfilled. To ensure this, we use a shooting method, which comprises varying the values of |$\delta \hat{Q}$| and |$\delta \hat{\psi }$| at the innermost radial point, integrating outwards, and checking whether the boundary conditions are fulfilled. We perform this iteration using a vectorial Newton–Raphson method, where the derivatives of the Jacobian are computed numerically using a stencil of the size given by the previous step. We also compute the eigenmodes using an alternative numerical method. The description of this method and the comparison with the first one can be found in Appendix B. In all tested cases, the differences in the eigenfrequencies between both methods were smaller than |$0.1{{\ \rm per\ cent}}$|⁠. In some cases (for very low frequencies) the alternative method shows numerical convergence problems that are not found in the first method. Therefore, all results presented in this paper have been obtained using the first method, which is more robust. The computation of |$\mathcal {G}$|⁠, defined in equation (13), involves some degree of arbitrariness because it can be computed either from the gradient of the pressure (⁠|${\mathcal {G}}_P \equiv \partial _r P / \rho h$|⁠) or from the gradient of the lapse (⁠|${\mathcal {G}}_\alpha \equiv -\partial _r \ln \alpha$|⁠). Unless stated otherwise, we use |$\mathcal {G} = \mathcal {G}_P$|⁠. We explore the effect of the definition of |$\mathcal {G}$| in the eigenmode calculation in Section 4.1. 3.3 Eigenmode classification The procedure we have just described allows us to compute the eigenvalues, along with their corresponding eigenfunctions, for a set of time slices of a simulation. However, it does not provide information about the nature of each eigenmode. Stellar oscillations can be classified according to the dominant restoring force giving rise to them, either pressure (p modes) or buoyancy (g modes). The local quantities determining the character of the modes are the Lamb frequency, |$\mathcal {L}$|⁠, and the Brunt–Väisälä frequency, |$\mathcal {N}$|⁠. Pressure-supported (sound) waves with frequency σ are possible in regions of the star in which |$\sigma ^2\, \gt\, \mathcal {L}^2, \mathcal {N}^2$|⁠, while buoyancy-supported (gravity) waves are possible in regions with |$\sigma ^2\, \lt\, \mathcal {L}^2 \mathcal {N}^2$| (see e.g. Cox 1980). The regions of the star where σ2 is between |$\mathcal {L}^2$| and |$\mathcal {N}^2$| are evanescent and no waves propagate in this region. Note that if |$\mathcal {N}^2\, \lt\, 0$| gravity waves are not possible because the system is convectively unstable. We want to stress that any classification scheme tries to label the modes according to their dominant nature. However, each mode extends through regions in the star with a local character (either pressure or buoyancy dominated) that can change from point to point. Therefore, many of the modes have a mixed p/g nature. Our aim is to find this dominant character. Using the properties described above several classification procedures are possible, which have been developed in the context of asteroseismology (see e.g. Unno et al. 1979; Cox 1980). 3.3.1 Cowling classification The first classification of non-radial oscillation modes of spherical stars was introduced by Cowling (1941).2 For stars with monotonically decreasing |$\mathcal {L}^2$| (note that it is proportional to r−2) and monotonically increasing |$\mathcal {N}^2$| (typical of simple stratified equilibrium models; see e.g. Cox 1980) there is a critical frequency above which only sound waves can propagate and below which only gravity waves are possible. This allows for a very simple classification purely based on the number of nodes of the radial part of the eigenfunction, ηr. The mode with zero radial nodes is the fundamental mode or f mode, denoted as lf. Modes with higher frequencies are p modes, denoted as lpn, with an increasing number of nodes, n, for increasing frequency. In a similar way, g modes, denoted as lgn, have frequencies lower than the f mode and have an increasing number of nodes, n, for decreasing frequency. A variant of this classification has been used in Torres-Forné et al. (2018). In that work, h modes (hybrid) were introduced to distinguish modes above or below the f mode, but with the same number of nodes. The upper panels of Fig. 4 below show the eigenmodes classified according to the Cowling classification (variant described in Torres-Forné et al. 2018). Although it can serve as a guide, the Cowling classification does not work properly in every case, in particular in those cases with no critical frequency. In those cases, it may happen that, for a given frequency, there are regions of the star supporting gravity waves at the same time as other regions support sound waves, sometimes separated by evanescent regions (see upper panel of Fig. 3). In those cases the ordering devised by Cowling may not apply and one has to rely upon a more general procedure (see Unno et al. 1979 and Cox 1980 for a deeper discussion). An indication that this is indeed a problem in our models is the necessity of introducing h modes in the classification. As we show below, this is an artefact of the classification scheme. Figure 3. View largeDownload slide Propagation and phase diagrams for l = 2 modes in model 35OC at 1.3 s post bounce. Upper panel: propagation diagram showing the radial profile of the Brunt–Väisälä frequency (⁠|$\mathcal {N}^2$|⁠, solid black line) and the Lamb frequency (⁠|$\mathcal {L}^2$|⁠, dashed black line). Note that in some regions |$\mathcal {N}^2\,\lt\, 0$| and no gravity waves are possible. Colours indicate the acoustic wave region (blue, |$\sigma ^2\, \gt \,\mathcal {N}^2, \mathcal {L}^2$|⁠), the gravity wave region (red, |$\sigma ^2\, \lt\, \mathcal {N}^2, \mathcal {L}^2$|⁠), and the evanescent region (white). Black dots indicate the location of the radial nodes of ηr for all the eigenmodes found at different σ. Lower panels: phase diagram for two eigenmodes at 3542 Hz (σ2 = 5.0 × 108 s−2) corresponding to a p mode (lower left) and at 237.8 Hz (σ2 = 2.2 × 106 s−2) corresponding to a g mode (lower right). Arrows indicate the direction of increasing r. The trajectory rotates clockwise in gravity wave regions and counterclockwise in sound wave regions. Figure 3. View largeDownload slide Propagation and phase diagrams for l = 2 modes in model 35OC at 1.3 s post bounce. Upper panel: propagation diagram showing the radial profile of the Brunt–Väisälä frequency (⁠|$\mathcal {N}^2$|⁠, solid black line) and the Lamb frequency (⁠|$\mathcal {L}^2$|⁠, dashed black line). Note that in some regions |$\mathcal {N}^2\,\lt\, 0$| and no gravity waves are possible. Colours indicate the acoustic wave region (blue, |$\sigma ^2\, \gt \,\mathcal {N}^2, \mathcal {L}^2$|⁠), the gravity wave region (red, |$\sigma ^2\, \lt\, \mathcal {N}^2, \mathcal {L}^2$|⁠), and the evanescent region (white). Black dots indicate the location of the radial nodes of ηr for all the eigenmodes found at different σ. Lower panels: phase diagram for two eigenmodes at 3542 Hz (σ2 = 5.0 × 108 s−2) corresponding to a p mode (lower left) and at 237.8 Hz (σ2 = 2.2 × 106 s−2) corresponding to a g mode (lower right). Arrows indicate the direction of increasing r. The trajectory rotates clockwise in gravity wave regions and counterclockwise in sound wave regions. 3.3.2 ESO classification To overcome the problems of the Cowling classification Eckart (1960), Scuflaire (1974) and Osaki (1975) developed a classification scheme (ESO scheme, hereafter) based not only on the number of radial nodes, but also on the character of each node. Using the radius r as a parameter, each mode can be plotted in the ηr versus η⊥ phase diagram. Nodes in the radial direction correspond to crossings of the ‘x’-axis (ηr). For g modes, the trajectory of the mode in the phase diagram (parametrized with r) is clockwise (see lower right-hand panel of Fig. 3 for an example), while p modes have counterclockwise trajectories (lower left-hand panel of Fig. 3). For frequencies in regions supporting gravity waves and regions supporting sound waves, the trajectory in the phase diagram is clockwise and counterclockwise in each of the regions. The ESO scheme is based on the number of clockwise turns minus the number of counterclockwise turns in the phase diagram, nESO. For nESO > 0 the mode is predominantly a g mode and it is classified as a lgn, with n = nESO. For nESO < 0 the mode is predominantly a p mode and it is classified as a lpn with n = −nESO. The mode with nESO = 0 is the f mode. The second row of panels of Fig. 4 shows the ESO classification method applied to our models. Figure 4. View largeDownload slide Evolution of the eigenfrequencies for l = 2 for models s20 (left-hand panels) and 35OC (right-hand panels) using the Cowling (upper panels), ESO (second row), and matching (third row) classification procedures. The bottom panels correspond to the decoupled computation of the modes. A selection of the classified modes (indicated in the legends) is plotted in colours and the rest of the modes are plotted in grey. Figure 4. View largeDownload slide Evolution of the eigenfrequencies for l = 2 for models s20 (left-hand panels) and 35OC (right-hand panels) using the Cowling (upper panels), ESO (second row), and matching (third row) classification procedures. The bottom panels correspond to the decoupled computation of the modes. A selection of the classified modes (indicated in the legends) is plotted in colours and the rest of the modes are plotted in grey. The ESO scheme significantly improves the classification of the modes, although it still has some drawbacks. One of the problems is the existence of trapped modes. If one looks at the propagation diagram in the upper panel of Fig. 3, at mid-range frequencies (σ2 ∼ 5 × 108 s−2) there is a region within 10 km of the centre where only gravity waves are possible and at the same time sound waves are only possible above ∼12 km. Both regions are disconnected by an evanescent region, so in principle it is possible to have trapped modes inside each of the two regions, corresponding to a g mode and a p mode, respectively, at the same (or similar) frequencies. In practice these two modes interact with each other channelling through the evanescent region and giving rise to a more complex mode which is a hybridization of both. This effect can be observed in both of our models when looking at the time evolution of the eigenfrequencies (Fig. 4; see ESO classification). As the frequency of any two modes becomes similar the phenomenon of the avoided crossing appears. To illustrate the phenomenon of the avoided crossing, let us consider the 2f and 2g1 modes at ∼750 ms in the model s20. According to the ESO classification scheme, the 2f mode has a higher frequency than the 2g1 mode during the avoided crossing (see left-hand panel in the second row of Fig 4). This produces an abrupt change of frequency near 750 ms, which may appear as an artefact of the classification scheme (it would look more natural if they crossed). Fig. 5 shows ηr for these two modes before (610 ms), during (757 ms), and after (860 ms) the crossing. Before the crossing (left-hand panel) the 2f mode appears more concentrated in the outer parts of the system (outside the PNS), while the 2g1 mode extends to the interior of the PNS. However, after the crossing (right-hand panel) the situation is reversed, appearing as if the ESO scheme would have misclassified the modes. During the avoided crossing, both modes hybridize (middle panel) and are actually quite similar (except for the number of nodes), which is the result of the process described in the previous paragraph. This phenomenon is well known in asteroseismology, and not a numerical artefact of the ESO scheme. However, the ESO classification scheme poses a problem for the purpose of this work. Our goal is to learn how the eigenmodes behave during the post-bounce evolution to try to devise ways, in the future, to infer properties of the PNS based on GW observations. For this purpose it is crucial to characterize the GW features seen in spectrograms (mostly described as rising arcs) and classify together modes with similar features (e.g. localized in the same part of the PNS), which are likely to be excited with similar energy and produce a similar GW output, during the post-bounce evolution. For this reason we need a method that is based on the similarity of the eigenfunctions and not on the number of nodes. Figure 5. View largeDownload slide Amplitude of the radial eigenfunction for the fundamental l = 2 2f mode (blue) and the 2g1 g mode (red) for the s20 model, classified according to the ESO scheme, around the time of the avoided crossing (t − tbounce ∼ 750 ms). Figure 5. View largeDownload slide Amplitude of the radial eigenfunction for the fundamental l = 2 2f mode (blue) and the 2g1 g mode (red) for the s20 model, classified according to the ESO scheme, around the time of the avoided crossing (t − tbounce ∼ 750 ms). 3.3.3 Matching classification In this work, we present a new classification procedure that is not based on the number of nodes but on the shape of the eigenfunctions. Our procedure traces the eigenmodes in time by finding the best match between the shape of the eigenfunction in each time step and those from the previous time-steps. We have found that this procedure works best when the matching is done backwards in time. The details of the matching algorithm can be found in Appendix C. Hereafter, we refer to this algorithm as the matching classification. Note that the algorithm does not give a proper classification in the sense that it does not tell g modes from p modes or the f mode; it just groups modes at different times in groups according to similarity. To tag each mode sequence with an appropriate class we use as a reference the ESO classification at the last time available, with some additional modifications that we discuss next. The third row panel of Fig. 4 shows the classification results for the matching scheme. With the new procedure the behaviour in time of the eigenmodes is smoother and there are mode crossings at places where avoided crossings appeared with the previous two classification methods. In particular, in the crossing of modes in the s20 model discussed above, the shapes are preserved (modes are swapped with respect to Fig. 5), but the number of nodes changes. This feature is consistent for all modes and a more detailed analysis of the eigenmode shape is given in the next section. Note that in the new classification scheme, the number of nodes indicated in the name of the class (e.g. 2g1, one node) is not indicative of the actual number of nodes. This new mode classification method is not perfect, and some modes are still clearly misclassified (specially high-order g modes). However, for low-order modes it gives consistent results. We will argue in the next sections that these modes are the most relevant for GW emission, and therefore our classification should be sufficient for this purpose. To better understand the classification of the modes we solve the eigenvalue problem decoupling p modes and g modes. This can be achieved by setting |$\mathcal {B}=0$| or |$c_\mathrm{ s}^2 \rightarrow \infty$|⁠, respectively, as described in Torres-Forné et al. (2018). The two lower panels of Fig. 4 show the decoupled modes in these limits. We use this information to retag some of the modes classified by the matching algorithm to better match the identification made in the decoupled computation. This result also confirms that the avoided crossings seen in the Cowling and ESO classifications are related to crossings of the decoupled modes and therefore our matching algorithm is unveiling these crossings properly. We note that many of the decoupled modes also differ significantly in frequency with respect to the corresponding modes computed with the full system. In particular, g modes tend to have higher frequencies in the full system, likely due to the presence of acoustic wave regions where those waves can propagate significantly faster than gravity waves. One of the consequences of the matching classification scheme is that it renders it unnecessary to introduce h modes to classify all modes. For example, what was misclassified as an 2h1 mode in the Cowling and ESO classification for the 35OC model, as well as in Torres-Forné et al. (2018), is actually classified as the f mode with the new scheme and all p modes have a value of n displaced in one unit with respect to our previous classifications. A more definitive proof that our matching procedure classifies the f mode and the p modes correctly is that the frequency of the 2pn modes is approximately an integer number (n + 1) the frequency of the f mode (See Fig. 6). This relation is significantly better than the one found in Torres-Forné et al. (2018), which misclassified the f mode. It also clarifies the intriguing feature found in Torres-Forné et al. (2018) of an h-mode with a frequency that was an integer fraction of the higher-order p modes. This mode was simply the f mode. Figure 6. View largeDownload slide Evolution of the frequencies of the 2pn modes divided by n + 1 compared to the 2f mode, for both the s20 model (left) and the 35OC model (right). Figure 6. View largeDownload slide Evolution of the frequencies of the 2pn modes divided by n + 1 compared to the 2f mode, for both the s20 model (left) and the 35OC model (right). Our classification scheme is however not perfect. At early times in the simulation, when g and p modes cross frequently, it is likely that some of the modes may have been misclassified (see e.g. the behaviour of f and p modes close to bounce in model s20). We note that the idea we follow in this work of using eigenfunction similarity to classify modes is not new and similar methods have been used e.g. in the context of rotating stars (Clement 1986; Yoshida & Eriguchi 2001). Finally, we would like to indicate that there are other classification schemes in the literature (e.g. Takata 2012) that tried to overcome the limitations of the ESO scheme using different approaches. How these classification schemes compare to our approach is something that could be explored in future work. 3.4 Eigenmode morphology Figs 7–9 show the time evolution of the Newtonian eigenmode energy density, defined by equation (59), for a selection of modes classified in three groups according to their shape: f mode and p modes: Fig. 7 shows a selection of modes (2f, 2p1, 2p2, 2p3) with a relatively large amount of energy density outside the PNS. These modes are basically trapped sound waves in the region between the PNS surface and the shock. However, due to the complex hybridization process mentioned above, they also couple with regions deep in the PNS interior, specially at frequencies where mode crossings appear. For higher values of n the number of nodes (here appearing as deeps in |$\mathcal {E}$|⁠) increases as expected for higher-order overtones of the f mode. What all these modes have in common is that their frequencies increase in time (except for some periods in the s20 model) and that the frequencies are integer multiples of the f-mode frequency. As mentioned in the previous section, some of the modes have been misclassified close to bounce. This is apparent e.g. in the 2p1 mode of the 35OC model. Core g-modes: Fig. 8 collects modes with the largest amplitude inside the PNS (2g1, 2g2, and 2g3). While lower n modes are more extended over the whole PNS, higher-order modes become more concentrated in the inner core region. These modes correspond to g modes associated with the stable region in the innermost part of the PNS (see Fig. 10). The frequency evolution of these modes is similar for both models. While the 2g1 mode shows an almost monotonic increase in frequency, higher-order overtones (2g2 and 2g3) decrease their frequency with time. Surface g modes: Fig. 9 collects modes with the largest amplitude at the surface of the PNS (e.g. 2g4). These modes result from the excitation of the buoyantly stable layer (⁠|$\mathcal {N}^2\gt 0$|⁠) at the surface of the PNS (see Fig. 10). These modes show an almost monotonic increase of their frequency during their evolution, although they are confined to low frequencies (<250 Hz). Figure 7. View largeDownload slide Time evolution of the radial profile of the logarithm of the energy density, |$\mathcal {E}(r)$|⁠, for a selection of modes classified as p modes and the f mode for models s20 (left-hand panels) and 35OC (right-hand panels). The solid black line indicates the position of the PNS surface and the dashed black line the surface of the inner core. The radius is normalized to the shock location. White stripes indicate times in which no modes were found. Figure 7. View largeDownload slide Time evolution of the radial profile of the logarithm of the energy density, |$\mathcal {E}(r)$|⁠, for a selection of modes classified as p modes and the f mode for models s20 (left-hand panels) and 35OC (right-hand panels). The solid black line indicates the position of the PNS surface and the dashed black line the surface of the inner core. The radius is normalized to the shock location. White stripes indicate times in which no modes were found. Figure 8. View largeDownload slide Same as Fig. 7, but for a selection of modes classified as core g-modes. Figure 8. View largeDownload slide Same as Fig. 7, but for a selection of modes classified as core g-modes. Figure 9. View largeDownload slide Same as Fig. 7, but for a selection of modes classified as surface g modes. Note that in model 35OC some modes have been misclassified by our algorithm and appear as sharp transitions. Figure 9. View largeDownload slide Same as Fig. 7, but for a selection of modes classified as surface g modes. Note that in model 35OC some modes have been misclassified by our algorithm and appear as sharp transitions. Figure 10. View largeDownload slide Logarithm of the Brunt–Väisälä frequency, |$\mathcal {N}$|⁠, for positive values of |$\mathcal {N}^2$|⁠, i.e. regions stable against convection. Convectively unstable regions (⁠|$\mathcal {N}^2\,\lt\, 0$|⁠) are plotted in black. The blue-black line indicates the position of the PNS surface and the dashed blue line the surface of the inner core. The radius is normalized to the shock location. Figure 10. View largeDownload slide Logarithm of the Brunt–Väisälä frequency, |$\mathcal {N}$|⁠, for positive values of |$\mathcal {N}^2$|⁠, i.e. regions stable against convection. Convectively unstable regions (⁠|$\mathcal {N}^2\,\lt\, 0$|⁠) are plotted in black. The blue-black line indicates the position of the PNS surface and the dashed blue line the surface of the inner core. The radius is normalized to the shock location. This interpretation is consistent with previous analysis (albeit more simplified) of the PNS structure (Cerdá-Durán et al. 2013; Andresen et al. 2017), which showed that the features observed in the GW spectrograms could be related to g modes in two different regions (surface and inner core). 4 COMPARISON WITH SIMULATIONS 4.1 GW spectrograms In this section we compare the eigenmode frequencies obtained applying the linear analysis to the spherical background with the GW frequencies computed in the numerical simulations, a result of the multidimensional dynamics of the model. Fig. 11 shows the GW signal (upper panels), the corresponding spectrograms (middle panels), and the spectrograms with a selection of modes overplotted (lower panels). As mentioned before, there is some uncertainty in the definition of |$\mathcal {G}$|⁠. To gauge the difference between using |$\mathcal {G}_P$| and |$\mathcal {G}_\alpha$|⁠, we plot in all cases the modes computed in both ways. The range between both options gives a measure of the error introduced by the definition of |$\mathcal {G}$|⁠. Figure 11. View largeDownload slide This figure shows the GW signal (upper panels), the corresponding spectrograms (middle panels), and the spectrograms with a selection of modes overplotted (lower panels) for models s20 (left) and 35OC (right). Solid lines and dashed lines are used to indicate that the calculations were made using |$\mathcal {G}_P$| and |$\mathcal {G}_\alpha$|⁠, respectively. Note that for model s20 those two lines overlap. Figure 11. View largeDownload slide This figure shows the GW signal (upper panels), the corresponding spectrograms (middle panels), and the spectrograms with a selection of modes overplotted (lower panels) for models s20 (left) and 35OC (right). Solid lines and dashed lines are used to indicate that the calculations were made using |$\mathcal {G}_P$| and |$\mathcal {G}_\alpha$|⁠, respectively. Note that for model s20 those two lines overlap. For the non-rotating s20 model (left-hand panels), both calculations lie on top of each other. This is a strong indication that the assumption of hydrostatic equilibrium inside the shock is indeed a very good approximation. However, in the rapidly rotating model 35OC (right-hand panels), there is a significant difference between both definitions of |$\mathcal {G}$| in many of the models (solid lines correspond to |$\mathcal {G}_P$| and dashed lines to |$\mathcal {G}_\alpha$|⁠). The difference may be due to the fact that we are neglecting rotation in our linear analysis. In fact, centrifugal forces play an important role in the equilibrium for this model and, thus, this result is not a surprise. What is reassuring (and somewhat striking) is that, even using a non-rotating analysis, the features observed in the spectrogram lie between our two possibilities for the eigenmode computation. In general, for both models (see lower panels of Fig. 11) there is a clear agreement of a selection of modes with the features in the spectrograms, which is a clear indication that these modes are responsible for the computed GW emission. In model s20 (bottom left-hand panel) the agreement for modes 2g1 and 2g2 is remarkable, with some hints of the presence of the 2f mode. Although p modes are not clearly visible, some excess power above the 2g1 mode may be an indication of their presence. In model 35OC (bottom right-hand panel), all the main features can be explained with a few modes. Most of them match well the spectrogram within the error produced by the definition of |$\mathcal {G}$|⁠. The only exception is the 2f mode, whose evolution runs parallel to the spectrogram feature but with a higher frequency. The main features can be explained by the 2g1 mode and the 2p1 mode. The f mode and all p modes up to order 5 are also clearly visible, albeit with lower amplitudes. We note in particular that our computation of the l = 0 mode is able to reproduce the characteristic feature of this mode close to black hole formation, namely that its frequency goes to zero at the onset of instability (Cerdá-Durán et al. 2013), as predicted by Chandrasekhar (1964). In addition to estimating the effect of the definition of |$\mathcal {G}$| in our mode comparison, we also test its effect on the expression for the Brunt–Väisälä frequency. In this work we first perform an angular average of the simulation data and then we compute the Brunt–Väisälä frequency as |$\mathcal {N}^2 = \mathcal {G} \mathcal {B}$|⁠, |$\mathcal {G}$| and |$\mathcal {B}$| being the radial component of the vectors |$\boldsymbol{ \mathcal {G}_i}$| and |$\boldsymbol{ \mathcal {B}_i}$|⁠. Alternatively one can compute |$\mathcal {N}^2 = \boldsymbol{ \mathcal {G}_i} \mathcal {B}^i$|⁠, on the 2D grid of the simulation, and then perform the angular average to obtain 1D profiles of |$\mathcal {N}^2$|⁠. For the fast-rotating case, the second procedure takes into account the non-radial components of |$\boldsymbol{ \mathcal {G}_i}$| and |$\boldsymbol{ \mathcal {B}_i}$|⁠, which are otherwise neglected in the first procedure. We have computed the eigenmodes using both definitions and the differences in the computed eigenfrequencies do not differ by more than |$1{{\ \rm per\ cent}}$|⁠. Regarding gravity, this work improves over previous work that considered no metric perturbations (Cowling, Torres-Forné et al. 2018) and only perturbations of the lapse function (Morozova et al. 2018). To compare with different approaches for the metric perturbations we have performed our analysis in Cowling and considering only perturbations of |$\delta \hat{Q}$|⁠, which corresponds to including only perturbations of the lapse function. Fig. 12 compares the three approaches. In general, the differences among them increase at later times. This is because the PNS becomes gradually more compact and GR effects become more relevant. Note in particular the 2g1 mode for the s20 model (yellow lines in the left-hand panel), in which only the approach followed in this work is able to reproduce the features observed in the spectrogram. For the rotating model 35OC, this is more difficult to discuss, because the differences with respect to the spectrogram are also strongly influenced by rotation. However, it is worth highlighting that our approach is the only one capable of reproducing the turning down of the 0f mode at the onset of black hole formation, while the other two approaches give qualitatively different behaviours. Therefore, we conclude that our approach is necessary whenever high accuracy in the eigenmode calculations is needed, specially close to black hole formation. Moreover, the high accuracy obtained in model s20 provides convincing evidence that perturbations of the shift vector, neglected in our analysis, are not important, at least for slowly rotating systems. Figure 12. View largeDownload slide Spectrograms with a selection of modes overplotted using space–time perturbations (solid lines), only metric perturbations of |$\delta \hat{Q}$| (dashed lines), and no metric perturbations (Cowling, dotted lines), for models s20 (left) and 35OC (right). Figure 12. View largeDownload slide Spectrograms with a selection of modes overplotted using space–time perturbations (solid lines), only metric perturbations of |$\delta \hat{Q}$| (dashed lines), and no metric perturbations (Cowling, dotted lines), for models s20 (left) and 35OC (right). 4.2 GW efficiency and mode energy To understand which modes are more efficient GW emitters, we compute the GW efficiency (see equation 65 for l = 2 and l = 0 modes). This calculation provides an estimation of the fraction of the energy stored in an eigenmode that is emitted per oscillation cycle. Model s20 is non-rotating, so the l = 0 mode will not contribute to the GW signal. Model 35OC corresponds to a rapidly rotating progenitor, which results in a PNS with a central period of about 1.5 ms at bounce. The PNS is considerably deformed, with a polar-to-equatorial radius ratio of 0.64 after bounce, and the GW amplitude is expected to be significantly affected by rotation (see Section 2.4.2). We use equation (87) to estimate the quadrupolar moment needed for the computation of the GW efficiency. In our numerical simulation the ratio re/rp ranges from ∼1.25 (in the PNS core) to ∼1.7 (at the PNS surface), which corresponds to ellipticities in the range e ∼ 0.6–0.8. The pre-factor to the integral of equation (87) is in the interval −0.06 to −0.14 for our ellipticity range. This means that, for similar eigenfunction structure and mode energy, the l = 0 modes are expected to emit GWs with an amplitude of about |$10{{\ \rm per\ cent}}$| the amplitude of the l = 2 modes. Fig. 13 shows the GW efficiency for both models. Similarly to the results of Torres-Forné et al. (2018), the efficiency grows with the frequency, with the p modes being in general more efficient than the g modes. This happens because the GW efficiency approximately scales with σ3. The implication is that to see low-frequency modes in the spectrogram one needs considerably large energy in those modes. In both models the most visible features correspond to the lowest-order modes, while high-order p modes are only observed in the 35OC model, likely related to the SASI activity (see the related discussion in Section 5). In neither model are high-order g modes observed, the highest-order mode observed being the 2g2 of the s20 model. Figure 13. View largeDownload slide GW efficiency for the l = 2 modes in model s20 (upper panel) and for the l = 2 and l = 0 modes in model 35OC (middle panel and lower panel, respectively). Figure 13. View largeDownload slide GW efficiency for the l = 2 modes in model s20 (upper panel) and for the l = 2 and l = 0 modes in model 35OC (middle panel and lower panel, respectively). Given that we know the complete eigenmode structure, we can use the spectrogram to estimate the energy stored in each mode and try to validate the previous claims. To this end, we extract the amplitude of the GW emission at some intervals along the main features of the spectrogram. Then, we rescale the amplitude of the eigenmodes so that their GW emission amplitude computed with equation (64) matches the value obtained from the spectrogram. The results of this calculation for model s20 are shown in the upper panel of Fig. 14. The modes that correspond to the main features on the spectrogram (Fig. 11, middle left-hand panel) are the 2g1 and the 2g2. The highest power in GW corresponds to the 2g1 mode, while the other mode presents a slightly lower amplitude. However, the energy stored in each mode (Fig. 14) exhibits the opposite behaviour. This is due to the GW efficiency (Fig. 13, upper panel). As the efficiency of the 2g2 mode is much lower than that of the 2g1 mode, the corresponding energy should be larger to radiate a similar energy in GW. In the case of model 35OC (Fig. 14, lower panel), we have analysed the four modes that show the most clear trace in the spectrogram (Fig. 11 middle right-hand panel). The energies are ordered as expected. The two fundamental modes, 0f and 2f, have the largest energies because their GW amplitudes are large (their traces on the spectrogram are clearly visible) but their efficiencies are low. Figure 14. View largeDownload slide Energy of the principal emission modes in logarithmic scale for l = 2 modes in model s20 (upper panel) and for l = 2 and l = 0 modes in model 35OC (lower panel). Figure 14. View largeDownload slide Energy of the principal emission modes in logarithmic scale for l = 2 modes in model s20 (upper panel) and for l = 2 and l = 0 modes in model 35OC (lower panel). 5 SUMMARY AND DISCUSSION In this work we have developed a numerical procedure to solve the eigenvalue problem of hydrodynamic and metric perturbations of a spherically symmetric self-gravitating system in general relativity. We have applied this method to compute the oscillation frequencies of the PNS–shock system formed during a core-collapse supernova explosion. Those frequencies have been compared to the time–frequency patterns observed in the GW templates from two CCSN numerical simulations. This work is an extension of our previous investigation (Torres-Forné et al. 2018) and brings forth significant improvements with respect to previous results. The numerical routines developed and used in this work to solve the eigenvalue problem are available at the GREAT (General Relativistic Eigenmode Analysis Tool) library that we have released as open source (https://www.uv.es/cerdupa/codes/GREAT/). The main results of this work can be summarized as follows: We have incorporated perturbations of the lapse function and of the conformal factor, improving the approach of Torres-Forné et al. (2018), based on the Cowling approximation, and of Morozova et al. (2018), which considered only perturbations of the lapse function. Our results show that it is necessary to consider both kinds of perturbations to accurately reproduce the l = 2 modes in the GW spectrograms. Regarding l = 0 modes, our approach is the only one able to trace the qualitative time–frequency behaviour of the 0f mode, for which the frequency decreases towards zero at the onset of black hole formation. This mode is of particular importance because it indicates the formation of a black hole and could be observed for rapidly rotating cores. Our approach has not accounted for perturbations of the shift vector. However, given the accuracy of our results, incorporating those perturbations does not seem necessary, at least for non-rotating models. For our non-rotating model s20, we have found an excellent agreement between the features observed in the spectrogram and the computed eigenmodes. For rapidly rotating cores, the assumed spherically symmetric background is not valid and the results are only qualitatively similar to the GW signal, due to the effect of centrifugal forces, not considered in our analysis. In both CCSN models, s20 and 35OC, the 2g1 mode has been identified as the main contributor to the GW signal. The 2f mode is also visible in the two models, along with a few overtones. Furthermore, we have estimated the eigenmode energy according to the amplitude of the GW signal for each mode. This analysis confirms that most of the energy is stored in the lowest-order eigenmodes. We have developed a formalism to estimate the contribution of quasi-radial modes (l = 0) to the quadrupolar component of the GW signal in the case of deformed backgrounds. This is of particular importance for rapidly rotating cores in which the PNS is deformed by the effect of centrifugal forces. Our improved analysis was possible thanks to a newly developed matching classification scheme, which allowed us to classify the eigenmodes in p, g, and f modes in a more accurate way than the preceding methods used in Torres-Forné et al. (2018). Despite our investigation being limited to only two numerical models, some conclusions about the mechanism producing GWs during the collapse of massive stars can be extracted from this work, in particular when comparing with previous results in the literature. The presence of a clear pattern of rising frequencies in the GW spectrograms of the post-bounce evolution of PNSs has been observed in a number of works (Murphy et al. 2009; Cerdá-Durán et al. 2013; Müller et al. 2013; Yakunin et al. 2015; Kuroda et al. 2016; Andresen et al. 2017). In all these works, this raising pattern was attributed to g modes of the PNS. Our work shows that this main feature is probably the 2g1 mode, i.e. the lowest-order g mode, associated with the buoyantly stable region in the innermost part of the PNS, and not to its surface. Although surface g modes are possible according to our analysis, we do not see their correspondence in the GW spectrograms of the two models computed. This is very likely due to the low GW efficiency of these modes, which displace a much smaller amount of matter when compared with the core g modes. Some works in the literature (Cerdá-Durán et al. 2013; Kuroda et al. 2016; Andresen et al. 2017) have related the presence of a low-frequency (∼100 Hz) component in the GW signal to the characteristic frequency of the SASI. In this work, however, we have identified the mode observed in Cerdá-Durán et al. (2013) as the fundamental 2f mode, and we suspect a similar result will hold for the results of Kuroda et al. (2016) and Andresen et al. (2017). The fundamental mode (and also higher-order p modes) seems to be excited in cases with strong SASI activity. Cerdá-Durán et al. (2013) showed that these features observed in the GW spectrogram match perfectly with the features observed in the spectrogram of the time evolution of the shock; i.e., the shock oscillates with frequencies matching the p modes. It is not completely clear why this is the case (see discussion in Cerdá-Durán et al. 2013). In principle, in the presence of the SASI the shock oscillates with frequencies corresponding to the unstable modes of the advective-acoustic cycle coupling the shock and the PNS surface (Foglizzo 2002; Foglizzo et al. 2007), which should not coincide with the frequencies of the purely acoustic cycle (p modes in our case). We are not sure if the matching of the frequencies observed in our analysis is generic or if it is a particular feature of the 35OC model (e.g. a resonance). The analysis presented in this work can be applied to other mutidimensional simulations presenting signatures of SASI to try to better understand the relationship between p modes and SASI. Our clear identification of the 2g1 and 2f modes is of great importance to devise strategies to infer PNS properties from a possible GW observation of a nearby supernova (or black hole formation) event. It follows from our work that most of the future efforts must focus on these two particular modes, in order to learn how they depend on the properties of the PNS. This will also allow development of data analysis methods aimed at detecting this kind of GW signals buried in detector noise. Our analysis may also simplify the GW templates from CCSNe, since most of the information (the time–frequency behaviour) encodes the general evolution of the PNS, and not the particular non-linear perturbations, which are to a large degree stochastic. Some of the information cannot be extracted from our analysis, most notably the energy stored in each of the eigenmodes, which does not allow us to provide complete templates of CCSNe. For this we still will have to rely on sophisticated core-collapse simulations. One of the main questions that have to be addressed is what is the evolution of the energy content of each eigenmode. The results presented here are only for two 2D simulations, but the energy distribution may change significantly in 3D simulations or in simulations accounting for more sophisticated microphysics or neutrino transport. Finally, we note that the results presented in this work strongly rely on our new eigenmode classification algorithm. Previous results by Cerdá-Durán et al. (2013) suffered from serious misclassification issues. In particular, the so-called SASI mode was not classified as the fundamental mode, which gave rise to confusion in explaining why it was an integer division of the higher-order modes. This is however explained naturally with its identification as an f-mode, possible with the new matching classification scheme. The new method still suffers from some problems. First, it sometimes misclassifies the modes at early times near core bounce. Secondly, the method is not completely automatic and requires some degree of manual intervention. We hope to improve our classification method in the future to allow for a fully automatic and robust eigenmode classification procedure. ACKNOWLEDGEMENTS Work supported by the Spanish MINECO (grant AYA2015-66899-C2-1-P), by the Generalitat Valenciana (PROMETEOII-2014-069), and by the European Gravitational Observatory (EGO-DIR-51-2017). P.C acknowledges the support from the Ramon y Cajal programme of the Spanish MINECO (RYC-2015-19074). A.P. acknowledges support from the European Union under the Marie Sklodowska Curie Actions Individual Fellowship, grant agreement no 656370. J.A.F. acknowledges support from the European Union’s Horizon 2020 RISE programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740. Footnotes 1 Note that in Torres-Forné et al. (2018) there is a missing r factor in the text (but not in the calculations). 2 We warn the reader not to confuse the Cowling approximation (static space–time approximation) and the Cowling classification procedure of modes, both introduced in Cowling (1941). 3 Linear Algebra PACKage library, http://www.netlib.org/lapack/. REFERENCES Abbott B. P. et al. , 2017a , Nature , 551 , 85 https://doi.org/10.1038/nature24471 Crossref Search ADS Abbott B. P. et al. , 2017b , ApJ , 848 , L12 https://doi.org/10.3847/2041-8213/aa91c9 Crossref Search ADS Abbott B. P. et al. , 2017c , ApJ , 848 , L13 https://doi.org/10.3847/2041-8213/aa920c Crossref Search ADS Abbott B. P. et al. , 2017d , ApJ , 850 , L39 https://doi.org/10.3847/2041-8213/aa9478 Crossref Search ADS Abbott B. P. et al. , 2018 , Phys. Rev. Lett. , 121 , 161101 Crossref Search ADS PubMed Abdikamalov E. , Gossan S. , DeMaio A. M. , Ott C. D. , 2014 , Phys. Rev. D , 90 , 044001 https://doi.org/10.1103/PhysRevD.90.044001 Crossref Search ADS Adsuara J. E. , Cordero-Carrión I. , Cerdá-Durán P. , Aloy M. A. , 2016 , J. Comput. Phys. , 321 , 369 https://doi.org/10.1016/j.jcp.2016.05.053 Crossref Search ADS Andresen H. , Müller B. , Müller E. , Janka H.-T. , 2017 , MNRAS , 468 , 2032 https://doi.org/10.1093/mnras/stx618 Crossref Search ADS Annala E. , Gorda T. , Kurkela A. , Vuorinen A. , 2018 , Phys. Rev. Lett. , 120 , 172703 https://doi.org/10.1103/PhysRevLett.120.172703 Crossref Search ADS PubMed Arfken G. B. , Weber H. J. , 1995 , Mathematical Methods for Physicists , 4th edn . Academic Press , San Diego, CA Banyuls F. , Font J. A. , Ibáñez J. M. , Martí J. M. , Miralles J. A. , 1997 , ApJ , 476 , 221 https://doi.org/10.1086/303604 Crossref Search ADS Blanchet L. , Damour T. , Schaefer G. , 1990 , MNRAS , 242 , 289 https://doi.org/10.1093/mnras/242.3.289 Crossref Search ADS Blondin J. M. , Mezzacappa A. , DeMarino C. , 2003 , ApJ , 584 , 971 https://doi.org/10.1086/345812 Crossref Search ADS Camelio G. , Lovato A. , Gualtieri L. , Benhar O. , Pons J. A. , Ferrari V. , 2017 , Phys. Rev. D , 96 , 043015 Crossref Search ADS Cerdá-Durán P. , Faye G. , Dimmelmeier H. , Font J. A. , Ibáñez J. M. , Müller E. , Schäfer G. , 2005 , A&A , 439 , 1033 https://doi.org/10.1051/0004-6361:20042602 Crossref Search ADS Cerdá-Durán P. , DeBrye N. , Aloy M. A. , Font J. A. , Obergaulinger M. , 2013 , ApJ , 779 , L18 https://doi.org/10.1088/2041-8205/779/2/L18 Crossref Search ADS Chandrasekhar S. , 1964 , ApJ , 140 , 417 https://doi.org/10.1086/147938 Crossref Search ADS Clement M. J. , 1986 , ApJ , 301 , 185 https://doi.org/10.1086/163886 Crossref Search ADS Cordero-Carrión I. , Cerdá-Durán P. , Dimmelmeier H. , Jaramillo J. L. , Novak J. , Gourgoulhon E. , 2009 , Phys. Rev. D , 79 , 024017 https://doi.org/10.1103/PhysRevD.79.024017 Crossref Search ADS Cowling T. G. , 1941 , MNRAS , 101 , 367 https://doi.org/10.1093/mnras/101.8.367 Crossref Search ADS Cox J. P. , 1980 , Theory of Stellar Pulsation , Princeton University Press , NJ De S. , Finstad D. , Lattimer J. M. , Brown D. A. , Berger E. , Biwer C. M. , 2018 , Phys. Rev. Lett. , 121 , 091102 Crossref Search ADS PubMed Dimmelmeier H. , Font J. A. , Müller E. , 2002 , A&A , 388 , 917 https://doi.org/10.1051/0004-6361:20020563 Crossref Search ADS Dimmelmeier H. , Novak J. , Font J. A. , Ibáñez J. M. , Müller E. , 2005 , Phys. Rev. D , 71 , 064023 https://doi.org/10.1103/PhysRevD.71.064023 Crossref Search ADS Dimmelmeier H. , Stergioulas N. , Font J. A. , 2006 , MNRAS , 368 , 1609 https://doi.org/10.1111/j.1365-2966.2006.10274.x Crossref Search ADS Doneva D. D. , Kokkotas K. D. , 2015 , Phys. Rev. D , 92 , 124004 https://doi.org/10.1103/PhysRevD.92.124004 Crossref Search ADS Eckart C. , 1960 , Hydrodynamics of Oceans and Atmospheres , Pergamon Press , New York Fattoyev F. J. , Piekarewicz J. , Horowitz C. J. , 2018 , Phys. Rev. Lett. , 120 , 172702 https://doi.org/10.1103/PhysRevLett.120.172702 Crossref Search ADS PubMed Ferrari V. , Miniutti G. , Pons J. A. , 2003 , MNRAS , 342 , 629 https://doi.org/10.1046/j.1365-8711.2003.06580.x Crossref Search ADS Ferrari V. , Gualtieri L. , Pons J. A. , Stavridis A. , 2004 , Class. Quantum Gravity , 21 , S515 https://doi.org/10.1088/0264-9381/21/5/019 Crossref Search ADS Foglizzo T. , 2002 , A&A , 392 , 353 https://doi.org/10.1051/0004-6361:20020912 Crossref Search ADS Foglizzo T. , Galletti P. , Scheck L. , Janka H.-T. , 2007 , ApJ , 654 , 1006 https://doi.org/10.1086/509612 Crossref Search ADS Friedman J. L. , Stergioulas N. , 2013 , Rotating Relativistic Stars, Cambridge Monographs on Mathematical Physics . Cambridge Univ. Press , Cambridge https://doi.org/10.1017/CBO9780511977596 Fryer C. L. , New K. C. B. , 2011 , Living Rev. Relativ. , 14 , 1 https://doi.org/10.12942/lrr-2011-1 Crossref Search ADS PubMed Fuller J. , Klion H. , Abdikamalov E. , Ott C. D. , 2015 , MNRAS , 450 , 414 https://doi.org/10.1093/mnras/stv698 Crossref Search ADS Gossan S. E. , Sutton P. , Stuver A. , Zanolin M. , Gill K. , Ott C. D. , 2016 , Phys. Rev. D , 93 , 042002 https://doi.org/10.1103/PhysRevD.93.042002 Crossref Search ADS Isenberg J. A. , 2008 , Int. J. Mod. Phys. D , 17 , 265 https://doi.org/10.1142/S0218271808011997 Crossref Search ADS Just O. , Obergaulinger M. , Janka H.-T. , 2015 , MNRAS , 453 , 3386 https://doi.org/10.1093/mnras/stv1892 Crossref Search ADS Kley W. , Schäfer G. , 1999 , Phys. Rev. D , 60 , 027501 https://doi.org/10.1103/PhysRevD.60.027501 Crossref Search ADS Klimenko S. et al. , 2016 , Phys. Rev. D , 93 , 042004 https://doi.org/10.1103/PhysRevD.93.042004 Crossref Search ADS Kokkotas K. D. , Schmidt B. G. , 1999 , Living Rev. Relativ. , 2 , 2 https://doi.org/10.12942/lrr-1999-2 Crossref Search ADS PubMed Krüger C. J. , Ho W. C. G. , Andersson N. , 2015 , Phys. Rev. D , 92 , 063009 https://doi.org/10.1103/PhysRevD.92.063009 Crossref Search ADS Kuroda T. , Kotake K. , Takiwaki T. , 2016 , ApJ , 829 , L14 https://doi.org/10.3847/2041-8205/829/1/L14 Crossref Search ADS Lattimer J. M. , Swesty F. D. , 1991 , Nucl. Phys. A , 535 , 331 https://doi.org/10.1016/0375-9474(91)90452-C Crossref Search ADS Lau H. K. , Leung P. T. , Lin L. M. , 2010 , ApJ , 714 , 1234 https://doi.org/10.1088/0004-637X/714/2/1234 Crossref Search ADS Marek A. , Dimmelmeier H. , Janka H.-T. , Müller E. , Buras R. , 2006 , A&A , 445 , 273 https://doi.org/10.1051/0004-6361:20052840 Crossref Search ADS McDermott P. N. , van Horn H. M. , Scholl J. F. , 1983 , ApJ , 268 , 837 https://doi.org/10.1086/161006 Crossref Search ADS Morozova V. , Radice D. , Burrows A. , Vartanyan D. , 2018 , ApJ , 861 , 10 Crossref Search ADS Müller B. , Janka H.-T. , Marek A. , 2013 , ApJ , 766 , 43 https://doi.org/10.1088/0004-637X/766/1/43 Crossref Search ADS Murphy J. W. , Ott C. D. , Burrows A. , 2009 , ApJ , 707 , 1173 https://doi.org/10.1088/0004-637X/707/2/1173 Crossref Search ADS , Obergaulinger M. , Just O. , Aloy M. Á. , 2018 , J. Phys. G: Nucl. Part. Phys. , 45 , 084001 https://doi.org/10.1088/1361-6471/aac982 Crossref Search ADS Osaki J. , 1975 , PASJ , 27 , 237 Ott C. D. , Dimmelmeier H. , Marek A. , Janka H.-T. , Zink B. , Hawke I. , Schnetter E. , 2007a , Class. Quantum Gravity , 24 , S139 https://doi.org/10.1088/0264-9381/24/12/S10 Crossref Search ADS Ott C. D. , Dimmelmeier H. , Marek A. , Janka H.-T. , Hawke I. , Zink B. , Schnetter E. , 2007b , Phys. Rev. Lett. , 98 , 261101 https://doi.org/10.1103/PhysRevLett.98.261101 Crossref Search ADS Passamonti A. , Bruni M. , Gualtieri L. , Sopuerta C. F. , 2005 , Phys. Rev. D , 71 , 024022 https://doi.org/10.1103/PhysRevD.71.024022 Crossref Search ADS Reisenegger A. , Goldreich P. , 1992 , ApJ , 395 , 240 https://doi.org/10.1086/171645 Crossref Search ADS Richers S. , Ott C. D. , Abdikamalov E. , O’Connor E. , Sullivan C. , 2017 , Phys. Rev. D , 95 , 063019 Crossref Search ADS Scuflaire R. , 1974 , A&A , 36 , 107 Shibata M. , Sekiguchi Y.-I. , 2004 , Phys. Rev. D , 69 , 084024 https://doi.org/10.1103/PhysRevD.69.084024 Crossref Search ADS Sotani H. , Takiwaki T. , 2016 , Phys. Rev. D , 94 , 044043 https://doi.org/10.1103/PhysRevD.94.044043 Crossref Search ADS Steiner A. W. , Hempel M. , Fischer T. , 2013 , ApJ , 774 , 17 https://doi.org/10.1088/0004-637X/774/1/17 Crossref Search ADS Takata M. , 2012 , PASJ , 64 , 66 https://doi.org/10.1093/pasj/64.4.66 Crossref Search ADS Thorne K. S. , 1969 , ApJ , 158 , 997 https://doi.org/10.1086/150259 Crossref Search ADS Thorne K. S. , Campolattaro A. , 1967 , ApJ , 149 , 591 https://doi.org/10.1086/149288 Crossref Search ADS Thrane E. , Coughlin M. , 2015 , Phys. Rev. Lett. , 115 , 181102 https://doi.org/10.1103/PhysRevLett.115.181102 Crossref Search ADS PubMed Torres-Forné A. , Cerdá-Durán P. , Passamonti A. , Font J. A. , 2018 , MNRAS , 474 , 5272 https://doi.org/10.1093/mnras/stx3067 Crossref Search ADS Unno W. , Osaki Y. , Ando H. , Shibahashi H. , 1979 , Nonradial Oscillations of Stars , University of Tokyo Press , Tokio Wilson J. R. , Mathews G. J. , Marronetti P. , 1996 , Phys. Rev. D , 54 , 1317 https://doi.org/10.1103/PhysRevD.54.1317 Crossref Search ADS Woosley S. E. , Bloom J. S. , 2006 , ARA&A , 44 , 507 https://doi.org/10.1146/annurev.astro.43.072103.150558 Crossref Search ADS Woosley S. E. , Heger A. , 2006 , ApJ , 637 , 914 https://doi.org/10.1086/498500 Crossref Search ADS Woosley S. E. , Heger A. , 2007 , Phys. Rep. , 442 , 269 https://doi.org/10.1016/j.physrep.2007.02.009 Crossref Search ADS Yakunin K. N. et al. , 2015 , Phys. Rev. D , 92 , 084040 https://doi.org/10.1103/PhysRevD.92.084040 Crossref Search ADS Yoshida S. , Eriguchi Y. , 2001 , MNRAS , 322 , 389 https://doi.org/10.1046/j.1365-8711.2001.04115.x Crossref Search ADS APPENDIX A: COEFFICIENTS OF THE PERTURBATIVE EQUATIONS The system of equations described in Sections 2.1 and 2.2 (for l ≠ 0 and l = 0, respectively) can be cast in a matrix form as \begin{eqnarray*} \partial _r \boldsymbol u = {\boldsymbol{\sf A}} \boldsymbol u. \end{eqnarray*} (A1) For l ≠ 0, |$\boldsymbol u \equiv (\eta _r,\eta _\perp , K,\delta \hat{Q}, \Psi ,\delta \hat{\psi })^{T}$| and the non-zero elements of |${\boldsymbol{\sf A}}$| are given by \begin{eqnarray*} A_{11} &=& - \frac{2}{r} - \frac{ \mathcal {G}}{c_\mathrm{ s}^2} - 6 \frac{ \partial _r \psi }{\psi }\, , \nonumber \\ A_{12} &=& \frac{\psi ^4 }{\alpha ^2 c^2_\mathrm{ s}} \left(\mathcal {L}^2 - \sigma ^2\right)\, , \nonumber \\ A_{14} &=& \frac{1}{c_\mathrm{ s}^2 Q}\, , \nonumber \\ A_{16} &=& - \left(6 + \frac{1}{c_\mathrm{ s}^2} \right)\frac{1}{\psi }\, , \end{eqnarray*} \begin{eqnarray*} A_{21} &=& \left(1 - \frac{ \mathcal {N}^2 }{\sigma ^2} \right)\, , \nonumber \\ A_{22} &=& -\partial _r \ln q + \mathcal {G}\left(1 + \frac{1}{c^2_\mathrm{ s}} \right)\, , \nonumber \\ A_{24} &=& \frac{\alpha }{\psi ^5\sigma ^2} \left[ \partial _r (\ln \rho h) - \left(1 + \frac{1}{c_\mathrm{ s}^2} \right) \mathcal {G}\right]\, , \nonumber \\ A_{26} &=& -\frac{\alpha ^2}{\psi ^5\sigma ^2} \left[ -\partial _r (\ln \rho h) + \left(1 + \frac{1}{c_\mathrm{ s}^2} \right) \mathcal {G}\right] = -\alpha A_{24}\, , \end{eqnarray*} \begin{eqnarray*} A_{31} &=& -2\pi \rho h \alpha \psi ^5 \mathcal {B}\, , \nonumber \\ A_{32} &=& 2\pi \rho h \alpha \psi ^5 \left(6 + \frac{1}{c_\mathrm{ s}^2}\right) \frac{\psi ^4 \sigma ^2}{\alpha ^2}\, , \nonumber \\ A_{33} &=& -\frac{2}{r}\, , \nonumber \\ A_{34} &=& \frac{l(l+1)}{r^2} - 2\pi \psi ^4\left(5 e + \frac{\rho h}{c_\mathrm{ s}^2}\right)\, , \nonumber \\ A_{36} &=& 2 \pi \psi ^4 \alpha \left[ 10 e + 30 P + \frac{\rho h}{c_\mathrm{ s}^2} \right]\, , \end{eqnarray*} \begin{eqnarray*} A_{43} = 1\, , \end{eqnarray*} \begin{eqnarray*} A_{51} &=& 2 \pi \rho h \psi ^5 \mathcal {B}\, , \nonumber \\ A_{52} &=& -2 \pi \rho h \psi ^5 \frac{\psi ^4 \sigma ^2}{\alpha ^2 c_\mathrm{ s}^2}\, , \nonumber \\ A_{54} &=& 2 \pi \psi ^4 \frac{\rho h}{\alpha c_\mathrm{ s}^2}\, , \nonumber \\ A_{55} &=& -\frac{2}{r}\, , \nonumber \\ A_{56} &=& \frac{l(l+1)}{r^2} - 2 \pi \psi ^4 \left(5 e + \frac{\rho h}{c_\mathrm{ s}^2} \right)\, , \nonumber \\ A_{65} &=& 1\, . \end{eqnarray*} For l = 0, |$\boldsymbol u \equiv (\eta _r,\delta \hat{P}, K,\delta \hat{Q}, \Psi ,\delta \hat{\psi })^{T}$| and the non-zero elements of |${\boldsymbol{\sf A}}$| are \begin{eqnarray*} A_{11} &=& - \frac{2}{r} - \frac{ \mathcal {G}}{c_\mathrm{ s}^2} - 6 \frac{ \partial _r \psi }{\psi }\, , \nonumber \\ A_{12} &=& - \frac{1}{\rho h c_\mathrm{ s}^2}\, , \nonumber \\ A_{16} &=& -\frac{6}{\psi }\, , \end{eqnarray*} \begin{eqnarray*} A_{21} &=& -q (\mathcal {N}^2 - \sigma ^2)\, , \nonumber \\ A_{22} &=& \mathcal {G}\left(1 + \frac{1}{c^2_\mathrm{ s}} \right)\, , \nonumber \\ A_{23} &=& -\rho h \alpha ^{-1} \psi ^{-1}\, , \nonumber \\ A_{24} &=& \rho h \alpha ^{-1} \psi ^{-1} (\partial _r \ln \psi - \mathcal {G})\, , \nonumber \\ A_{25} &=& \rho h \psi ^{-1}\, , \nonumber \\ A_{26} &=& -\rho h \psi ^{-1} (\partial _r \ln \psi)\, , \end{eqnarray*} \begin{eqnarray*} A_{31} &=& -2\pi \rho h \alpha \psi ^5 \mathcal {B}\, , \nonumber \\ A_{32} &=& 2 \pi \alpha \psi ^5 \left(6 + \frac{2}{c_\mathrm{ s}^2} \right)\, , \nonumber \\ A_{33} &=& -\frac{2}{r}\, , \nonumber \\ A_{34} &=& 2 \pi \psi ^4 (\rho h + 5 P)\, , \nonumber \\ A_{36} &=& 8 \pi \alpha \psi ^4 (\rho h + 5 P)\, , \end{eqnarray*} \begin{eqnarray*} A_{43} = 1\, , \end{eqnarray*} \begin{eqnarray*} A_{51} &=& 2 \pi \rho h \psi ^5 \mathcal {B}\, , \nonumber \\ A_{52} &=& -2 \pi \psi ^5 \frac{1}{c_\mathrm{ s}^2}\, , \nonumber \\ A_{55} &=& -\frac{2}{r}\, , \nonumber \\ A_{56} &=& -10 \pi \psi ^4 e\, , \nonumber \\ A_{65} &=& 1\, . \end{eqnarray*} APPENDIX B: ALTERNATIVE NUMERICAL METHOD We have also computed the eigenmodes using an alternative numerical method. Instead of solving the system of six coupled ODEs, we solve a system of two ODEs for the fluid variables (ηr and η⊥ for the case l ≠ 0, or ηr and |$\delta \hat{P}$| for l = 0) coupled with two elliptic equations for the metric perturbations |$\delta \hat{\psi }$| and |$\delta \hat{Q}$| (equations 37 and 38 for l ≠ 0 or equations 47 and 48 for l = 0). Each of the metric equations, as well as the corresponding boundary conditions given by equations (54) and (55), is discretized to second-order accuracy and written as two linear systems of equations that are solved using the LAPACK library.3 Details on the implementation and tests of the elliptic solver can be found in Adsuara et al. (2016). Since the metric equations and the fluid equations are coupled, we obtain the solution of the system of four equations (for each value of σ) in an iterative way. We first integrate the fluid variables considering |$\delta \hat{Q}=\delta \hat{\psi }=0$| (Cowling), then we use the values of their values to compute the metric perturbations, and then we continue with the iteration until the residual of all four quantities, computed as the L2 norm of the difference between two consecutive iterations, is below 10−4. For most of the values of σ, this procedure converges to a solution in less than ∼10 iterations. Below a certain threshold in the value of σ, the iterative procedure becomes unstable and no convergence is achieved, even using a small relaxation factor in the iteration. The reason is that towards lower values of σ, there appear g modes with an increasing number of nodes. For sufficiently low values of σ, the g modes have a number of nodes comparable to the number of grid points and this triggers point-to-point numerical noise, which spoils the solution. However, this is not a limitation of our method, because it only affects the calculation of very high-order g modes (typically above order 10), which are not relevant to GW observations (see discussion in Section 5.3 in Torres-Forné et al. 2018) and are not well resolved in simulations, anyway. We have compared the eigenmodes computed with this alternative method with the ones given in the main text. The discrepancy in the eigenfrequencies is in all cases below |$0.1{{\ \rm per\ cent}}$|⁠. Given that the alternative method becomes numerically unstable in some cases and does not improve the solution, we only present results obtained with our main method in this work. APPENDIX C: EIGENFUNCTION CLASSIFICATION PROCEDURE We describe in this appendix the steps that we carry out to classify the eigenmodes of our linear analysis. This classification is based on the similarity with the modes from previous time-steps. For each mode at each time, we interpolate linearly the values of ηr and η⊥ (and all necessary quantities) to an equally spaced grid of 300 points between the centre and the shock location. Instead of the radial coordinate we use a rescaled coordinate x ∈ [0, 1], which maps the interval r ∈ [0, rcore] to x ∈ [0, 1/3], the interval r ∈ [rcore, rpns] to x ∈ [1/3, 2/3] and the interval r ∈ [rpns, rshock] to x ∈ [2/3, 1], where rcore, rpns, and rshock are the radial location of the core surface, the PNS surface, and the shock. We normalize the eigenfunctions such that they all have the same mean energy density, E/V = 1, where E is the energy of the mode, given by equation (58), and V is the volume of the region inside the shock. In this way, the eigenfunctions at different times are easier to compare with each other. Next, we count the number of interior nodes (not counting the ones at the centre and at the shock) by searching for changes in the sign of ηr. This information can be used to classify the modes according to the Cowling or ESO classification schemes, and serves as a guide here. We compute the energy density of each mode, |$\mathcal {E}(x)$|⁠, using equation (59). Our matching algorithm is based on the similarity of this function at different times. This quantity has some advantages with respect to the eigenfunctions. First, it is a combination of the radial and perpendicular parts, and secondly it is a positive function and does not suffer from the sign ambiguity that the eigenfunctions have. Our matching algorithm serves to identify a series of modes at different times as members of the same class. However, it does not give a name for the class. We use the ESO classification at the starting point of the algorithm to give a preliminary tag for the modes. In most of the cases we applied our matching procedure backwards in time, using as starting point the last time output (the exception being the l = 0 modes). For each of the identified modes, lmn, with m = {f, p, g, h} denoting the possible mode classes, we create a template |$\mathcal {E} (r; ^l m_l)$| as a basis for comparison. We proceed to the next time output and we compare all the templates with the energy density of the new eigenmodes |$\mathcal {E}(r;\sigma)$|⁠, where σ belongs to all possible eigenvalues of the new time output. For all possible values of σ and lmn (for the same l) we compute the L2 norm of the difference \begin{eqnarray*} L2 (\sigma | ^l m_n) \equiv \sum _{i=1}^{N} \left(\mathcal {E} (x_i; ^l m_n) - \mathcal {E} (x_i; \sigma) \right)^2 \, . \end{eqnarray*} (C1) Values of L2(σ|lmn) ≪ 1 indicate a good match between the eigenfunction corresponding to σ and the template for lmn. We restrict our comparison to frequencies σ that are close to the sequence corresponding to the template lmn. To do this, we extrapolate the sequence of the modes already classified as lmn to the new time and compare this value with σ. If the relative difference is larger than a certain threshold (⁠|$10-20{{\ \rm per\ cent}}$|⁠) we reject this combination as possible. We use linear extrapolation in this procedure for model 35OC and constant extrapolation for s20. The extrapolation function is a least-squares fit to the last 10−20 points already classified. We order all possible matching combinations of modes in ascending order (better to worse matching) and assign sequentially to each unclassified mode a class given by the corresponding matching template. Modes already classified are removed from the sequence to avoid repetitions. If there are unclassified modes, we use them to create new templates using the ESO criterion. We update the templates incorporating the information of the newly classified eigenfunction and repeat the process for the next time output. For p- and f-mode sequences, the template is a mean of the last 10 classified modes. This allows for smooth variations in the form of the eigenfunction over time, in particular the location of the nodes. For g modes, which evolve more slowly in time, we use all previously classified modes. Once all modes have a preliminary classification, we perform some modifications to improve the matching: For f and p modes, we reorder the frequencies according to the number of radial nodes at each frequency such that higher-order modes have more nodes. This solves some misclassification issues of high-order p modes. We retag manually some of the low-order modes such that the time evolution of the profiles of |$\mathcal {E}$| is similar to the corresponding modes in the decoupled case (see Section 3.3.3) and that the p modes are approximately integer multiples of the f mode. Note that the classification is not fully automatic as it requires the adjustment of a few parameters and thresholds and the manual retag at the end, for each of the models. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Towards asteroseismology of core-collapse supernovae with gravitational wave observations – II. Inclusion of space–time perturbations
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty2854
DA - 2019-01-21
UR - https://www.deepdyve.com/lp/oxford-university-press/towards-asteroseismology-of-core-collapse-supernovae-with-XL1L9nwf3F
SP - 3967
VL - 482
IS - 3
DP - DeepDyve
ER -