© ESO, 2012

1. Introduction

The electromagnetic signature of core-collapse supernovae has been exploited comprehensively and thoroughly by countless observations during the past decades, providing only indirect information about the explosion mechanism, however. The up to now only recorded neutrino signal of a core-collapse supernova (SN1987A) confirmed the idea that the collapse of the core of a massive star to neutron star densities provides the necessary energy for the explosion (Baade & Zwicky 1934). Because gravitational waves (GW), the only other means to probe the supernova engine besides neutrinos, are yet to be detected, supernova modelers are preparing for this prospective measurement by predicting the gravitational wave signature of core-collapse supernovae with ever increasing realism (for reviews, see e.g., Kotake et al. 2006; Ott 2009; Fryer & New 2011).

According the Einstein’s theory of general relativity (GR), gravitational waves will be generated by any mass-energy distribution whose quadrupole (or higher) moment varies non-linearly with time. In core-collapse supernovae this criterion is satisfied by time-dependent rotational flattening particularly during collapse and bounce, prompt post-shock convection, non-radial flow inside the proto-neutron star and in the neutrino-heated hot bubble, the activity of the standing accretion shock instability (SASI), asymmetric emission of neutrinos, and by asymmetries associated with the effects of magnetic fields (for a recent review see, e.g., Ott 2009, and references therein). While significant rotational deformation and dynamically relevant magnetic fields require particular progenitors that possess some (considerable) amount of rotational and magnetic energy or that must rotate fast and differentially (additional differential rotation develops during collapse) to amplify an initially weak magnetic field, all other processes are genuinely operative in any core-collapse supernova.

Obviously, an adequate modeling of these effects and an accurate prediction of the GW signal ultimately requires three dimensional (3D) hydrodynamic simulations covering the collapse, bounce, and post-bounce evolution (at least) until a successful launch of the explosion and including a proper treatment of neutrino transport and the relevant microphysics. However, most studies of the past thirty years were either concerned with the collapse and bounce signal only (Müller 1982; Finn & Evans 1990; Mönchmeyer et al. 1991; Yamada & Sato 1994; Zwerger & Müller 1997; Rampp et al. 1998; Dimmelmeier et al. 2001, 2002; Kotake et al. 2003; Shibata 2003; Shibata & Sekiguchi 2004; Ott et al. 2004; Cerda-Duran et al. 2005; Saijo 2005; Shibata & Sekiguchi 2005; Kotake et al. 2006; Dimmelmeier et al. 2007; Ott et al. 2007; Dimmelmeier et al. 2008), or were restricted to axisymmetric (2D) models (Müller et al. 2004; Ott et al. 2006; Kotake et al. 2007; Marek et al. 2009; Murphy et al. 2009; Yakunin et al. 2010). Several authors also investigated the influence of magnetic fields on the GW signal during the collapse and early post-bounce evolution assuming axisymmetry (Kotake et al. 2004; Yamada & Sawai 2004; Kotake et al. 2005; Obergaulinger et al. 2006a,b) and no symmetry restriction at all (Scheidegger et al. 2008, 2010). The GW signal caused by aspherical neutrino emission was first studied by Epstein (1978) and subsequently by Burrows & Hayes (1996), Müller & Janka (1997), and Kotake et al. (2007, 2009a,b, 2011), where the investigations by Müller & Janka (1997) and Kotake et al. (2009b, 2011) considered also 3D, i.e., non-axisymmetric models.

Concerning 3D post-bounce models, the topic of the study presented here, Müller & Janka (1997) were the first to analyze the GW signature of 3D non-radial flow and anisotropic neutrino emission from prompt post-bounce convection in the outer layers of a proto-neutron star during the first 30 ms after supernova-shock formation. Because of smaller convective structures and slower overturn velocities, the GW wave amplitudes of their 3D models are more than a factor of 10 lower than those of the corresponding 2D ones, and the wave amplitudes associated with anisotropic neutrino emission are a factor of 10 higher than those caused by non-radial gas flow. Fryer (2004) and Fryer et al. (2004) presented gravitational wave signals obtained from 3D core-collapse simulations that extend up to 150 ms past bounce and were performed with a Newtonian smoothed particle hydrodynamics code coupled to a three-flavor gray flux-limited diffusion neutrino transport scheme. Gravitational wave emission occurs owing to matter asymmetries that arise from perturbations caused by precollapse convection, core rotation, and low-mode convection in the explosion engine itself, and owing to anisotropic neutrino emission. Kotake et al. (2009b) simulated 3D mock-up models that mimic neutrino-driven explosions aided by the SASI, and computed the GW signal resulting from anisotropic neutrino emission by means of a ray-tracing method in a post-processing step. They pointed out that the gravitational waveforms of 3D models vary much more stochastically than those of axisymmetric ones, i.e., in 3D the GW signals do not possess any template character. However, when considering rotating models, Kotake et al. (2011) argue that the GW waveforms resulting from neutrino emission exhibit a common feature, which results from an excess of neutrino emission along the spin axis due to the growth of spiral SASI modes. Scheidegger et al. (2008) simulated the collapse of two rotating and magnetized cores in 3D until several 10 ms past bounce, applying a parametrized deleptonization scheme (Liebendörfer 2005) that is a good approximation until a few milliseconds past bounce. Scheidegger et al. (2010) extended their study by systematically investigating the effects of the equation of state, the initial rotation rate, and both the initial magnetic field strength and configuration on the GW signature. They also simulated three representative models until  ~200 ms of post-bounce accretion (no development of explosions) incorporating a treatment for neutrino transport based on a partial implementation of the isotropic diffusion source approximation (Liebendörfer et al. 2009).

In the following we present the GW analysis of a set of parametrized 3D models of neutrino-powered supernova explosions covering the evolution from core bounce until  ~1.4 s later, where the high-density inner core of the proto-neutron star (PNS) is excised and replaced by a time-dependent boundary condition and a central point mass. The neutrino transport is treated by an approximate solver based on a ray-by-ray treatment of the multi-dimensional effects (Scheck et al. 2006). Because we analyze the GW signal arising from both non-radial mass flow and anisotropic neutrino emission, we discuss the neutrino emission of these 3D models as well, and particularly address its multidimensional properties, some of which have previously been investigated in 2D by Janka & Mönchmeyer (1989a,b), Ott et al. (2008), Kotake et al. (2009a), Marek & Janka (2009), Marek et al. (2009), and Brandt et al. (2011).

Based on 2D simulations of rotational core-collapse, Janka & Mönchmeyer (1989a,b) pointed out that the energy output in neutrinos seen by an observer may vary as much as a factor of 3 with his inclination angle relative to the rotation axis, while for the neutrino energy much smaller angular variations are to be expected. Marek et al. (2009) and Marek & Janka (2009) found that neutrino luminosities and mean energies show quasi-periodic time variability in their 2D simulations of non-rotating and slowly rotating 15 M_⊙ progenitors covering up to  ~700 ms past bounce. Caused by the expansion and contraction of the shock in the course of SASI oscillations, the level of the fluctuations is  ≲50% for the luminosities and roughly 1 MeV for the mean neutrino energies. The luminosity fluctuations are somewhat bigger for ν_e and  than for heavy-lepton neutrinos. The neutrino quantities also depend on polar angle as a consequence of the preference of the SASI motions along the symmetry axis of the 2D models. Additional short-wavelength spatial variations of the average radiated energies and of the neutrino fluxes are caused by local downdrafts of accreted matter. This is in accordance with the findings of Müller & Janka (1997), who inferred from a post-processing analysis of the neutrino emission anisotropy that features in the GW signal of their 2D models of convection in the hot-bubble region are well-correlated with structures in the neutrino signal. The features are associated with sinking and rising lumps of matter and with temporal variations of aspherical accretion flows toward the proto-neutron star. Kotake et al. (2009a) computed neutrino anisotropies with a ray-tracing scheme by post-processing their 2D SASI models and derived analytical expressions for evaluating GW signals for neutrinos in 3D models, too. A generalization of these expressions will be presented below. Brandt et al. (2011) performed 2D multi-group, multi-angle neutrino transport simulations for both a non-rotating and rapidly rotating 20 M_⊙ model extending  ~400 ms beyond bounce. Their simulations predict that the neutrino radiation fields vary much less with angle than the matter quantities in the region of net neutrino heating because most neutrinos are emitted from deeper radiative regions and because the neutrino energy density combines the specific intensities as integrals over sources at many angles and depths. The rapidly rotating model exhibits strong, flavor-dependent asymmetries in both peak neutrino flux and light curves, the peak flux and decline rate having pole-equator ratios  ≲3 and  ≲2, respectively. Brandt et al. (2011) also provide estimates of the detectability of neutrino fluctuations in IceCube and Super-Kamiokande as previously done by Lund et al. (2010) on the basis of the Marek et al. (2009) non-rotating models.

