Effects of LESA in Three-Dimensional Supernova Simulations
with Multi-Dimensional and Ray-by-Ray-plus Neutrino Transport
A set of eight self-consistent, time-dependent supernova (SN) simulations in three spatial dimensions (3D) for 9 and 20 progenitors is evaluated for the presence of dipolar asymmetries of the electron lepton-number emission as discovered by Tamborra et al. and termed lepton-number emission self-sustained asymmetry (LESA). The simulations were performed with the Aenus-Alcar neutrino/hydrodynamics code, which treats the energy- and velocity-dependent transport of neutrinos of all flavors by a two-moment scheme with algebraic M1 closure. For each of the progenitors, results with fully multi-dimensional (FMD) neutrino transport and with ray-by-ray-plus (RbR+) approximation are considered for two different grid resolutions. While the 9 models develop explosions, the 20 progenitor does not explode with the employed version of simplified neutrino opacities. In all 3D models we observe the growth of substantial dipole amplitudes of the lepton-number (electron neutrino minus antineutrino) flux with stable or slowly time-evolving direction and overall properties fully consistent with the LESA phenomenon. Models with RbR+ transport develop LESA dipoles somewhat faster and with temporarily higher amplitudes, but the FMD calculations exhibit cleaner hemispheric asymmetries with a far more dominant dipole. In contrast, the RbR+ results display much wider multipole spectra of the neutrino-emission anisotropies with significant power also in the quadrupole and higher-order modes. Our results disprove speculations that LESA is a numerical artifact of RbR+ transport. We also discuss LESA as consequence of a dipolar convection flow inside of the nascent neutron star and establish, tentatively, a connection to Chandrasekhar’s linear theory of thermal instability in spherical shells.
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany \move@AU\move@AF\@affiliationPhysik Department, Technische Universität München, James-Franck-Straße 1, D-85748 Garching, Germany
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany \move@AU\move@AF\@affiliationPhysik Department, Technische Universität München, James-Franck-Straße 1, D-85748 Garching, Germany
Astrophysical Big Bang Laboratory, RIKEN Cluster for Pioneering Research, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan \move@AU\move@AF\@affiliationMax-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany
Tamborra et al. (2014a) and Janka et al. (2016) reported the discovery of a large hemispheric asymmetry of the total electron lepton-number emission from the newly formed neutron star in a set of three-dimensional (3D) core-collapse supernova (CCSN) simulations, which made use of the ray-by-ray-plus (RbR+) approximation to treat the multi-dimensionality of the neutrino transport, thereby ignoring the nonradial components of the neutrino flux vector. They termed this stunning, new phenomenon “lepton-number emission self-sustained asymmetry” (LESA). It is characterized by a substantial amplitude of the dipole component (in some cases and phases even exceeding the monopole) of the lepton-number flux, which is defined by the difference of the radial number flux of electron neutrinos and antineutrinos at spatial point ,
The radial number flux of any neutrino species is given by
with the radial energy-flux density and the neutrino energy . Both the dipole amplitude and the direction of this dipole were observed to exhibit basically stationary or slowly migrating behavior over time intervals of at least several hundred milliseconds covered by the simulations.
In contrast, Dolence et al. (2015) did not observe any evidence for LESA when carrying out 2D CCSN simulations with a multi-group, flux-limited diffusion scheme for neutrino transport. For this reason, Sumiyoshi et al. (2015) speculated that LESA may be an artifact of the RbR+ approximation. In a later study, Vartanyan et al. (2018) also did not find indications of LESA when conducting 2D simulations with a fully multi-dimensional two-moment scheme. However, Tamborra et al. (2014a) argued that the behavior of LESA was considerably different between their 2D and 3D models. While in the 2D case the LESA dipole direction, which can only be aligned with the artificial axis of symmetry, flips between both polar directions on a time scale of some 10 ms, strongly affected by SASI mass motions around the proto-neutron star, only their 3D models showed a persistent and directionally stable LESA dipole.
Roberts et al. (2016) also investigated their 3D simulations with a fully multi-dimensional two-moment neutrino-transport scheme for effects of LESA, but could not trace any evidence of this phenomenon. Although their models did not take into account velocity-dependent terms in the neutrino-moment equations and ignored inelastic neutrino scattering on electrons, this group interpreted the absence of LESA in contrast to the SN models of Tamborra et al. (2014a) not necessarily as a consequence of these differences in the neutrino treatment. Instead, Roberts et al. (2016) hypothesized that it could be linked, for example, to the fact that neutron-star convection showed up rather late in their simulations.
The first independent confirmation of the LESA phenomenon in 3D SN models with detailed neutrino transport different from those of the Garching group has recently been reported by O’Connor & Couch (2018b) and has possibly be seen in a 3D simulation by Vartanyan et al. (2018), too. Both groups also applied a fully multi-dimensional two-moment scheme, for which reason their results provide support for LESA not being an artifact of the RbR transport approximation. Although O’Connor & Couch (2018b) observed a substantial LESA dipole only in their single full-sphere calculation including velocity-dependent terms in the transport, they argued that the appearance of LESA is not directly linked to this improvement but appears to be a consequence of the stronger proto-neutron star convection seen in the model with velocity dependence. This conclusion may offer an explanation why LESA was absent in the models of Roberts et al. (2016), and it is consistent with the finding by Janka et al. (2016) that LESA is weaker in rotating and further diminished in rapidly rotating models, because angular momentum gradients in the neutron star suppress convective activity.
Here we present an evaluation of a set of eight new 3D simulations published in a companion paper by Glas et al. (2018, Paper I) for the LESA phenomenon. The simulations were performed with the Aenus-Alcar code for an exploding 9 progenitor and a nonexploding 20 star with two chosen grid resolutions. The transport treatment for all three neutrino flavors was based on an energy- and velocity-dependent two-moment scheme (combining energy and momentum equations with an analytic “M1” Eddington closure relation), which was applied in a fully multi-dimensional (FMD) version as well as an RbR+ mode. In all of these eight simulations (four for the low-mass and four for the high-mass progenitor), we find clear evidence for the LESA with all of its characteristic properties described by Tamborra et al. (2014a) and Tamborra et al. (2014b).
Our paper is structured as follows. In Section 2 we briefly summarize the numerical methods and computed models, refering the reader for more detailed information to Paper I. In Section 3 we present the results from our LESA analysis of the 3D simulations with FMD and RbR+ transport treatments. In Section 4 we provide a more detailed description of the proto-neutron star convection layer in the phase where LESA with a dominant dipole mode is fully developed. In Section 5 we construct a tentative link between the growth and initial evolution of different LESA multipoles and the onset of thermally driven convection discussed via a linear stability theory by Chandrasekhar (1961). In Section 6 we wrap up by an assessment of our results and conclusions.
2 Numerics, Simulation Setup,
and 3D Models
This study of LESA is based on the 3D simulations that were performed in Paper I. We will briefly summarize the numerical methods and the simulation setup in the following, and refer the reader to Paper I for detailed descriptions.
We employ the Aenus-Alcar code (Obergaulinger, 2008; Just et al., 2015, 2018), which solves the hydrodynamics equations and the velocity- and energy-dependent neutrino-moment equations in spherical polar coordinates (given by radius , polar angle , and azimuth angle ) with a Godunov-type, directionally unsplit finite-volume scheme and high-order reconstruction methods with an approximate Riemann solver to obtain cell-interface fluxes. We make use of the microphysical equation of state (SFHo) from Steiner et al. (2013) to close the system of hydrodynamics equations and include spherically symmetric gravitational self-interaction with general relativistic corrections (see Marek et al., 2006, case A). In the “M1” two-moment neutrino-transport scheme the energy density and energy-flux density (with ) of neutrinos are evolved for electron neutrinos , anti-electron neutrinos , and a third species that represents all four heavy-lepton neutrinos. The system of moment equations is closed by an algebraic relation for non-evolved moments like the radiation pressure tensor that depends on and . In contrast to the FMD two-moment transport scheme, the non-radial neutrino-flux components are set to zero in calculations with the RbR+ approximation. In the 3D simulations discussed here, we used 15 bins for the neutrino energy , that are logarithmically spaced in the interval for all neutrino species. The employed neutrino interactions are specified in Paper I.
In total, we conducted eight self-consistent, time-dependent CCSN simulations in 3D for two different non-rotating, solar-metallicity progenitor models (see Table 2): the progenitor model from Woosley & Heger (2007, labeled s20 in the following), and the solar-metallicity progenitor model from Woosley & Heger (2015, labeled s9.0 in the following). For simulations of both progenitor models, we changed the numerical grid resolution by varying the number of grid cells (, , and , respectively). Simulations with “low” (“L”) resolution use , , and zones, whereas models with “high” (“H”) resolution use twice the number of cells, i.e., , , and . In all cases, the radial grid ranges from 0 to with logarithmically increasing cell sizes, and the angular grids cover the entire solid angle of the full sphere with cell sizes of (low) and (high) except for two bigger lateral zones at each of both poles of the spherical coordinate system (for details, see Paper I). In contrast, the azimuthal grid is equidistant. For each progenitor model and grid resolution, we performed simulations with the RbR+ approximation and the FMD neutrino-transport scheme.
The small cell sizes near the center of the polar coordinate system lead to severe constraints on the numerical time step, which we circumvent by evolving the innermost in spherical symmetry. Thus, we neglect angular variations in the (spherically symmetric) innermost region of the proto-neutron star. This is an acceptable approach as long as the central core of remains convectively stable. Practically this is only fulfilled until after bounce but can become problematic in some of our simulations at later times (see end of Section 4 and Section 5).
3 LESA in FMD and RbR+ Models
In order to shed more light on the rather unclear situation pictured by the diverse (but not necessarily conflicting) results of the presence or absence of LESA in models reported in the literature, we also evaluate our 3D simulations for signatures of LESA. We can exploit the unique possibility to directly asses the influence of the neutrino-transport schemes, FMD vs. RbR+. To this end we first perform a multipole analysis and decompose the radial lepton-number flux in the laboratory frame (i.e., the rest frame of the stellar center),
into real spherical harmonics, , of degree and order . The radial number fluxes of neutrinos in the lab frame are obtained by the transformation
where is the radial fluid velocity and is the number density of neutrinos in the co-moving frame. The multipole coefficients at a given radius are calculated by integrating the flux over the whole spherical surface:
The moments of the lepton-number flux are then obtained as
Here, we use a normalization for which the monopole moment, , is identified with the total lepton-number flux through the surface of a sphere with radius 333We point out that the normalization of the spherical harmonics coefficients in the present case differs from what we used in the decomposition of the shock radius in Equations (11) and (12) of Paper I. There, all coefficients are rescaled by a factor of compared to . With this normalization the monopole moment stands for the mean shock radius..
Figure 3 displays the time-dependent evolution of the amplitudes of the monopole (, dotted lines) and of the dipole (, solid lines, scaled by a factor of 3 compared to Equation (6)) of the total electron lepton-number flux (evaluated as lab-frame quantities at a radius of ) for all of our 3D simulations. To visualize the evolution of the dipole direction, the polar and azimuthal angles of the dipole vector444For the dipole direction of the lepton-number flux, we use the multipole coefficients as components of the direction vector, which we then transform into spherical polar coordinates. This procedure is equivalent to the analysis of the dipole direction of the shock radius (see Equations (11) and (12) and the associated explanations in Paper I). are shown as functions of time after bounce in Figure 3.
All of our simulations, both with RbR+ and FMD neutrino transport, exhibit dipole components that grow with time and reach between 30% and more than 100% of the monopole amplitude at the end of the calculations. While the individual characteristics of the evolution differ from model to model, a clear correlation can be observed in all cases between the time when the dipole amplitude reaches a substantial fraction of the monopole and the time when the direction of the dipole approaches a fairly stable state with only slow and modest subsequent migration. This can be most clearly seen in the SN runs for the s9.0 progenitor, for example around 200 ms p.b. in models ‘s9.0 FMD L’ and ‘s9.0 RbR+ H’, and at 350 ms in ‘s9.0 FMH H’, whereas in model ‘s9.0 RbR+ L’ the latitudinal angle reaches a well-defined value after about 100 ms while the azimuthal angle continues to drift until 350 ms p.b. But this motion takes place in a rather small region around the pole of the grid since (compare Figures 3 and 3).
In the s20 runs the dipole direction also has the tendency to evolve more slowly after an initial phase of basically random motion during the period when the amplitude still grows. However, different from the convection-dominated post-shock layer in the s9.0 models, the post-bounce dynamics of the s20 runs is characterized by repeated episodes of violent SASI activity. As shown by Tamborra et al. (2014a), the LESA dipole is mostly formed in the convective layer interior of the proto-neutron star and further enhanced by the lepton-number loss from the accretion layer around the neutrinosphere, where the electron fraction develops a pronounced hemispheric asymmetry. For this reason Tamborra et al. (2014a) and Tamborra et al. (2014b) observed that SASI-induced variations of the accretion flow between stalled shock and neutron star and the associated time- and direction-dependent modulations of the accretion luminosity (mostly plus ) can have a considerable impact on the evolution of the LESA dipole. This influence of SASI on LESA is strongest when the LESA direction coincides with the plane or even the direction of the SASI sloshing or spiral motions. The interplay between SASI and LESA does not depend on the transport method employed.
The corresponding effects can be clearly spotted in the s20 runs, see Figures 3 (upper panel) and 3, where phases exist when both the dipole amplitude and the direction angles exhibit pronounced oscillations. These oscillations correlate with episodes of strong SASI activity in the post-shock layer (see Section 4.1 and Figure 6 of Paper I). The overlap of SASI and LESA directions is particularly evident in model ‘s20 RbR+ H’, where between about 200 ms and 450 ms after bounce the LESA dipole and the SASI direction are roughly aligned as visible in the upper right panel of Figure 3. Exactly during this time interval the LESA dipole amplitude (Figure 3) shows large quasi-periodic modulations and both direction angles of the LESA dipole exhibit a major shift by 90 instead of the usual temporal fluctuations, which seem to be a consequence of the time- and directional variability of the massive accretion flows to the neutron star in the nonexploding s20 models (Figure 3). While the ubiquitous temporal fluctuations of the direction angles are also present in ‘s20 FMD H’, this case does not display any quasi-monotonic large-scale migration of the LESA direction at ms, except for local motions around the pole at , signalled by large excursions of between 350 ms and 500 ms after bounce. This is compatible with the fact that in model ‘s20 FMD H’ the LESA dipole direction is far off the main plane of SASI activity during most of the time (see upper left panel in Figure 3). For the Aitoff projections of the SASI direction in the latter figure we consider the dipole vector associated with the shock deformation555We use the multipole coefficients , , and from the decomposition of the angle-dependent shock radius into real spherical harmonics (see Equations (11) and (12) and the associated explanations in Paper I) as components of the direction vector, which we then transform into spherical polar coordinates., following Tamborra et al. (2014a). Alternatively, Tamborra et al. (2014b) employed the dipole of the neutrino-energy flux, summed over all six neutrino species. The spatial motions of both direction vectors are closely correlated, and both can therefore be used equally well as tracers of the SASI motions and of the associated accretion-modulated anisotropic neutrino emission, which then has an impact on the LESA dipole.
In contrast, the runs of the s9.0 progenitor explode rather early ( ms) and do not possess any preceding SASI activity. Moreover, convective fluctuations of the mass-accretion rate also have a much reduced influence on the neutrino emission in these simulations because of the low mass-accretion rate of this low-mass progenitor. For these reasons the LESA amplitude (Figure 3, lower panel) and direction angles (Figure 3, two right columns and Figure 3, bottom panels) do not show any phases of quasiperiodic oscillations as in the s20 models, and the LESA vector is more stationary or shifts only slowly on longer timescales.
Figure 3 presents the variations of the local electron lepton-number flux densities in all emission directions for all high-resolution 3D simulations at the times when the runs were terminated. The dominance of a dipolar mode is visible in all cases, but also contributions from higher multipolar models can be recognized. This is confirmed by Figure 3, which provides “spectrograms” for the temporal evolution of the moments from order up to for the same set of simulations. Both of the models with FMD transport exhibit an early phase between p.b. and p.b. in which only the multipole modes around are excited. The growth of these asymmetries coincides with the time when convective activity inside of the proto-neutron star sets in (–50 ms), whereas nonradial instabilities in the gain layer start only some 10 ms later (at 70–80 ms p.b., see the nonradial kinetic energies given in Figure 6 of Paper I). A dominant dipole begins to appear in ‘s20 FMD H’ at ms, at which time ‘s9.0 FMH H’ develops most strength in the quadrupole mode, before after 450 ms of post-bounce evolution also in this model the dipole becomes the clearly dominant mode, as prominently visible in the lower left panel of Figure 3. At the end of both runs the dipole has by far the highest amplitude of all modes with .
Also in the simulations with RbR+ transport the modes of higher orders between and appear first, the quadrupole mode reaches a dominant amplitude several 10 ms afterwards, and the dipole mode even 150–200 ms later. And also similar to the FMD runs, as time goes on there is an obvious tendency that the power in the mode spectrum shifts to lower and lower modes until the dipole has the highest amplitude. Different from the FMD runs, for which the lepton-flux asymmetry at –300 ms is distinctly concentrated in the and modes, the RbR+ models possess a broader distribution of power in the spectrograms. The higher moments, , always yield visible contributions and their power decreases gradually with higher values of . This is fully compatible with the flux variations seen in Figure 3, where the RbR+ cases exhibit a fragmented pattern of larger and smaller spots (i.e., of lower and higher multipole orders), superimposed on the hemispheric (dipolar) asymmetry.
This observation, however, is not really astonishing in view of the fact that the RbR+ approximation ignores nonradial components of the neutrino fluxes and therefore variations in the directions perpendicular to the radius vector remain more localized on smaller scales. In contrast, the nonlocal behavior of radiative transfer with the FMD scheme, which couples spatial volumes over a neutrino mean free path, leads to a “smearing” of angular variations by the nonradial flux components. This produces much smoother radiation characteristics around the neutrinosphere and above it, as discussed in Section 4.1 and as visible in Figures 9 and 10 of Paper I. Therefore large-scale or global directional variations of the neutrino emission such as accretion-induced hemispheric asymmetries and the LESA phenomenon can be captured in their basic properties by both FMD and RbR+ transport treatments, whereas small-scale variations are stronger with RbR+ transport.
A more detailed comparison reveals further differences between the results from the two treatments, though none of these differences is of fundamental nature. For example, the simulations with RbR+ show a somewhat earlier rise of the LESA dipole component for both progenitor models than in the corresponding FMD cases. Moreover, in the s9.0 models the dipole amplitudes reach somewhat higher maximal values in the RbR+ calculations. Although this does not apply to the s20 runs, the SASI-induced modulations of the LESA dipole seem to be more extreme in the RbR+ cases, in particular in model ‘s20 RbR+ H’ (Figure 3). But these modulations depend sensitively on the orientation of the LESA dipole relative to the SASI direction, which is stochastic and differs from model to model and can also vary during the post-bounce evolution of a single run. The migration of the dipole (or its directional stability), although difficult to compare in detail, does not reveal any qualitatively different behavior between RbR+ and FMD simulations.
Overall, we conclude that the LESA phenomenon with its characteristic features is shared by all of our 3D models. This disproves speculations that LESA is a numerical artifact of the RbR+ approximation. Our results are consistent with those reported by Tamborra et al. (2014a), Tamborra et al. (2014b), and Janka et al. (2016), where LESA was discussed as a lepton-number emission self-sustained asymmetry first. Since a key feature of LESA is a dipole emission component that grows relative to the monopole, it shall be noted that during the phases when the dipole dominates the higher-order multipoles, the absolute value of the total electron lepton-number flux (the monopole) in our s20 models is higher by a factor of than in the simulations with the Vertex code presented by Tamborra et al. (2014a). This discrepancy might stem from differences in the neutrino-interaction rates, which are considerably more elaborate in many aspects in the calculations with Vertex. It could also be a consequence of a different contraction behavior of the neutron star in response to different energy and lepton-number loss through neutrino radiation. Moreover, proto-neutron star convection is substantially stronger in calculations with the Alcar code as discussed by Just et al. (2018). All of these aspects can directly or indirectly influence the production and the convective and radiative transport of lepton number inside of the neutron star and can thus have an effect on the relative strength of the monopole and dipole of the electron lepton-number emission.
The detailed temporal evolution of LESA strongly depends on the simulation setup, including the neutrino-transport method, the grid resolution, and the progenitor model. Although the onset of lepton-emission asymmetries is seen at roughly the same time in all of our models, with –6 modes showing the fastest growth, the subsequent evolution that leads to the emergence of a dominant LESA dipole could depend also on stochastic initial seeds and a stochastically-triggered complex feedback mechanism between asymmetric neutron-star convection and mass accretion of lepton-rich material as described by Tamborra et al. (2014a) and Janka et al. (2016). Further analysis is needed, based on a future, larger pool of 3D simulations with varied physics inputs and numerical treatments of hydrodynamics and neutrino transport. In this context it is assuring that O’Connor & Couch (2018b) and Vartanyan et al. (2018) also observed the LESA phenomenon in their 3D simulation with two-moment transport including velocity-dependent terms.
4 LESA and Proto-Neutron Star Convection
A detailed description of the characteristic empirical features of the LESA phenomenon was provided by Tamborra et al. (2014a) and Janka et al. (2016), whose results were confirmed by the analysis of our current models in Section 3 (see also O’Connor & Couch, 2018b; Vartanyan et al., 2018). As mentioned above and discussed in the previous works, the major contribution to the dipole emission originates from the convective layer inside of the neutron star, whose hemispheric differences in the lepton-number loss rate can be further amplified by an asymmetry in the mass-accretion flow from the SN shock to the neutron star. Figure 3 displays a zoom of the interior of the neutron-star and its immediate surroundings in our 3D model ‘s20 FMD H’. The cross-sectional plane chosen for the plots is close to the axis of the LESA dipole.
The upper left panel shows the radial component of the lepton-number flux in the laboratory frame (i.e., the rest frame of the stellar center), (see Equation (3)), scaled with the surface area , at ms, which is about the time when the LESA dipole reaches its maximum (see Figure 3, upper left panel). The higher lepton-number flux in the northern hemisphere is correlated with an excess of the electron fraction, , compared to its angular average at each radius (lower left panel). This well-known LESA pattern grows within the convective layer, which extends from the edge of the white circular region indicating the 10 km core that is computed in 1D, to the outer boundary of the convective layer, which roughly coincides with the effective energy sphere of muon and tau neutrinos (dotted green circle in the lower left panel). The LESA -emission and asymmetries continue within the overlying radiative shell, where faster neutrino diffusion carries the leptons to the neutrinospheres of and (solid and dashed green circles, respectively), and persist into the accretion layer exterior to the neutrinospheres.
The effective flow pattern that is responsible for these hemispheric differences in the lepton-number transport out from the deep inner core of the proto-neutron star is displayed in the right panels of Figure 3, where the component of the mass flux in the north-south () direction, , is visualized by color coding. Superimposed on smaller-scale convective variations (see the pattern of convective cells in the upper left panel), a global dipolar flow can be identified, indicated by the white arrows in the lower right panel: Around the north pole lepton-rich plasma rises through the convective layer, transporting electron-lepton number from the inner core to regions closer to the neutrinospheres, where the high electron abundance causes enhanced emission. Matter with decreasing streams along the outer edge of the convective layer from the north pole towards the south pole, further losing lepton number by radiating . Near the south pole the flow of neutron-rich, specifically (i.e., per nucleon) heavier and thus denser matter submerges again to deeper layers, where it is channelled into currents that stream back towards the northern hemisphere along the inner border of the convective layer. On its way from south to north, the plasma is refuelled with fresh leptons by absorbing neutrinos diffusing out from the high-density core.
A deficiency of the present models in connection to this analysis of LESA should be mentioned here. The 1D core with a chosen radius of 10 km, which is employed in our 2D and 3D simulations for reasons of computational efficiency, turned out to be too large during the later post-bounce evolution of our models. At times beyond about 200 ms after bounce, the inner boundary of the convective layer hits the radius of the 1D core (see upper left panel of Figure 3), and correspondingly, the spherically symmetric treatment of the 1D core prevents the convective shell to grow deeper. The 2D calculations of Just et al. (2018), where the s20 progenitor was run with the same input physics but with a considerably smaller 1D core, demonstrate that the inner boundary of the convective shell reaches a radius of 10 km already after about 150 ms and becomes as small as 7 km after 500 ms. It is therefore clear that the use of the spherical core has constrained the later growth of the proto-neutron star convection layer in the simulations of the present paper. Potentially, this might even have prevented a growth of the LESA dipole beyond the amplitudes obtained in our models. Nevertheless, the basic features of the LESA phenomenon agree well with the results reported by Tamborra et al. (2014a) and Janka et al. (2016) from calculations with the Prometheus-Vertex code, where the 1D core was as small as 1.5 km and had no impact on the development of neutron-star convection. Evaluating the available set of 3D models from the Prometheus-Vertex runs, Stockinger (2015) also observed a dipolar flow pattern that encompasses the convective layer in the neutron star and transports electron lepton-number from one hemisphere to the other, very similar to the one visible in Figure 3.
5 Physics Origin of the LESA Phenomenon
What is the physical mechanism that leads to the LESA lepton-emission asymmetry? Why does the dipole mode become the dominant multipole component after a few 100 ms of post-bounce evolution? In the following we intend to further illuminate the development of convection in the proto-neutron star, where the major contribution to LESA is produced. For this purpose we attempt, tentatively, to establish a link to the theoretical discussion of the onset of thermal instability in spherical shells by Chandrasekhar (1961). Although this cannot provide a satisfactory explanation of the properties of the system in the fully nonlinear stage, we will be able to demonstrate that the different -components in the mode spectrograms of the lepton-number flux (Figure 3) exhibit a growth behavior compatible with basic expectations from Chandrasekhar’s linear theory. Besides the fact that this connection may be considered as a check of physical consistency, our analysis also allows for speculations about the subsequent evolution of the system in view of the trends displayed by its characteristic parameters.
Building up on the seminal work by Lord Rayleigh, Chandrasekhar (1961) determined the critical conditions for the onset of thermal instability in fluid spheres and spherical shells by linear stability analysis of a variety of situations, invoking simplifying assumptions (e.g., the Boussinesq approximation, simple radial dependences of gravity and temperature, and constant medium-dependent coefficients such as those for volume expansion, viscosity, specific heat, or thermal conductivity). He constitutes that, provided the principle of the exchange of stabilities is valid (i.e., that the transition from stability to instability occurs via a stationary pattern of fluid motions instead of oscillatory motions), instability will set in for disturbances of the spherical harmonics mode with the lowest value of the critical Rayleigh number. This means that at marginal stability the convective pattern will manifest itself in this -mode when the Rayleigh number of the system,
exceeds the critical Rayleigh number of the mode. In Equation (7), is the gravitational potential, the temperature, the coefficient of volume expansion, the thermometric conductivity, and the kinematic viscosity. Since in our case of application viscosity is dominated by numerical effects, we use the subscript “N” for the kinematic viscosity. The radii and are the inner and outer boundaries of the convective layer, respectively.
With this definition of the Rayleigh number the critical condition means that instability sets in at the minimum temperature gradient at which a balance can be maintained between the kinematic energy dissipation by viscosity and the internal energy liberated by the buoyancy force (Chandrasekhar, 1961). It further means that for Rayleigh numbers lower than the critical Rayleigh number of a mode , disturbances of this mode are expected to be stable, while such disturbances will become unstable when the Rayleigh number exceeds .
We attempt now to verify these aspects from our numerical results. In Figure 3, left panels, we provide information on the evolution of the outer boundary, , and inner boundary, , of the neutron-star convection layer in our 3D model ‘s20 FMD H’, on the central radius of this layer, , and on the radius ratio . Here, is defined as the radius at which the angle-averaged, non-radial velocities exceed , whereas is determined as the radial position where the angle-averaged, non-radial velocity drops to roughly 40% of its maximum value in the convective shell. The latter condition allows for an unambiguous identification of the outer shell boundary also at times when accretion downflows create substantial velocity perturbations in the convectively stable layer that separates the convective shell inside of the neutron star from the turbulent flows in its surroundings.
The lower left panel of Figure 3 also shows the Rayleigh number , evaluated at , as a function of post-bounce time. is calculated according to Equation (7), using angular averages for all quantities needed. The coefficient of volume expansion, , can be derived from the equation of state,
with , , and being pressure, density, and temperature. For the thermometric conductivity we use, as a proxy666The heat flux , with being the coefficient of heat conduction, can be identified with the total neutrino flux density (accounting suitably for all neutrino species), , where is the diffusion coefficient and the neutrino chemical potential. The latter relation holds under the assumption of equilibrium diffusion, in which case , and for negligible chemical potential gradients. Therefore . Using the definition (Chandrasekhar, 1961) and with being the specific heat (per unit of mass) at constant volume and the internal energy density including the contributions from stellar plasma and neutrinos, one gets: . We assume that the factor is roughly a constant, and thus we can identify , which, for the radial component, yields ., an effective diffusion coefficient for the total energy flux carried by neutrinos (summed over all neutrino species ):
where and are measured in the co-moving frame of the stellar plasma. In the kinematic viscosity appearing in Equation (7), we consider the contributions from neutrino viscosity as well as numerical viscosity associated with the hydrodynamics solver. The latter is usually dominant by a factor of 20–100, depending on the time and position in the convective shell. The neutrino viscosity, , is evaluated by following Keil et al. (1996) and Guilet et al. (2015, see Equation (8) there) and is found to range between some cm s and nearly cm s, in agreement with estimates in the mentioned references. In contrast, the numerical Reynolds number during the relevant time is around 200–250 (employing Equation (6.36) of Melson 2016, in rough agreement with crude estimates by Keil et al. 1996). Using
with km and maximal turbulent velocities in the convective layer of cm s, we obtain cm s. Therefore we get within the convective layer of the proto-neutron star.
The Rayleigh number (Figure 3, lower left panel) exhibits a steep, monotonic rise during the time interval between 60 ms and 130 ms after bounce, when convection in the proto-neutron star of model ‘s20 FMD H’ gains strength and asymmetry modes of the electron lepton-number flux with multipole orders from to develop significant amplitudes for the first time (see Figure 3). For defining the onset times we choose a threshold value for the multipole coefficients of , which corresponds to the level where the initial multipole structures can be identified in Figure 3. The start times as function of are given in the upper right panel of Figure 3. Now picking the Rayleigh numbers at times from the lower left panel of Figure 3, we can plot the black (solid and dashed) lines in the lower right panel of this figure.
These values of can be confronted with the critical Rayleigh numbers derived by Chandrasekhar (1961) for the two shell-radius ratios of and , which are in the ballpark of our model results (see the blue line in the lower left panel of Figure 3). We refer the reader to Chandrasekhar (1961) for explicit expressions and plotted as well as tabulated values. We are only interested in comparing the shapes of the curves and as functions of mode-number , because Chandrasekhar’s linear analysis of thermal instability considered largely different physical conditions adequate for thermal instability in laboratory and terrestrial environments, and because our evaluation of disregards potentially relevant numerical factors (e.g., when using ). Therefore we scale Chandrasekhar’s critical Rayleigh numbers such that the value for coincides with the minimum obtained for our Rayleigh numbers at this multipole mode, indicated by black asterisks in the lower right panel of Figure 3.
We consider two of the four cases of boundary conditions employed by Chandrasekhar (1961), namely on the one hand free surfaces (where tangential viscous stresses are assumed to vanish) at and (case A displayed by dashed colored lines in Figure 3 and shifted by a factor of 3 for better visibility) and on the other hand a rigid surface at (where the nonradial components of the velocity vanish) and a free surface at (case B, shown with solid colored lines in the lower right panel of Figure 3). Although not motivated by the physical conditions in the proto-neutron star, we also picked the second case because the 1D core used in our simulations has a non-negligible influence on the development of the convective layer by damping or even suppressing nonradial fluid motions close to the 3D/1D boundary.
The lower right panel of Figure 3 thus allows us to compare our Rayleigh numbers at the times when the different multipole modes of the lepton-number flux between and begin to grow (black lines) with the (rescaled) critical Rayleigh numbers when disturbances of such -modes are expected to become convectively unstable (blue and orange lines for and , respectively) according to the linear analysis by Chandrasekhar (1961). Interestingly, the shapes of the corresponding curves exhibit good overall resemblance.777The location of the minimum of the critical Rayleigh line depending on the shell thickness, i.e. the value of , can be roughly understood by a simple argument. For a shell thickness the preferred wavelength of a convective cell is . With one obtains a preferred wavelength of , which yields and . On the other hand, the wavelength of an -mode () on a circle of circumference corresponds to . Requesting , one finds that fits better for and matches better . In particular case B for follows the trend of amazingly well, with a minimum in both cases being located at , where a substantial value of can be observed first in Figure 3. The value of is indeed very close to the boundary-radius ratio for the convective layer in time interval , and the better match for case B makes perfect sense in view of our numerical constraints at the inner boundary radius .
These results therefore suggest that the -mode pattern in the electron lepton-number emission, emerging in model ‘s20 FMD H’ around 100 ms p.b. (Figure 3), is compatible with Chandrasekhar’s linear theory for the onset of thermal convection in spherical shells applied to the convective layer in the neutron star. While initially convective modes around are favored in their growth, because their critical Rayleigh numbers are lowest, the monotonic rise of the Rayleigh number within the contracting neutron star also permits dipolar and quadrupolar modes to appear after 100 ms. Since the Rayleigh number increases further by about an order of magnitude, the later dominance of the convective mode (Figure 3) and of the corresponding LESA dipole (Figures 3 and 3) is certainly not in conflict with the critical behavior of thermally unstable shells as discussed by Chandrasekhar (1961).
However, Chandrasekhar’s theoretical assessment does not yield any prediction of the evolution of the convection layer in the nonlinear phase and, moreover, it does not provide any growth rates for the different convective -modes. Therefore, it only allows one to conclude that the conditions at ms post bounce are favorable also for the presence of the mode, but it does not give an argument why the mode indeed becomes the dominant convective and lepton-emission asymmetry during long phases of the evolution. In fact, follow-up work points to the need of a nonlinear treatment to answer the question of the pattern of convection in spherical shells (e.g., Busse, 1975).
Our connection to Chandrasekhar’s analysis can only be considered as tentative, because we invoked a number of shortcuts. Convection inside nascent neutron stars is driven by entropy and lepton-fraction gradients, captured by the Ledoux criterion, and our evaluation in Section 4 highlights the importance of variations of the lepton fraction in the stellar plasma and in the neutrino diffusion fluxes. A rigorous analytic assessment must take this into account and must thus generalize the mathematical treatment by Chandrasekhar (1961), which was focussed on the problem of thermal instability driven by temperature gradients and the associated heat conduction. The critical -mode-dependent Rayleigh numbers for the more complex, more general situation of proto-neutron star convection still have to be derived.
An interesting implication of our discussion is the flashlight on the role of numerical viscosity. With the resolution used in our simulations, it is clearly dominant over the neutrino viscosity by a significant factor (one to two orders of magnitude). Therefore the numerical viscosity determined our estimated values of the Rayleigh numbers for neutron-star conditions. This suggests that higher values of the numerical viscosity might artificially damp or suppress the appearance of dipolar and quadrupolar flows, because the onset of these modes is constrained by higher values of the critical Rayleigh number. This fact should be kept in mind when LESA is searched for in the results from numerical calculations.
Nonlinear effects in the context of LESA beyond Chandrasekhar’s mathematical theory were described by Tamborra et al. (2014a) and Janka et al. (2016). In order to explain the development of a dominant dipolar mode in neutron-star convection and lepton-emission, and the long-time stability of this phenomenon, Tamborra et al. (2014a) proposed a coupling to a hemispheric difference of the accretion flow between shock and neutron star. Because of the harder spectra of compared to , the LESA dipole leads to stronger neutrino heating in the direction of the higher flux, i.e., opposite to the LESA dipole direction. The stronger heating pushes the SN shock to larger radii and channels the mass-accretion flow preferentially to the opposite hemisphere. The higher accretion on this side of the neutron star brings fresh leptons mainly to the hemisphere that is already lepton-rich, and thus it strengthens the LESA dipole emerging from the neutron-star convection layer.
Because LESA is also found in neutron stars from low-mass progenitors, where post-bounce accretion plays a relatively unimportant role, Janka et al. (2016) suggested another feedback mechanism that could be responsible for the steep growth of the LESA dipole. These authors discussed an amplification effect caused by a changing sign of the thermodynamic derivative (with being the entropy per baryon) that appears in association with the gradient of the lepton fraction in the Ledoux criterion for convective instability. The sign changes from negative to positive when the medium deleptonizes, with the threshold value of being higher for higher density and temperature. When drops below the critical value, a negative lepton-fraction gradient becomes stabilizing instead of destabilizing, thus damping convective activity instead of driving it (because only a negative entropy gradient remains as a driving term for the Ledoux convection). As a consequence, the convective flow is attenuated especially in the hemisphere opposite to the LESA direction, where the electron fraction and the total lepton fraction have dropped to lower values (see lower left panel in Figure 3). This weakens the convective transport of lepton number out from the diffusive inner core in this hemisphere, thus enhancing the LESA dipole asymmetry. Such a positive feedback cycle could gradually amplify the asymmetry and therefore might explain why at some point in the evolution the dipole amplitude of LESA increases monotonically. Since the electron fraction of the stellar plasma as a function of density and temperature plays a role (through its influence on the sign of ), this picture might also explain why differences in the neutrino treatment (e.g., FMD vs. RbR+ or different neutrino opacities) as well as the nuclear equation of state (which determines the thermodynamic derivative) can have an impact on the growth of LESA (see Janka et al., 2016, for a report on equation-of-state dependent differences).
We presented an in-depth analysis of the effects of LESA in a set of eight 3D SN runs for two progenitors published in Paper I. LESA was witnessed as a dipolar electron lepton-number emission asymmetry first by Tamborra et al. (2014a) and was recently also seen in one 3D model by O’Connor & Couch (2018b) and, possibly, by Vartanyan et al. (2018). We observed the appearance of a strong LESA dipole mode after some 100 ms p.b. in all of our models, for exploding 9 and nonexploding 20 cases, both with lower and doubled grid resolution, and with fully multi-dimensional (FMD) as well as ray-by-ray-plus (RbR+) neutrino transport, which neglects the nonradial components of the neutrino fluxes.
All characteristic features of LESA reported by Tamborra et al. (2014a) and Tamborra et al. (2014b) were also diagnosed in the present models: (1) phases where the dipole component of the lepton-number flux clearly dominates all higher-order multipoles and can reach even the strength of the monopole; (2) a relatively stable position or only slow migration of the dipole direction during the phases of a strong dipole; (3) corresponding hemispheric differences in the strength of proto-neutron star convection and in the electron and lepton fractions in the convection layer and above it; (4) a hemispheric asymmetry of the accretion layer between deformed shock and neutron star with greater accretion on the side of the higher lepton-number flux; (5) SASI-induced modulations of the LESA dipole in periods of strong SASI activity, when the LESA dipole is located close to the direction of SASI sloshing or close to the plane of SASI spiral motions.
Our results clearly prove that LESA is not an artifact of the RbR+ approximation. Nevertheless, there are interesting differences between the FMD and RbR+ results in certain aspects. The LESA dipole amplitudes with RbR+ neutrino transport tend to grow noticeably faster, reach higher amplitudes earlier, and their maximum amplitudes can, at least temporarily, be higher than those obtained with FMD transport. However, the dipole, once fully developed, dominates the multipole spectra of the lepton-number flux in the FMD cases in a much cleaner way than in the RbR+ runs (although it can be somewhat lower compared to the monopole). The spectrograms of the RbR+ models typically show prominent quadrupoles as well, and significant power over a wider range of higher-order multipoles. This patchiness of the directional variations of the lepton-number emission is expected as a natural consequence of the more localized neutrino effects in the RbR+ approximation (see Paper I). In contrast, the neutrino field leaving the convective layer is smoothened in the surrounding radiative shell by non-local radiation-transport effects of FMD.
Closely inspecting the mass flows in and around the convection layer interior to the proto-neutron star during phases of a strong LESA dipole, we reveal a global, volume-filling dipolar pattern that streams around the central, convectively stable core. These currents effectively transport lepton-rich matter from the hemisphere opposite to the LESA dipole direction to the hemisphere in the LESA dipole direction, and from deeper regions to layers closer to the neutrinospheres. An important effect is the gain or loss of electron lepton-number by these currents through exchange with diffusion fluxes (via neutrino absorption or emission) when the gas flow moves along the inner and outer boundaries between the convective layer and the diffusive central core on the one side and the enveloping shell on the other side. This exchange process fuels and defuels the lepton content of the currents.
The presence of a dipolar flow pattern through and around the convection zone of the neutron star motivated us to investigate whether this finding is compatible with predictions by Chandrasekhar’s linear stability analysis of thermal instability in spherical shells (Chandrasekhar, 1961). Our assessment can, of course, only be tentative because of the substantial differences between the complex physical conditions in the Ledoux-unstable neutron-star layer compared to the idealized terrestrial setups investigated by Chandrasekhar (1961). Despite this caveat the growth behavior of the different multipoles of the lepton-number flux, which we relate to multipoles of the flow in the convective layer, follows amazingly well the expectations from Chandrasekhar’s theory. We conclude this from comparing the Rayleigh numbers of the system, which we estimate at the center of the convective layer, at the times when we identify the initial rise of an -mode emission asymmetry, with the critical Rayleigh number for this -mode derived by Chandrasekhar. His theory predicts that an -mode perturbation should reach marginal stability when the Rayleigh number grows to the critical value of the mode. In agreement with the mathematical analysis we witness the onset of modes around first, and the growth of the other multipole modes setting in at times when the Rayleigh numbers for the conditions of the system follow the critical curve adopted from Chandrasekhar (1961). Since the Rayleigh number of the convective layer continues to rise monotonically also after the dipole mode, , appears for the first time, one can expect its presence also at later times.
Calculating the Rayleigh number requires the specification of the kinematic viscosity of the environment. In current models viscous effects are far dominated (by a factor 10–100) by numerical viscosity. This needs to be kept in mind when numerical simulations are evaluated for the growth of LESA multipole modes.
Our reported work has attempted to achieve progress towards a better understanding of the puzzling LESA phenomenon. Despite the assuring result that the existence of a dipolar flow pattern in the neutron-star convection layer seems compatible with Chandrasekhar’s theory, the detailed reason for the dominance of the mode remains unclear. Previous speculations that it is driven by feedback effects with a hemispheric accretion asymmetry around the neutron star (Tamborra et al., 2014a), or by an amplifying feedback cycle connected to -dependent differences of Ledoux convection in opposite neutron star hemispheres (Janka et al., 2016), require corroboration by future investigations.
We thank Irene Tamborra for her valuable comments on the manuscript. At Garching, this project was supported by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA, and by the Deutsche Forschungsgemeinschaft through Sonderforschungbereich SFB 1258 “Neutrinos and Dark Matter in Astro- and Particle Physics” (NDM) and the Excellence Cluster Universe (EXC 153; http://www.universe-cluster.de/). OJ acknowledges support by the Special Postdoctoral Researchers (SPDR) program and the iTHEMS cluster at RIKEN. Computer resources for this project have been provided by the Leibniz Supercomputing Centre (LRZ) under grant pr62za, and by the Max Planck Computing and Data Facility (MPCDF) on the HPC system Hydra.
- Busse  Busse, F. H. 1975, Journal of Fluid Mechanics, 72, 67
- Chandrasekhar  Chandrasekhar, S. 1961, International Series of Monographs on Physics, Oxford: Clarendon, 1961
- Dolence et al.  Dolence, J. C., Burrows, A., & Zhang, W. 2015, ApJ, 800, 10
- Glas et al.  Glas, R., Just, O., Janka, H.-T., Obergaulinger, M. 2018, submitted to ApJ (Paper I)
- Guilet et al.  Guilet, J., Müller, E., & Janka, H.-T. 2015, MNRAS, 447, 3992
- Hunter  Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
- Janka et al.  Janka, H.-T., Melson, T., & Summa, A. 2016, Annual Review of Nuclear and Particle Science, 66, 341
- Just et al.  Just, O., Obergaulinger, M., & Janka, H.-T. 2015, MNRAS, 453, 3386
- Just et al.  Just, O., Bollig, R., Janka, H.-T., et al. 2018, arXiv:1805.03953
- Keil et al.  Keil, W., Janka, H.-T., & Mueller, E. 1996, ApJL, 473, L111
- Marek et al.  Marek, A., Dimmelmeier, H., Janka, H.-T., Müller, E., & Buras, R. 2006, A&A, 445, 273
- Melson  Melson, T. 2016, PhD thesis, Technische Universität München
- Obergaulinger  Obergaulinger, M. 2008, PhD thesis, Technische Universität München
- O’Connor & Couch [2018b] O’Connor, E., & Couch, S. 2018b, arXiv:1807.07579
- Oliphant  Oliphant, T. E. 2007, Computing in Science and Engg., 9, 10
- Perez & Granger  Perez, F., & Granger, B. E. 2007, Computing in Science and Engg., 9, 21
- Roberts et al.  Roberts, L. F., Ott, C. D., Haas, R., et al. 2016, ApJ, 831, 98
- Steiner et al.  Steiner, A. W., Hempel, M., & Fischer, T. 2013, ApJ, 774, 17
- Stockinger  Stockinger, G. 2015, Master thesis, Technische Universität München
- Sumiyoshi et al.  Sumiyoshi, K., Takiwaki, T., Matsufuru, H., & Yamada, S. 2015, ApJS, 216, 5
- Tamborra et al. [2014a] Tamborra, I., Hanke, F., Janka, H.-T., et al. 2014a, ApJ, 792, 96
- Tamborra et al. [2014b] Tamborra, I., Raffelt, G., Hanke, F., Janka, H.-T., & Müller, B. 2014b, PhRvD, 90, 045032
- Vartanyan et al.  Vartanyan, D., Burrows, A., Radice, D., Skinner, M. A., & Dolence, J. 2018, MNRAS, 477, 3091
- Vartanyan et al.  Vartanyan, D., Burrows, A., Radice, D., Skinner, A., & Dolence, J. 2018, arXiv:1809.05106
- Woosley & Heger  Woosley, S. E., & Heger, A. 2007, PhR, 442, 269
- Woosley & Heger  Woosley, S. E., & Heger, A. 2015, ApJ, 810, 34