The paper is organized as follows: in Sect. 2 we discuss the numerical methods, the input physics, and the properties of the progenitor models and the set of 3D simulations that we analyzed. Section 3 contains a description of the formalism we used to extract the observable neutrino signal of our 3D models, and a discussion of some of its properties relevant for the corresponding GW signal. In Sect. 4 we give the formalism necessary to calculate the GW signature of 3D non-radial flow and anisotropic neutrino emission, and discuss the GW signature of the investigated 3D models. Finally, in Sect. 5 we summarize our results and discuss shortcomings and possible implications of our study.

2. Model setup

2.1. Code and computational grid

The 3D supernova models we analyzed for their neutrino and GW signature have been simulated with the explicit finite-volume, Eulerian, multi-fluid hydrodynamics code Prometheus (Fryxell et al. 1991; Müller et al. 1991a,b). This code integrates the multidimensional hydrodynamic equations using the dimensional splitting method of Strang (1968), the piecewise parabolic method of Collela & Woodward (1984), and a Riemann solver for real gases proposed by Colella & Glaz (1985). Inside grid cells with strong grid-aligned shocks fluxes computed from the Riemann solver are replaced by the AUSM+ fluxes of Liou (1996) in order to prevent odd-even decoupling (Quirk 1994). Nuclear species are advected using the consistent multi-fluid advection (CMA) scheme of Plewa & Müller (1999).

The simulation code employs an axis-free overlapping “Yin-Yang” grid (Kageyama & Sato 2004) in spherical polar coordinates, which was recently implemented into Prometheus, for spatial discretization (Wongwathanarat et al. 2010a). The Yin-Yang grid relaxes the CFL-timestep condition and avoids numerical artifacts near the polar axis. Concerning the ray-by-ray neutrino transport no special procedure needs to be applied for the Yin-Yang grid. The (scalar) quantities involved in the transport algorithm are computed on both the Yin and Yang grid and are then linearly interpolated as any other scalar quantity.

The grid consists of 400(r)  ×  47(θ)  ×  137(φ)  ×  2 cells corresponding to an angular resolution of 2° and covers the full 4π solid angle. The radial grid has an equidistant spacing of 0.3 km from the inner grid boundary out to r = 80 km (models W15 and N20; see Table 1) or 115 km (model L15; see Table 1), respectively. Beyond this radius the radial grid is logarithmically spaced. The outer radial grid boundary R_ob is at 18 000 km, which is sufficiently far out to prevent the supernova shock from leaving the computational domain during the simulated epoch. This radial resolution suffices the requirement that there are always more than 15 radial zones per decade in density.

A central region, the dense inner core of the proto-neutron star (PNS) at ρ ≳ 10^12...13 g cm^-3, is excised from the computational domain and replaced by an inner time-dependent radial boundary condition and a point mass at the coordinate origin. The radius of the inner radial boundary shrinks according to Eq. (13) of Scheck et al. (2008). For the W15 and N20 models the initial and final (asymptotic) boundary radii are  km and  km, respectively. For the L15 models the corresponding radii are 82 km and 25 km, accounting for a less extreme contraction of the neutron star within the simulation time. The timescale for the contraction is t_ib = 1 s for all models. This choice of parameters implies R_ib ≈ 19 km at 1.3 s for models W15 and N20, and R_ib ≈ 30 km at 1.4 s for the L15 models. Hydrostatic equilibrium is assumed at the inner radial grid boundary R_ib, which is thus a Lagrangian (co-moving) position, while a free outflow boundary condition is employed at the outer radial grid boundary (for more details, see Wongwathanarat 2011; Wongwathanarat et al. 2010a).

2.2. Input physics

Self-gravity is fully taken into account by solving Poisson’s equation in integral form using an expansion into spherical harmonics as in Müller & Steinmetz (1995). The monopole term of the potential is corrected for general relativistic effects as described in Scheck et al. (2006) and Arcones et al. (2007). The cooling of the PNS is prescribed by neutrino emission properties (luminosities and mean neutrino energies) as boundary condition at the inner radial grid boundary (for details, see Scheck et al. 2006). The contraction of the PNS is mimicked by the movement of the inner radial grid boundary (see Sect. 2.1). “Ray-by-ray” neutrino transport and neutrino-matter interactions are approximated as in Scheck et al. (2006) by radial integration of the one-dimensional (spherical), gray transport equation for all angular grid directions (θ, φ) independently. This approach allows for angular variations of the neutrino fluxes. The tabulated equation of state (EoS) of Janka & Müller (1996) is used to describe the stellar fluid. It includes arbitrarily degenerate and arbitrarily relativistic electrons and positrons, photons, and four predefined nuclear species (n, p, α, and a representative Fe-group nucleus) in nuclear statistical equilibrium.

2.3. Models

We have analyzed a set of 3D simulations (Wongwathanarat 2011; Wongwathanarat et al. 2010b, and in prep.) based on two 15 M_⊙ progenitor models (W15 and L15), and a 20 M_⊙ progenitor model (N20). The W15 model is obtained from the non-rotating 15 M_⊙ progenitor s15s7b2 of Woosley & Weaver (1995), the L15 model from a star evolved by Limongi et al. (2000), and the N20 model from a SN 1987A progenitor of Shigeyama & Nomoto (1990). The progenitor models were evolved through collapse to 15 ms after bounce with the Prometheus-Vertex code in one dimension (Marek & Buras, priv. comm.) providing the initial models for the 3D simulations. To break the spherical symmetry of the initial models, random seed perturbations of 0.1% amplitude are imposed on the radial velocity (v_r) field at the beginning of the 3D simulations. Explosions are initiated by neutrino heating at a rate that depends on suitable values of the neutrino luminosities imposed at the lower boundary chosen such that the desired value of the explosion energy is obtained. The evolution is followed until 1.3 ṡ after bounce for the W15 and N20 progenitor models, while the L15 models are simulated until 1.4 s post-bounce. The GW analysis presented below comprises five models (see Table 1), where models W15-2 and W15-4 differ only by the initial seed perturbations. The explosion energies, E_exp, given in Table 1 are instantaneous values at the end of the simulations (t = t_f), adding up the total energies (kinetic + internal + gravitational) in all zones where the sum of these energies is positive. The explosion time, t_exp, is defined as the time when this sum exceeds a value of 10⁴⁸ erg, roughly corresponding to the time when the average shock radius is 400 to 500 km (see, however, Pejcha & Thompson 2012, for an alternative definition of the time of the onset of the explosion).

3. Neutrino signal

The non-radial motions caused by the SASI and convection in the neutrino-heated hot-bubble as well as by convection inside the proto-neutron star (driven by Ledoux unstable lepton gradients) give rise to a time-dependent, anisotropic emission of neutrinos of all flavors, and thus to the emission of gravitational waves (Epstein 1978; Burrows & Hayes 1996; Müller & Janka 1997; Kotake et al. 2007, 2009a,b), as discussed in Sect. 4. We have analyzed this emission for the 3D models discussed above (see Sect. 2.3), particularly addressing its multidimensional properties.

3.1. Formalism

To derive observable luminosities of an emitting source we consider an observer located at far distance D from that source (see Fig. 1). According to definition the flux measured by the observer is given by the following integral at the position of the observer: (1)where μ is the cosine of the angle between the direction of the radiation and the line of sight (between the observer and the center of the source), ω denotes the radiation direction at the observer’s location (defined by a pair of angles), and dω is the solid angle element around the radiation direction ω. The intensity I adopts non-zero values within the opening angle subtended by the emitting surface (not necessarily a sphere). We note that here and in the following we suppress the dependence of the intensity on the neutrino energy and assume that energy-integrated quantities are considered (the outlined formalism, however, is valid also for an energy-dependent treatment). The integration over ω at the observer’s location can be substituted by an integration over the emitting surface of the source, because the radiation intensity is constant along rays, i.e., (2)for any ray arriving at the observer from the source (and zero otherwise), where R_o denotes the position of a surface element of the emitting surface in the coordinate frame of the source and ω_o the direction of the radiation field at that position toward the observer. Note that we ignore in Eq. (2) the trivial effect that the time t for I_o relative to the time for I is subject to a retardation.

Moreover, in the following we disregard spectral and angular corrections that may be relevant when the emitting surface is in relative motion to the observer or sitting deep in the gravitational potential of a compact star (in which case general relativistic (GR) effects like redshifting and ray bending would be important). Considering the source to be at rest is a good assumption for the neutrinospheric region in the supernova core after bounce (the velocities of mass motions in this layer are unlikely to be higher than some 1000 km s^-1, i.e., at most around one percent of the speed of light), while GR energy redshift is certainly of relevance on the  ~10−20% level during the proto-neutron star cooling phase (t ≳ 1 s after bounce), but much lower in the accretion and shock revival phases, when the forming neutron star is still considerably less compact than the final remnant.

Fig. 1 Sketch illustrating the various quantities involved when deriving the observable luminosity of a radiating source, whose visible surface is shaded in blue. When the observer is at infinity, the emitting surface is a sphere, and Ro is measured from the origin of this sphere, one obtains α = γ.

For a distant observer D ≫ max {|R_o|}  holds, i.e., the value of μ is very close to one for the whole emitting surface. Denoting the solid angle subtended by a surface element of the emitting surface by dω, we have dω = dA_⊥/D², where dA_⊥ = cosγ   dA is the projected area of this surface element perpendicular to the line of sight, when γ is the angle between the normal unit vector n_A of the emitting surface element dA and the line of sight (see Fig. 1), but specified to the case D ≫ |R_o|. Hence, we obtain for the observable luminosity the expression (3)In order to evaluate the integral in Eq. (3), one needs to know the intensity I_o as a function of energy, emission direction, and time at every point of the radiating surface of the source. Determining I_o(R_o,ω_o,t) in general requires calculating full-scale neutrino transport. With this quantity as the solution of the transport problem at hand, Eq. (3) can be evaluated directly (in general with energy dependence) by performing the integration over any surface that encloses the radiation (neutrino) emitting source and that lies outside the volume where radiation interacts with matter (i.e., the intensity I_o in all points on the chosen surface must be given in the reference frame of the observer and must be preserved on the way from the emission point to the observer).

Our ray-by-ray transport approximation, however, yields only the local neutrino energy density E(R_o,t) and the neutrino flux density F(R_o,t). To estimate the neutrinos radiated from every point of the neutrinosphere to the observer, we therefore have to develop an approach that yields a reasonable representation of the direction-dependent intensity as function of the quantities delivered by our transport approximation¹. To this end, we assume that the neutrino distribution is axisymmetric around the normal vector n_A at all points R_o. This implies that the direction dependence of the intensity I_o is described by the direction angle γ only (see Fig. 1), i.e., I_o = I_o(R_o,γ,t), and that the flux direction is given by n_A. Assuming further that I_o(γ) can be approximated by the lowest two terms of an expansion in spherical harmonics (as in the diffusion approximation), one can write (4)Because the radiation flux density is defined as the first angular moment of the intensity, one can easily verify that the numerical coefficient 3/2 of the dipole term ensures that the flux density F_o(R_o,t) (normal to the emitting surface element dA) is given by , if the radiating surface does not receive any incident neutrinos from outside (i.e., I_o(R_o,γ,t) = 0 for cosγ < 0)².

Inserting Eqs. (4) into (3), we find for the observable neutrino luminosity the expression (5)We further define an observable mean neutrino energy according to (6)where (7)is the observable neutrino number flux with ϵ being the mean energy of the neutrino energy spectrum radiated from point R_o.

Our 3D radiation hydrodynamics code computes the time-dependent neutrino energy flux density, F_o(R_o,t), and neutrino number flux density, F_n,o(R_o,t), through a sphere of radius R_o = |R_o| in dependence of the angular position Ω ≡ (θ,φ), but actually stores the related quantities (8)and (9)because these quantities are constant in the ray-by-ray approximation of the free streaming region. In this approximation, both the neutrino energy flux and the neutrino number flux are purely radial.

Using Eqs. (8) and (9), and the fact that with dΩ = sinθdθdφ for the special case of an emitting sphere of radius R_o, we can rewrite the general expression for the observable neutrino luminosity given in Eq. (5) in the form (10)and that of the observable neutrino number flux given in Eq. (7) in the form (11)where in both cases the integration is performed over the visible hemisphere³.

For the evaluation of the gravitational wave amplitude in Sect. 4.1.2 we will also need the quantity (12)and the corresponding angle-integrated quantity (13)For the later discussion of the results we finally define the surface-averaged neutrino flux density (14)where (15)is the total energy loss rate at time t from the supernova core to all directions, which (because of the flux variations over the sphere) is no directly observable quantity.

We have also analyzed the evolution of the neutrino flux asymmetry by calculating the angular pseudo-power spectrum of the neutrino energy flux variation (16)where Λ(Ω,t) and Λ(t) are defined in Eqs. (8) and (13), respectively. The pseudo-power spectrum is given by the decomposition of δ_Λ(Ω,t) in spherical harmonic coefficients (17)where is the respective (complex conjugate) spherical harmonics. For our mode analysis we actually used the pseudo-power coefficients C₀ ≡ |Λ₀₀|² and (18)for l > 0, respectively.

3.2. Results

The evolution of our models can be divided into four distinct phases (Figs. 2, 3).

• (1)

  The first phase, the quasi-spherical shock-expansion phase(Fig. 2, top row), lasts from shock formationshortly after core bounce to 80−150 ms, when convection sets in. During this phase the shock rapidly propagates out to a radius of  ~200 km, where its expansion comes to a halt.

• (2)

  The second phase, the hydrodynamically vigorous pre-explosion phase, comprises the growth of post-shock convection and of the standing accretion shock instability, SASI (Fig. 2, second row from top).

• (3)

  The post-explosion accretion phase begins when energy deposition by ν-heating in the post-shock layers becomes sufficiently strong to launch the explosion, and the total energy in the post-shock region ultimately becomes positive (see Sect. 2.3 for a definition). During this phase the shock accelerates outward while gas is still accreted onto the PNS. This process is commonly called “shock revival” (Fig. 2, third row from top). Non-radial instabilities during the latter two stages cause considerable temporal and angular fluctuations of the neutrino energy flux density as illustrated in Figs. 3−5. Besides the evolution of the shock radius, the figures show the surface-averaged neutrino light curve Λ(t), i.e., the total energy loss due to neutrinos versus time (Eqs. (13), (15)), together with the time evolution of the maximum and minimum values of Λ(Ω,t) (Eq. (8); the numerical evaluation is performed on an arbitrarily chosen sphere of 500 km radius). Distinct and high-amplitude spikes in Λ_max(Ω,t) are visible for several 100 ms and reflect violent post-shock convection, possible SASI activity, and anisotropic accretion fluctuations after the onset of the explosion. We have marked the explosion time t_exp (see Sect. 2.3, and Table 1) by a vertical dashed line in Figs. 3 and 4. The post-explosion accretion phase lasts until  ~500 ms (models W15-4 and N20-2) or  ~700 ms (model L15-3) depending on the progenitor.

• (4)

  During the post-accretion phase, the fourth and final phase characterizing the evolution of our models (Fig. 2, bottom row), gas infall to the proto-neutron star has come to an end and the newly formed neutron star looses mass at a low rate in a nearly spherical neutrino-driven wind. We find considerably less temporal variability and a lower level of angular variation (≲10%) of the neutrino emission during this fourth phase (Figs. 3−5). While in model L15-3 the amplitudes of the neutrino emission fluctuations decrease continuously, the other two models exhibit growing temporal emission variations (though at a lower level than the earlier variability) during a later stage (notice the decrease/increase in Λ_max − Λ_min in Figs. 3 and 4), which might be considered as a fifth evolutionary phase. This phase is associated with growing convective activity below the neutrinosphere. This PNS convection develops more or less strongly in the different models depending on the location of the convectively unstable region relative to the inner radial boundary of our computational domain.

Fig. 2 Snapshots of models W15-4 (left) and L15-3 (right) illustrating the four phases characterizing the evolution of our 3D models (see text for details). Each snapshot shows two surfaces of constant entropy marking the position of the shock wave (gray) and depicting the growth of non-radial structures (greenish). The time and linear scale are indicated for each snapshot.

We have evaluated the time evolution of the neutrino energy flux asymmetry by producing 4π-maps that show the relative angular variation ΔF_o/⟨F_o⟩ of the total (i.e., sum of all neutrino flavors) neutrino energy flux density across a sphere (normalized to the surface-averaged flux density; Eq. (14)). Several snapshots of this evolution are shown for model W15-4 in Fig. 5. The evolution of the typical angular scales of the fluctuations is reflected by the pseudo-power spectrogram of the electron neutrino energy flux variation (Eq. (16)) in Fig. 6, top panels, which give the color-coded pseudo-power coefficient distribution normalized to the maximum value versus time. The variation of the pseudo-power coefficients with angular mode number is shown in Fig. 7 at selected times of 200 ms (blue), 400 ms (red), and 1000 ms (black).

During the quasi-spherical shock expansion phase the level of angular fluctuations of F is low (≲10^-2), while the fluctuation amplitudes of the total neutrino energy flux density reach a level of several 10% during the hydrodynamically vigorous second phase and the post-explosion accretion phase, where a few distinct regions or even single spots with an angular size of 10° to 20° dominate the emission (Fig. 5, panels 2 and 3). The mode number l of the dominant angular perturbation scale is of no relevance during the first phase, as the maximum pseudo-power coefficient  (see Eq. (18)) is tiny  ≲10^-6 (Fig. 6, middle panels), i.e., the dominating l = 2 and l = 4 modes visible in the upper panels of Fig. 6 only reflect tiny angular perturbations imprinted presumably by the computational grid. When neutrino heating eventually causes significant non-radial flow during the second and third phases,  rises sharply to a level of  ~10^-3 (Fig. 6, middle panels), and the relative angular variations of the electron neutrino flux density grow to the several ten percent level (Fig. 6, bottom panels). The latter quantity gives the maximum minus the minimum flux density on the sphere divided by the angle-averaged flux density in percent. Compared to the total neutrino emission in Figs. 3−5, the temporal and angular variations in different directions are even more pronounced when considering the energy flux of the electron neutrinos or electron anti-neutrinos alone (Fig. 6), where angular variations can exceed 100% in all models during the pre-explosion and accretion phases, and peak values are close to 200% during short episodes (Fig. 6, lower panels).

Fig. 3 Shock radius (top) and total (i.e., summed over all flavors) energy loss rate due to neutrinos (bottom) as functions of time for model W15-4. In the upper panel, the black curve shows the angle-averaged mean shock radius, the blue (red) curve gives the maximum (minimum) shock radius, and the vertical dashed line marks the time of the onset of the explosion as defined in Sect. 2.3. In the lower panel, the blue and red curves show the time evolution of Λmax(Ω,t) and Λmin(Ω,t), the maximum and minimum value of Λ(Ω,t) (Eq. (8)) on a sphere of 500 km radius, respectively. The black line gives the corresponding surface-averaged value Λ(t) (Eq. (13)). Note that the luminosities imposed at the inner radial grid boundary are kept constant during the first second and later are assumed to decay like t−2/3.

Fig. 4 Same as Fig. 3 but for models L15-3 (uppermost two panels) and N20-2 (lowermost two panels), respectively.

Fig. 5 Neutrino flux asymmetry at 170 ms, 200 ms, 342 ms, 600 ms, and 1.3 s (from top to bottom), respectively. The 4π-maps show the relative angular variation ΔFo/⟨Fo⟩ of the total (i.e., sum of all neutrino flavors) neutrino energy flux density over a sphere (normalized to its angular average) for model W15-4. The maximum value is given in the lower right corner of each panel. Regions of higher emission are shown in bright yellow, while orange, red, green, and blue colors indicate successively less emission. Note that the color scale of each panel is adjusted to the maximum and minimum values at the corresponding time. The total energy loss rate due to neutrinos is given in the lower left corner.

During the vigorous pre-explosion phase including the post-explosion accretion stage, electron neutrinos and antineutrinos dominate the angular flux variations, while muon and tau neutrinos (accounting for roughly 50% of the total luminosity) exhibit essentially isotropic emission in all directions. The reason of this finding is that ν_e and  are produced almost exclusively by efficient charged-current reactions in the accretion region perturbed by non-radial fluid flows. The spectrogram of the two phases is characterized by initially very small-scale angular variations with l ≳ 12, which are associated with the onset of the Rayleigh-Taylor overturn activity, and which merge to continuously larger angular structures that correspond to l ≈ 1...4 modes toward the end of the accretion period at 0.4−0.6 s (depending on the model). This evolution is accompanied by a steady decrease of  to a level of  ~10^-5 and a reduction of the electron neutrino flux density variations from values well beyond 100% to a level of  ~10%, only (see Fig. 6, left panels).

When neutrino-energy deposition in the post-shock layers becomes sufficiently strong and the explosion is eventually launched at about 250 to 500 ms (depending on the model; Table 1), subsequent radial shock expansion rapidly diminishes the activity of the SASI and freezes post-shock convection. Single, longer lasting downdrafts of accretion flows are associated with isolated hot spots, where the variations of the total flux density can reach peak amplitudes of up to  ~70% (Fig. 5, panel 3). When accretion has ended, the amplitude of the angular variations of the total neutrino energy flux reduces to a level of a few percent (Figs. 3, 4), and the angular pattern of the emission becomes more uniform over the sphere, consisting of many spots with an angular size of  ~30° (Fig. 5, panel 4).

In the early post-accretion phase of model W15-4, 0.6   s ≲ t ≲ 0.8 s, the spectrogram indicates the presence of low-amplitude (), small-angular size (l ≳ 10) perturbations in the electron neutrino energy flux caused by some low-amplitude turbulent flow in and below the neutrinospheric region. When strong convection inside the PNS is encountered for t ≳ 0.8 s the spectrogram drastically changes, being dominated by angular modes with l = 4, but still with . The electron neutrino flux density variations rise somewhat to a level of 10% to 20%, and become manifest in the total energy loss rate, too (Figs. 3, 4).

Model N20-2 exhibits quite a similar behavior as model W15-4 except for the appearance of even larger (l ~ 3) angular structures clearly recognizable in the pseudo-power spectrogram between 1.0 s and 1.2 s (Figs. 6, 7). This differs from the behavior of model L15-3, where the amplitudes and angular size of the energy flux density variations remain small and even decrease in the post-explosion phase (Figs. 4, 6).

Fig. 6 Pseudo-power spectrogram of the electron-neutrino energy flux density (top row) for models W15-4 (left), L15-3 (middle), and N20-2 (right), respectively. The panels in the middle row show the corresponding maximum pseudo-power coefficient  as a function of time, and the panels in the lower row give the relative angular variation of the electron neutrino flux density (maximum minus minimum flux density on the sphere divided by the angle-averaged flux density in percent) with time.

The reason for the fluctuation behavior of the neutrino emission during the vigorous pre-explosion and post-explosion accretion phases has been discussed, but what causes the spatial and temporal variations during the post-accretion phase? Because the explosion is well on its way at this time, neither post-shock convection nor the SASI nor accretion can be responsible. Hence, there only remains non-radial gas flow in the outer layers of the proto-neutron star. Ledoux convection in the proto-neutron star thus may become visible eventually, i.e., its presence in the inner parts of the computational domain may become dominant in observable signals. This happens in models W15-4 and N20-2, where the level of the non-radial specific kinetic energy , volume-averaged over the computational domain, inside the neutrinosphere shows a steep rise at  ~0.8 s and  ~0.9 s, respectively (Fig. 8). These non-radial flows that develop in models W15-4 and N20-2 at late times also become manifest in all discussed quantities: Λ_max(Ω,t), Λ_min(Ω,t), , the dominant low l-modes (2 ≲ l ≲ 4), and relative angular flux-density variations. In contrast, no such effect is present in model L15-3 (see Fig. 8), where we find a steady decrease of Λ_max(Ω,t) − Λ_min(Ω,t), higher l-modes (l ≳ 10), smaller , and lower flux-density variation amplitudes than in models W15-4 and N20-2 (see Figs. 3, 4, and 6).

Simulations with fully self-consistent treatment of the PNS interior show the presence of convection inside the PNS, i.e., below the neutrinosphere (see Keil et al 1996; Buras et al. 2006; Dessart et al. 2006) more or less from the early post-bounce phase on. With the use of our inner radial grid boundary excising the inner parts of the PNS, and imposing neutrino luminosities at this boundary, convective activity is triggered only when the neutrino energy (or lepton number) inflow into the layers close to the grid boundary is faster than neutrino transport can carry away this energy (or lepton number). Then convectively unstable gradients develop and convective flows begin to carry the energy and lepton-number outward. Whether this happens or not depends on the boundary luminosities as well as on the location of the grid boundary within the density and temperature profiles of the PNS layers below the neutrinosphere. That location determines the efficiency of the neutrino transport and varies with the stellar progenitor, whose mass-infall rate decides how much mass accumulates in the near-surface layers of the PNS outside the inner grid boundary. The relative strength of the artificially imposed inflow of neutrino energy and lepton number compared to the efficiency of the neutrino transport on the grid, both sensitive to the location and contraction of the grid boundary on the one hand and the chosen values of the boundary luminosities on the other, therefore decides about when, where, and how strongly convective activity develops below the neutrinosphere.

Because the position of and the conditions imposed at the inner boundary can thus influence the neutrino emission properties, in particular during the post-accretion phase, our respective model predictions must be considered with care. While they do not allow us to make any definite statements concerning the neutrino signal of a particular progenitor model because of the neglected treatment of the inner parts of the proto-neutron star, the models nevertheless show that convective flows below the neutrinosphere are likely to imprint themselves on the neutrino emission, and hence also on the GW signal of core-collapse supernovae. A measurement of these signals may actually provide some insight into the conditions inside proto-neutron stars.

Because the neutrino energy flux density varies in our models both with latitude and longitude, the observable neutrino luminosity L_o(t) is obtained by an integral over the hemisphere visible to an observer (Eq. (10)). In Fig. 9 we show the observable electron neutrino and electron anti-neutrino luminosities for one chosen viewing direction for the three models W15-4, L15-3, and N20-2, respectively. The results for other directions look very similar with all characteristic features being independent of the observer position. We provide these quantities in addition to the total neutrino energy loss rate (Eqs. (13) and (15); Figs. 3 and 4), because their temporal evolutions are the ones expected to be measurable in the IceCube and Super-Kamiokande detectors. These detectors (mainly for ) will be sensitive to a combination of the observable neutrino luminosity L_o and the observable mean neutrino energy  ⟨ E ⟩ _o. Thus, we also provide in Fig. 9 the time evolution of the observable mean neutrino energy and of the combination , which (roughly) enters the IceCube detection rate of Cherenkov photons originating from the dominant inverse beta decay reaction (Lund et al. 2010)⁴. Again one can recognize the different evolution stages, and in particular the post-shock convection and SASI phase, during which the quantity  exhibits rapid low-amplitude variations for all three models. The level of the variations is a few percent (Fig. 9), which is considerably lower than that of the angular fluctuation amplitudes of the flux density, which reaches almost 100% for the total neutrino flux density (Fig. 5) and almost 200% for the electron neutrino and electron antineutrino flux densities (Fig. 6, lower panels). However, because the flux density variations are caused by a few individual hot spots covering only angular areas of size  ~(π/9)², the observable fluctuations (of L_o and  ⟨ E ⟩ _o) are lower by a factor of roughly (π/9)²/(2π) ~ 1/50. Some of this activity is also present at late times in the two models W15-4 and N20-2, where Ledoux convection develops in the simulated outer parts of the proto-neutron star (see discussion above).

Fig. 7 Pseudo-power coefficients of the electron-neutrino flux density as functions of angular mode number l at 200 ms (blue), 400 ms (red), and 1000 ms (black) for models W15-4 (top), L15-3 (middle), and N20-3 (bottom), respectively.

From the results presented above we conclude that the signals carry clear information about the postshock hydrodynamic activity, and about the duration and decay of the accretion period. Composition-shell interfaces present in the progenitor star can also leave an imprint. In model W15-4 the transition from the Fe-core to the Si-shell manifests itself in fast drops of the luminosities of ν_e and at  ~150 ms, when the mass accretion rate decreases steeply at the time the interface between the Fe-core and the Si-shell of the 15 M_⊙ progenitor falls through the shock.

4. Gravitational wave signature

Non-radial mass motions caused by gravity waves in the near-surface layers of the PNS, which are caused by the SASI and convection in the post-shock region as well as by convective activity inside the proto-neutron star (Murphy et al. 2009; Marek et al. 2009) (driven by Ledoux unstable lepton or entropy gradients) result in a time-dependent, aspherical density stratification that produces gravitational radiation. The anisotropic emission of neutrinos associated with the non-radial mass flow (see Sect. 3) contributes to the gravitational wave signal, too. We computed and analyzed the signature of this gravitational radiation for the 3D models discussed in Sect. 2.3.

4.1. Formalism

4.1.1. Non-radial mass flow

If a source is of genuine three-dimensional nature, as it is the case for our models, it is common to express the gravitational quadrupole radiation tensor, h^TT, in the transverse traceless gauge in the following tensorial form (19)(see, e.g., Misner et al. 1993). R denotes the distance between the observer and the source, and the unit linear-polarization tensors are given by with e_θ and e_φ being the unit polarization vectors in θ and φ-direction of a spherical coordinate system, and  ⊗  denoting the tensor product.

Fig. 8 Evolution of the non-radial specific kinetic energy  volume averaged over the computational domain inside the neutrinosphere for models W15-4 (solid), L15-3 (dashed), and N20-2 (dashed-dotted), respectively.

Fig. 9 Observable luminosity Lo (top row), observable mean energy  ⟨ E ⟩ o (middle row), and normalized quantity (bottom row) of electron neutrinos (left column) and electron anti-neutrinos (right column) as a function of time for three of our models. Although we only present the results for one particular observer direction here, the global behavior and characteristics are very similar for all viewing directions.

The wave amplitudes A₊ and A_× represent the only two independent modes of polarization in the TT gauge (Misner et al. 1993). In the slow-motion limit, they are obtained from linear combinations of the second time derivatives (evaluated at retarded time, and denoted by a double dot accent) of the components of the transverse traceless mass quadrupole tensor (Misner et al. 1993) We computed the latter using a post-Newtonian approach whereby the numerically troublesome second-order time derivatives of the mass quadrupole tensor components are transformed into much better tractable spatial derivatives. Following Nakamura & Oohara (1989) and Blanchet et al. (1990), the second-order time derivatives read in a Cartesian orthonormal basis (the spatial indices i and j run from 1 to 3) (24)where G is Newton’s gravitational constant, c the speed of light in vacuum, Φ_eff the effective Newtonian gravitational potential including the general relativistic “case A” correction of the monopole term according to Marek et al. (2006), ρ the mass-density, v_i the Cartesian velocity components, and ∂_i the partial derivative with respect to the coordinate xⁱ of a Cartesian basis.

We note that the integrand in Eq. (24) has compact support and is known to the (2nd order) accuracy level of the numerical scheme employed in the hydrodynamics code. It can easily be shown that evaluating the integral of Eq. (24) by an integration scheme (of at least 2nd order) is by one order of accuracy superior to twice applying numerical time-differentiation methods to quadrupole data given at discrete points of time (Finn & Evans 1990; Mönchmeyer et al. 1991).

Exploiting the coordinate transformation between the orthonormal Cartesian basis xⁱ and the orthonormal basis in spherical coordinates  (with ), the wave amplitudes A₊ and A_× (Eqs. (22) and (23)) are obtained from the following second time derivatives of the spherical components of the mass quadrupole tensor (Oohara et al. 1997; Scheidegger et al. 2008) where we used the abbreviation (28)Choosing φ = 0 one obtains the polarization modes (see, e.g., Misner et al. 1993) for θ = 0, and for θ = π/2, respectively. These expressions were already discussed in earlier investigations concerned with the evaluation of the gravitational wave signature of 3D core-collapse supernova models (Müller & Janka 1997; Fryer et al. 2004; Scheidegger et al. 2008, 2010).

The total energy radiated in the form of gravitational waves due to nonspherical mass flow is given in the quadrupole approximation by (see, e.g., Misner et al. 1993) (33)with , and the corresponding GW spectral energy density is given by (where ν denotes the frequency) (34)where (35)is the Fourier transform of .

4.1.2. Anisotropic neutrino emission

To determine the gravitational wave signal associated with the anisotropic emission of neutrinos, we follow Müller & Janka (1997) and use Eq. (16) of Epstein (1978) in the limit of a distant source, R → ∞, together with the approximation that the gravitational wave signal measured by an observer at time t is caused only by radiation emitted at time t′ = t − R/c. Hence, we take t − t′ = const. = R/c, i.e., we assume that only the neutrino pulse itself causes a gravitational wave signal, whereas memory effects, which prevail after the pulse has passed the observer, are disregarded.

With these simplifications, the dimensionless gravitational wave amplitudes of the two polarization modes are given in the transverse-traceless gauge for an observer located at a distance R along the z-axis of the observer frame by Müller & Janka (1997)(36)and (37)respectively. Here dΛ(Ω′,t′)/dΩ′ is given by Eq. (12) and denotes the total neutrino energy radiated at time t′ per unit of time into a solid angle dΩ′ in direction (θ′,φ′). Except for position-dependent factors the gravitational wave amplitudes are simply a function of this quantity provided by the ray-by-ray transport approximation (note that in Müller & Janka 1997, we used the symbol L_ν instead of Λ).

The angular integration, dΩ′ = −d(cosθ′)   dφ′, in Eqs. (36) and (37) extends over all angles θ′ and φ′ in the coordinate frame of the source (x′,y′,z′) that we identify with the (arbitrarily chosen) spherical polar coordinate frame to which the hydrodynamic results were mapped from the Yin-Yang grid employed in the simulations. For the evaluation of the polarisation modes we used the (asymptotic) values of dΛ(Ω,t)/dΩ extracted at a radius of 500 km from our 3D models.

Fig. 10 Relation between the source coordinate system (x′, y′, z′) and the observer coordinate system (x,y,z). Changing from the observer system to the source system involves a rotation by an angle α about the z′-axis to an intermediate coordinate system (x∗,y∗,z∗), followed by a rotation by an angle β about the y∗-axis (which thus is also the y-axis).

The angles θ and φ in Eqs. (36) and (37) are measured in the observer frame (x,y,z), while the neutrino luminosity is measured in the source frame (x′,y′,z′). To allow for an arbitrary orientation of the observer relative to the source, we introduce two viewing angles α ∈  [−π, + π]  and β ∈  [0,π]  (see Fig. 10). The coordinates measured in the observer frame are then related to the coordinates in the source frame by the following coordinate transformations and With these coordinate transformations and the relations between Cartesian coordinates (x,y,z) and spherical polar coordinates (r,θ,φ), we obtain These expressions relate the angular coordinates in the observer frame (θ,φ) to those in the source frame (θ′,φ′). For the special case α = 0 they were already presented by Kotake et al. (2009a). Using Eq. (49) and the equalities derived from Eqs. (44) to (46), the two polarization modes (Eqs. (36) and (37)) are given by (52)where S ∈ (+, ×) and Λ(t) is the angular integral of the neutrino energy radiated at time t per unit of time given in Eqs. (13) and (15). (53)are anisotropy parameters, which provide a quantitative measure of the time-dependent anisotropy of the emission in both polarization modes. Note that the evaluation of the anisotropy parameter α(t) defined in Eq. (29) of Müller & Janka (1997), which should not be confused with the observer angle α introduced in Fig. 10, does neither involve a dependence on observer angles (α,β) nor on the polarization mode.

The angular weight functions appearing in the above expression for the anisotropy parameters are given by (54)where Choosing α = 0 and β = π/2 the observer is located in the equatorial plane of the source (i.e., perpendicular to the source’s z′-axis) at the azimuthal position φ′ = 0. In that case one obtains simpler expressions for the angular functions (see also Kotake et al. 2009a) Note that for axisymmetric sources h_× = 0.

In general, the total energy E_GW(t) radiated to infinity by a source in form of gravitational waves until time t is given by (see, e.g., Misner et al. 1993; Greek indices run from 0 to 3, and repeated indices are summed over) (60)where the angular integration is performed over a two-sphere at spatial infinity , and n^μ = (0,1,0,0) is a unit spacelike vector in polar coordinates  {ct,r,θ,φ}  normal to . Denoting by  ⟨...⟩  an average over several wavelengths, the gravitational-wave energy-momentum tensor τ_μν is given in transverse-traceless gauge by (61)Thus, Eq. (60) can be rewritten as (62)where we have used the facts that , , and for radially outgoing gravitational radiation. Evaluating the double sum in Eq. (62) and using the relations and (see, e.g., Misner et al. 1993), we finally find (63)Inserting the expressions for h₊ and h_× given in Eqs. (52) into (63), we obtain for the energy E_N(t) radiated in form of gravitational waves until time t due to anisotropic neutrino emission (64)with dΩ_αβ = sinβ   dβ   dα and (65)The corresponding spectral energy density is given by (66)where is the Fourier transform of (67)For completeness we also provide an expression for the total energy radiated in form of gravitational waves until time t, i.e., due to anisotropic mass flow and neutrino emission. It is obtained by inserting the total GW amplitude, i.e., the sum of the amplitudes given by Eqs. (19) and (52) into (63), which leads to (68)

4.2. Results

Fig. 11 The four panels show the gravitational wave amplitudes (top) and spectrograms of dEM/dν (bottom; normalized to the absolute maximum) arising from non-spherical mass flow of models W15-2 (top left), W15-4 (top right), L15-2 (bottom left), and L15-3 (bottom right), respectively. Blue curves give the amplitude A+ at the pole (solid) and the equator (dotted), while red curves show the other independent mode of polarization A× from the same directions.

Fig. 12 Gravitational wave amplitudes (blue) and (red) due to anisotropic mass flow and neutrino emission as a function of time for models W15-4 (left), L15-3 (middle), and N20-3 (right), respectively. The solid curves show the amplitudes for an observer located above the north pole (α = β = 0; see Fig. 10) of the source, while the other curves give the amplitudes at the equator (α = 0, β = π/2).

Fig. 13 Asymmetry parameter αS of the neutrino emission (Eq. (53)) as a function of time for models W15-4 (left), L15-3 (middle), and N20-2 (right), respectively. The panels in the upper row show α+ (blue) and α× (red) for a particular observer direction, while the panels in the lower row give for both parameters the maximum and minimum values in all directions.

Although an observer can only measure the total gravitational wave amplitude, i.e., that due to the combined effect of non-radial flow and anisotropic neutrino emission, we will first discuss the GW signal of non-radial mass flow only, because it reflects the various phases of the post-bounce evolution already introduced in the discussion of the neutrino signal above.

Until post-shock convection and the SASI are eventually mature at around 150 ms, the GW signal is very small (Fig. 11). Note that our models are not able to follow the GW emission that is caused by prompt post-shock convection because of the excised inner region of the PNS (Marek et al. 2009). Later on, sizable g-mode activity is instigated in the outer layers of the proto-neutron star by convective overturn and the SASI during the hydrodynamically vigorous pre-explosion phase, and by the impact of anisotropic accretion flows during the subsequent post-explosion accretion phase (Marek et al. 2009). This g-mode activity is the cause of GW signals (Marek et al. 2009; Murphy et al. 2009; Yakunin et al. 2010), whose maximum amplitudes are on the order of a few centimeters centered around zero.

The GW frequency distribution possesses a very broad maximum in the range of 100 Hz to 500 Hz, and the frequency corresponding to this maximum slowly increases with time (Fig. 11). Partially already during the post-explosion accretion phase, but at latest when the shock wave starts to rapidly propagate to large radii between  ~0.4 s and  ~0.7 s (see Figs. 3 and 4), the GW amplitudes start to grow by about a factor of ten until approximately asymptoting at  ~0.9 s in the case of the models based on the progenitor W15, and at  ~1 s in the case of models based on the progenitor L15, respectively (Fig. 11). This growth of the amplitude is associated with the anisotropic expansion of the shock wave, and a positive/negative wave amplitude indicates a prolate/oblate explosion, respectively (Murphy et al. 2009).

While the GW amplitudes grow, the GW energy distribution dE_M/dν becomes narrower and dimmer, and the frequency at maximum power continues to increase. The latter effect was also observed in the 2D models of Murphy et al. (2009). At late times, the GW signal of the W15 models clearly signifies the convective activity inside the proto-neutron star through low-amplitude, high-frequency fluctuations around the asymptotically roughly constant mean GW amplitudes, while no such fluctuations are present in the case of the L15 models (see discussion of the neutrino signal above). This model discrepancy is also evident from the energy spectrograms, which do exhibit a pronounced broad maximum (between  ~350 Hz and  ~550 Hz) at t > 0.8 s in the case of the W15 models, but none for the L15 ones. Furthermore note that until the end of the simulations the frequency of the maximum of dE_M/dν has increased from around 100 Hz to almost 500 Hz for the former models (owing to the increasing speeds of mass motions in the postshock region at times  ≲0.5 s and because of the increased compactness of the proto-neutron star at times  ≳0.6 s, respectively).

The behavior of the total (matter plus neutrinos) GW amplitudes is significantly different from that of the flow-only GW amplitudes for models that exhibit PNS convection below the neutrinosphere, i.e., for the models based on the progenitors W15 and N20. Particularly at late times, anisotropic neutrino emission causes a continuing growth of the GW amplitudes (instead of a saturation) in these models, while this is not the case for the L15 models (see Fig. 12, and compare with Fig. 11). The latter behavior is also reflected in the time evolution of the asymmetry parameter α_S (Eq. (53)) of the neutrino emission (Fig. 13). The asymmetry parameter is practically zero in model L15-3 at late times, while it remains, after having temporarily grown to values beyond about 0.4−0.5%, at the level of  ~0.3% until the end of the simulations in models W15-4 and N20-2.

The final GW amplitudes are up to a factor of two to three higher when taking the contribution of anisotropic neutrino emission into account. In contrast, the amount of energy radiated in the form of GW, which is proportional to the third time-derivative of the quadrupole moment and hence proportional to the time derivative of the GW amplitude is only insignificantly changed, and is practically constant for all simulated models after the onset of the explosion (see Fig. 14). The integral value of the GW energy radiated by neutrinos is low (≲1%) compared to that emitted by matter through the slow variation of the GW neutrino amplitude with time, i.e., its time derivative is much smaller than that of the GW matter amplitude. For this reason we also abstained from evaluating the total energy radiated in form of GW (Eq. (68)). It differs little from that caused by anisotropic matter flow alone (Eq. (33)), because the mixed term in Eq. (68), resulting from the square of the sum of the matter and neutrino parts, contributes  ≲10% to the total radiated GW energy, and the pure neutrino term  ≲1%. Figure 14 also shows that the (small) contribution of anisotropic neutrino emission to the radiated GW energy is enhanced at late times when proto-neutron star convection occurs below the neutrinosphere, as it is the case for models W15-4 and in particular N20-2.

Fig. 14 Energy emitted in form of gravitational waves due to anisotropic mass flow (top panel) and due to anisotropic neutrino emission (bottom panel) as a function of time for models W15-4 (solid), L15-3 (dashed), and N20-2 (dash-dotted), respectively.

Fig. 15 Gravitational wave amplitudes due to anisotropic mass flow and neutrino emission, (top) and (bottom), as functions of the observer angles (see Fig. 10) for model W15-4 at 1.3 s past bounce. The white contours give the locations, where the amplitudes are zero. Yellow and red areas indicate positive amplitudes, green and blue negative ones.

The variation of the total GW amplitudes with observer angle is illustrated in Fig. 15 for model W15-4 at 1.3 s (when the simulation was stopped). Both the amplitude variations and the typical angular size of the speckled GW emission are similar for all other simulated models. The model-independent level of the amplitude variations is also supported by Fig. 12 when comparing various amplitudes at any given (late) time.

The (normalized) amplitude spectrograms of the total gravitational wave amplitudes (d(A _+ , ×  + Rh _+ , × )/dν; Figs. 16 and 17) illustrate two model-independent findings. Firstly, during the hydrodynamically vigorous pre-explosion and post-explosion accretion phases (0.2 ≲ t ≲ 0.5−0.7 s) the spectra of all models are characterized by some power at low frequencies (≲100 Hz) and a broad power maximum at frequencies  ~200 Hz and another weak one at  ~800 Hz. The latter broad maximum at high frequency is more pronounced in the models based on the W15 and N20 progenitors and in the cross polarization GW mode. Secondly, during the post-accretion phase (t ≳ 0.7 s) the spectra of all models are dominated by a low-frequency (≲40 Hz) contribution peaked toward the lower end of the spectrogram. In the models where PNS convection occurs below the neutrinosphere (models W15 and N20) we also find a double-peaked high-frequency contribution decreasing/increasing from  ~700 Hz (400 Hz) at t ~ 0.8 s, and eventually merging into a single power maximum at  ~500 Hz at t ~ 1.2 s. Again, this contribution is more pronounced for the cross polarization GW mode.

The spectra of the total GW amplitudes are dominated by the contribution from non-isotropic neutrino emission at low (≤100 Hz) frequencies (Figs. 16 and 17). At higher frequencies (≥100 Hz) the spectra of model W15-4 show two pronounced maxima (at 100−200 Hz and 600−800 Hz, respectively) at all times. These maxima are also present in model L15-3 at times  ≤0.7 s, the high-frequency one being, however, much less pronounced. The lower maximum (at 100−200   Hz) results from g-mode activity in the PNS surface instigated by non-radial flow (SASI, accretion) in the post-shock region until  ~0.5−0.7 s. At later times PNS convection is responsible for the peak between 300 and 500 Hz. We have proposed this explanation already for the corresponding maxima present in the GW energy spectrograms arising from non-spherical mass flow (Fig. 11), and discussed why the frequencies of these maxima increase with time. The source of the high-frequency maximum (600−800 Hz) is unclear, but a further detailed analysis shows that (i) the maximum is solely caused by non-radial gas flow, i.e., it is not connected to neutrinos; (ii) it does not result from stellar layers below the neutrino sphere but from those close to or slightly above it; and (iii) does not depend on the position of the observer.

Note that the high-frequency maximum present in the amplitude spectrograms is strongly suppressed in the corresponding energy spectrograms (Fig. 11), because the latter involve the squared time derivatives of the amplitudes. Thus, the already high ratio of the low- and high-frequency maxima in the amplitude spectrograms (about two orders of magnitude) translates into an even higher ratio for the energy spectrogram maxima, rendering the high-frequency maximum practically invisible.

Fig. 16 Normalized (to the absolute maximum) amplitude spectrograms of the total gravitational wave amplitudes A+ + Rh+ (left panels) and A× + Rh× (right panels) at the pole for model W15-4. The lower panels show the spectrograms in the frequency range 5 Hz to 100 Hz, and the upper ones in the frequency range 100 Hz to 1 kHz.

Fig. 17 Same as Fig. 16, but for model L15-3.

5. Discussion and conclusions

Based on a set of three-dimensional (3D) parametrized neutrino-driven supernova explosion models of non-rotating 15 and 20 M_⊙ stars, employing a neutrino transport description with a gray spectral treatment and a ray-by-ray approximation of multi-dimensional effects (the scheme is applicable in the regime outside the dense neutron star core, i.e., around and outside the neutrinosphere), we evaluated both the time-dependent and direction-dependent neutrino and gravitational-wave emission of these models. To this end we presented the formalism necessary to compute both the observable neutrino and gravitational wave signals for a three-dimensional, spherical source. For the neutrino signal we presented formulas that allow one to estimate the apparent luminosity when the local flux density F on a sphere is known and the ansatz of Eq. (4) about the angular distribution of the intensity can be made. While in general the location of the neutrino-decoupling surface has to be suitably defined (e.g., by a criterion based on the optical depth), our ray-by-ray transport scheme implies that r²F  =  const. in the free-streaming limit in every direction. Thus, it spares us such a definition of the neutrinosphere for evaluating the angular integration for the observable luminosity. Concerning the gravitational-wave analysis, we extended and generalized previous studies, where the source was either assumed to be axisymmetric or where the formulas for the signals of a 3D source were only given for special observer directions.

Our models followed the evolution from shortly after core bounce up to more than one second into the early cooling evolution of the PNS without imposing any symmetry restrictions and covering a full sphere. The extension over such a relatively long evolution time in 3D was possible through the usage of an axis-free overset grid (the Yin-Yang grid) in spherical polar coordinates, which considerably eases the CFL time-step restriction and avoids axis artifacts. A central region, the dense inner core of the proto-neutron star, was excised from the computational domain and replaced by an inner, time-dependent radial boundary condition and a gravitating point mass at the coordinate origin. Explosions in the models were initiated by neutrino heating at a rate that depends on suitably chosen values of the neutrino luminosities imposed at the inner radial boundary.

The post-bounce evolution of our models can be divided into four distinct phases (Fig. 3). The first phase, the quasi-spherical shock-expansion phase, lasts from shock formation shortly after core bounce to 80−150 ms, when convection sets in. The second phase, the hydrodynamically vigorous pre-explosion phase, comprises the growth of post-shock convection and of the standing accretion shock instability (SASI). The post-explosion accretion phase begins when energy deposition by ν-heating in the post-shock layers becomes sufficiently strong so that the total energy in the post-shock region ultimately becomes positive. During this phase the shock accelerates outward while gas is still accreted onto the PNS. This process is commonly called “shock revival”. The duration of the latter two phases depends on the progenitor. During the post-accretion phase, the fourth and final phase characterizing the evolution of our models, accretion ends and the proto-neutron star develops a nearly spherical neutrino-driven wind.

The neutrino emission properties (fluxes and effective spectral temperatures) of our 3D models exhibit the generic time-dependent features already known from 2D (axisymmetric) models (e.g., Buras et al. 2006; Scheck et al. 2006; Marek et al. 2009; Brandt et al. 2011), showing fluctuations over the neutron star surface on different spatial and temporal scales. We found that non-radial mass motions caused by the SASI and convection in the neutrino-heated hot-bubble region as well as by PNS convection below the neutrinosphere give rise to a time-dependent, anisotropic emission of neutrinos, particularly of electron neutrinos and anti-neutrinos, and thus also to the emission of gravitational waves. Because very prominent, quasi-periodic sloshing motions of the shock due to the standing accretion-shock instability as visible in 2D simulations are absent and the emission from different surface areas facing an observer adds up incoherently, the modulation amplitudes of the measurable neutrino luminosities and mean energies are significantly lower than predicted by 2D models (for 2D results see Marek et al. 2009; Brandt et al. 2011).

During the quasi-spherical shock expansion phase shortly after bounce the level of temporal and angular fluctuations of the neutrino emission is low (≲10^-2). In contrast, the fluctuation amplitudes reach a level of several 10% of the average values during the hydrodynamically vigorous pre-explosion phase and the post-explosion accretion phase, where a few distinct, highly time-variable regions or even short-lived single spots with an angular size of 10° to 20° are responsible for the brightest emission maxima. As the outward shock expansion is well on its way in the post-explosion accretion phase, still existing accretion downdrafts can be responsible for similar fluctuations in the neutrino emission, though the number of corresponding hot spots decreases with diminishing accretion. When accretion has ended and the post-accretion phase has started, directional variations can be caused by the occurrence of Ledoux convection in the outer layers of the proto-neutron star, which we indeed observe in models based on two of our three progenitors (see also the discussion of the influence of the inner radial boundary condition below). The temporal and angular variations of the emission in different directions are even more pronounced when considering the energy flux of the electron neutrinos or electron anti-neutrinos alone (instead of the emission in all neutrino flavors). In that case the angular variations of local flux densities can exceed 100% in all models during the pre-explosion and post-explosion accretion phases, and the peak values can be close to 200% during short episodes. The total energy loss rates in neutrinos and the observable luminosities as surface-integrated quantities, however, are much smoother in time during all phases, showing fluctuation amplitudes of at most several percent.

The gravitational wave emission also exhibits the generic time-dependent features already known from 2D (axisymmetric) models, but the 3D wave amplitudes are considerably lower (by a factor of 2−3) than those predicted by 2D models (Müller et al. 2004; Marek et al. 2009; Murphy et al. 2009; Yakunin et al. 2010) owing to less coherent mass motions and neutrino emission. Note in this respect that the GW quadrupole amplitudes, which are usually quoted for 2D models (), have to be multiplied by a geometric factor (which is equal to  ≈0.27 for θ = 90°). Violent, non-radial hydrodynamic mass motions in the accretion layer and their interaction with the outer layers of the proto-neutron star give rise to a GW signal with an amplitude of  ~2−4 cm in the frequency range of  ~100 Hz to  ~400 Hz, while anisotropic neutrino emission is responsible for a superimposed low-frequency evolution of the wave amplitude, which thus can grow to maximum values of 10−20 cm. Variations of the mass-quadrupole moment caused by convective activity inside the nascent neutron star contribute a high-frequency component (300−600 Hz) to the GW signal during the post-accretion phase. The GW signals exhibit strong variability between the two polarizations, different explosion simulations and different observer directions, and besides common basic features do not possess any template character.

Finally we would like to reflect on some of the deficiencies of the presented 3D models. Because of the ray-by-ray treatment of the ν-transport, the directional variations of the neutrino emission in response to local inhomogeneities in the star may be overestimated (Ott et al. 2008; Brandt et al. 2011). However, when we evaluate observable signals, these artificial effects are mostly compensated for by integrations of the neutrino flux densities over the surface areas visible to observers from different viewing directions (see Eqs. (5), (7), (10), (11)), or by the integration of the neutrino energy loss in all directions (Eqs. (13), (53)).

Another deficiency concerns the usage of the inner radial grid boundary, because of which our simulations do not fully include (either in space or time) the convective flow occurring in the PNS interior after core bounce (Keil et al 1996; Buras et al. 2006; Dessart et al. 2006). Moreover, convective activity in the simulated outer layers of the PNS occurs for special conditions: it is triggered only when the artificially imposed inflow of neutrino energy and lepton number through the inner radial boundary into the adjacent layers is faster than the neutrinos can carry away this energy or lepton number. Whether this is the case sensitively depends on the employed neutrino-transport approximation, but also on the location and contraction of the grid boundary, the chosen values of the boundary luminosities, and on the stellar progenitor. Its mass-infall rate decides how much mass accumulates in the near-surface layers of the PNS outside the inner grid boundary. Because the position of and the conditions imposed at the inner boundary can thus influence the neutrino emission properties, in particular during the post-accretion phase, our respective model predictions must be considered with care. While they do not allow us to make any definite statements concerning the detailed neutrino signal of a particular progenitor model due to the neglected treatment of the inner parts of the proto-neutron star, the models nevertheless show that convective flows below the neutrinosphere are likely to imprint themselves on the neutrino emission, and hence also on the GW signal of core-collapse supernovae. A measurement of these signals may actually provide some insight into the conditions inside proto-neutron stars.

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft through the Transregional Collaborative Research Centers SFB/TR 27 “Neutrinos and Beyond” and SFB/TR 7 “Gravitational Wave Astronomy” and the Cluster of Excellence EXC 153 “Origin and Structure of the Universe” (http://www.universe-cluster.de).

References

All Tables

All Figures

	Fig. 1 Sketch illustrating the various quantities involved when deriving the observable luminosity of a radiating source, whose visible surface is shaded in blue. When the observer is at infinity, the emitting surface is a sphere, and Ro is measured from the origin of this sphere, one obtains α = γ.
In the text

	Fig. 2 Snapshots of models W15-4 (left) and L15-3 (right) illustrating the four phases characterizing the evolution of our 3D models (see text for details). Each snapshot shows two surfaces of constant entropy marking the position of the shock wave (gray) and depicting the growth of non-radial structures (greenish). The time and linear scale are indicated for each snapshot.
In the text

Fig. 3 Shock radius (top) and total (i.e., summed over all flavors) energy loss rate due to neutrinos (bottom) as functions of time for model W15-4. In the upper panel, the black curve shows the angle-averaged mean shock radius, the blue (red) curve gives the maximum (minimum) shock radius, and the vertical dashed line marks the time of the onset of the explosion as defined in Sect. 2.3. In the lower panel, the blue and red curves show the time evolution of Λmax(Ω,t) and Λmin(Ω,t), the maximum and minimum value of Λ(Ω,t) (Eq. (8)) on a sphere of 500 km radius, respectively. The black line gives the corresponding surface-averaged value Λ(t) (Eq. (13)). Note that the luminosities imposed at the inner radial grid boundary are kept constant during the first second and later are assumed to decay like t−2/3.

In the text

	Fig. 4 Same as Fig. 3 but for models L15-3 (uppermost two panels) and N20-2 (lowermost two panels), respectively.
In the text

Fig. 5 Neutrino flux asymmetry at 170 ms, 200 ms, 342 ms, 600 ms, and 1.3 s (from top to bottom), respectively. The 4π-maps show the relative angular variation ΔFo/⟨Fo⟩ of the total (i.e., sum of all neutrino flavors) neutrino energy flux density over a sphere (normalized to its angular average) for model W15-4. The maximum value is given in the lower right corner of each panel. Regions of higher emission are shown in bright yellow, while orange, red, green, and blue colors indicate successively less emission. Note that the color scale of each panel is adjusted to the maximum and minimum values at the corresponding time. The total energy loss rate due to neutrinos is given in the lower left corner.

In the text

Fig. 6 Pseudo-power spectrogram of the electron-neutrino energy flux density (top row) for models W15-4 (left), L15-3 (middle), and N20-2 (right), respectively. The panels in the middle row show the corresponding maximum pseudo-power coefficient  as a function of time, and the panels in the lower row give the relative angular variation of the electron neutrino flux density (maximum minus minimum flux density on the sphere divided by the angle-averaged flux density in percent) with time.

In the text

	Fig. 7 Pseudo-power coefficients of the electron-neutrino flux density as functions of angular mode number l at 200 ms (blue), 400 ms (red), and 1000 ms (black) for models W15-4 (top), L15-3 (middle), and N20-3 (bottom), respectively.
In the text

	Fig. 8 Evolution of the non-radial specific kinetic energy  volume averaged over the computational domain inside the neutrinosphere for models W15-4 (solid), L15-3 (dashed), and N20-2 (dashed-dotted), respectively.
In the text

Fig. 9 Observable luminosity Lo (top row), observable mean energy  ⟨ E ⟩ o (middle row), and normalized quantity (bottom row) of electron neutrinos (left column) and electron anti-neutrinos (right column) as a function of time for three of our models. Although we only present the results for one particular observer direction here, the global behavior and characteristics are very similar for all viewing directions.

In the text

	Fig. 10 Relation between the source coordinate system (x′, y′, z′) and the observer coordinate system (x,y,z). Changing from the observer system to the source system involves a rotation by an angle α about the z′-axis to an intermediate coordinate system (x∗,y∗,z∗), followed by a rotation by an angle β about the y∗-axis (which thus is also the y-axis).
In the text

Fig. 11 The four panels show the gravitational wave amplitudes (top) and spectrograms of dEM/dν (bottom; normalized to the absolute maximum) arising from non-spherical mass flow of models W15-2 (top left), W15-4 (top right), L15-2 (bottom left), and L15-3 (bottom right), respectively. Blue curves give the amplitude A+ at the pole (solid) and the equator (dotted), while red curves show the other independent mode of polarization A× from the same directions.

In the text

	Fig. 12 Gravitational wave amplitudes (blue) and (red) due to anisotropic mass flow and neutrino emission as a function of time for models W15-4 (left), L15-3 (middle), and N20-3 (right), respectively. The solid curves show the amplitudes for an observer located above the north pole (α = β = 0; see Fig. 10) of the source, while the other curves give the amplitudes at the equator (α = 0, β = π/2).
In the text

	Fig. 13 Asymmetry parameter αS of the neutrino emission (Eq. (53)) as a function of time for models W15-4 (left), L15-3 (middle), and N20-2 (right), respectively. The panels in the upper row show α+ (blue) and α× (red) for a particular observer direction, while the panels in the lower row give for both parameters the maximum and minimum values in all directions.
In the text

	Fig. 14 Energy emitted in form of gravitational waves due to anisotropic mass flow (top panel) and due to anisotropic neutrino emission (bottom panel) as a function of time for models W15-4 (solid), L15-3 (dashed), and N20-2 (dash-dotted), respectively.
In the text

	Fig. 15 Gravitational wave amplitudes due to anisotropic mass flow and neutrino emission, (top) and (bottom), as functions of the observer angles (see Fig. 10) for model W15-4 at 1.3 s past bounce. The white contours give the locations, where the amplitudes are zero. Yellow and red areas indicate positive amplitudes, green and blue negative ones.
In the text

	Fig. 16 Normalized (to the absolute maximum) amplitude spectrograms of the total gravitational wave amplitudes A+ + Rh+ (left panels) and A× + Rh× (right panels) at the pole for model W15-4. The lower panels show the spectrograms in the frequency range 5 Hz to 100 Hz, and the upper ones in the frequency range 100 Hz to 1 kHz.
In the text

	Fig. 17 Same as Fig. 16, but for model L15-3.
In the text