One-, Two-, and Three-dimensional Simulations of Oxygen Shell Burning Just Before the Core-Collapse of Massive Stars
We perform two- (2D) and three-dimensional (3D) hydrodynamics simulations of convective oxygen shell-burning that takes place deep inside a massive progenitor star of a core-collapse supernova. Using one dimensional (1D) stellar evolution code, we first calculate the evolution of massive stars with the initial mass of 9–40 . Four different overshoot parameters are applied, and similar CO core mass trend with previous works is obtained in the 1D models. Selecting eleven 1D models that have a silicon and oxygen coexising layer, we perform 2D hydrodynamics simulations of the evolution for 100 s until the onset of core-collapse. We find that strong and large-scale convection with the Mach number 0.1 is obtained in the models having a Si/O layer with the scale of 10 cm, whereas most models that cast an extended O/Si layer up to a few cm do not exhibit strong convective activity. Our results indicate that the supernova progenitors that possess a thick Si/O layer could provide a preferable condition for perturbation-aided explosions. We perform 3D simulation of a 25 model, which exhibits strong convective activity in the 2D models. The 3D model develops large-scale () convection similar to the 2D model, however, the 3D convection is weaker than the 2D convection. By estimating the neutrino emission properties of the 3D model, we point out that a time modulation of the event rates, if observed in KamLAND and Hyper-Kamiokande, would provide an important information about structural changes in the presupernova convective layer.
0000-0002-8967-7063]Takashi Yoshida \move@AU\move@AF\@affiliationDepartment of Astronomy, Graduate School of Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
0000-0003-0304-9283]Tomoya Takiwaki \move@AU\move@AF\@affiliationDivision of Theoretical Astronomy, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
0000-0003-2456-6183]Kei Kotake \move@AU\move@AF\@affiliationDepartment of Applied Physics, Fukuoka University, Fukuoka 814-0180, Japan
0000-0002-6705-6303]Koh Takahashi \move@AU\move@AF\@affiliationArgelander-Institute für Astronomie, Universitäte Bonn, D-53121 Bonn, Germany
0000-0002-8734-2147]Ko Nakamura \move@AU\move@AF\@affiliationDepartment of Applied Physics, Fukuoka University, Fukuoka 814-0180, Japan
Department of Astronomy, Graduate School of Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
From theory and observations, it is almost certain that the explosions of massive stars as core-collapse supernovae (CCSNe) are generically multi-dimensional (multi-D) phenomena (see foglizzo_15; janka16_rev; patat_book for reviews). To faciliate the neutrino-driven mechanism of CCSNe (bethe85), multi-D hydrodynamics instabilities such as neutrino-driven convection and the standing accretion shock instability (Blondin03) play a pivotal role in enhancing the neutrino heating efficiency to trigger the onset of the explosion. In fact, a growing number of self-consistent models in two or three spatial dimensions (2D, 3D) now report revival of the stalled bounce shock into explosion for a wide mass range of progenitors (see, e.g., vartanyan19; oconnor18b; bernhard17; Roberts16; Nakamura16; melson15a; Summa16; Lentz15; Takiwaki14; Hanke13 for collective references therein).
These successes, however, provide further motivation for exploring missing ingredients in the neutrino mechanism, partly because the estimated explosion energies obtained in the multi-D models generally do not reach the typical observed value (e.g., 10erg, Tanaka09). Various possible candidates to obtain more robust explosions have recently been proposed, including multi-D effects during the final stage of the presupernova evolution (see couch_book for a review), general relativity (GR, e.g., BMuller12a; Ott13; KurodaT12; KurodaT16), rapid rotation (e.g., Marek09; Suwa10; Takiwaki16; Summa18; harada18) and/or magnetic fields (e.g., martin06; moesta14; jerome15; masada15; martin17), and sophistication in the neutrino opacities (melson15b; bollig17; Burrows18; kotake18) and in the neutrino transport schemes (e.g., sumi12; richers17; nagakura18; just18). In this work, we focus on the first item listed in the above list.
couch_ott13 were the first to demonstrate that the inhomogenities seeded by convective shell burning fosters the onset of the neutrino-driven explosion (see also rodorigo14; couch_ott15; bernhard15; Burrows18). This is because the infalling perturbation that could be amplified in the supersonic accretion (takahashi14; nagakura13; nagakura19) enhances turbulence behind the postshock material leading to the reduction of the critical neutrino luminosity for shock revival (e.g., bernhard15; ernazar16). In these studies, the non-spherical structures in the burning shells, albeit physically motivated, were treated in a parametric manner due to the paucity of the multi-D stellar evolution models covering the whole lifespan of massive stars up to the iron core-collapse. Currently one-dimensional (1D) stellar evolution calculations are the only way to accomplish this (woosley02; Woosley07; Sukhbold18), where multi-D hydrodynamics effects are taken into accout by means of the mixing length theory (MLT, e.g., kippenhahn12).
The truely multi-D hydrodynamics stellar evolution calculations have been done generally over several turnover timescales of convection (limited by the affordable computational resources) in selected burning shells (e.g., viallet13; simon15; cristini17 for different burning shells, and see arnett16 for a review). Pushed by the observation of SN1987A, 2D and 3D stellar evolution simulations focusing on the late burning stages have been extensively carried out since the 1990s (arnett94; bazan94; bazan98; asida00; kuhlen03; meakin06; meakin07; arnett11; chatz14; chatz16; jones17).
More recently, ground-breaking attempts to evolve convective shells in 3D prior to the onset of collapse have been firstly reported by Couch15 for silicon shell burning in a 15 star and by bernhard16_prog for oxygen shell burning in an 18 star (and also in 11.8, 12, and 12.5 stars by bernhard18_prog). Couch15 obtained earlier onset of a neutrino-driven explosion for the 3D progenitor model of the star compared to that in the corresponding 1D progenitor. By performing 3D GR simulations with more advanced neutrino transport scheme, bernhard17 obtained a neutrino-driven explosion with the seed perturbations, whereas the shock was not revived in the corresponding 1D progenitor model. These studies clearly show that convective seed perturbations could potentially have a favorable impact on the neutrino-driven explosions. In order to clarify the criteria to make grow the precollapse seed perturbations, Collins18 recently reported a detailed analysis on the convective oxygen and silicon burning shells by performing a broad range of 1D presupernova calculations. Using the prescription of the MLT theory in 1D, they pointed out that the extended oxygen burning shells between 16 and 26 are most likely to exhibit large-scale convective overturn with high convective Mach numbers, leading to the most favorable condition of perturbation-aided explosions. In fact the 3D progenitor model of the 18 star (bernhard16_prog) is exactly in the predicted mass range.
Joining in these efforts, we investigate in this study how the asphericities could grow, particularly driven by the convective oxygen shell burning in the O and Si-rich layer. First we perform a series of 1D stellar evolution calculation with the zero-age main sequence (ZAMS) masses between 9 and 40 with the HOSHI code developed by Takahashi16; Takahashi18. Based on the 1D results, we select ten 1D progenitors that cast extended and enriched O and Si layers, presumably leading to vigorous convection. At a time of s before the onset of collapse, the 1D evolution models are mapped to multi-D hydrodynamics code (a branch of 3DnSNe, e.g., Takiwaki16; Nakamura16; kotake18). We perform axisymmetric (2D) simulations for the selected progenitors having an extended O and Si-rich layer and investigate their convection activities. We then move on to perform a 3D simulation by choosing one of the progenitors that exhibits strong convective activity in 2D. We make an analysis to investigate how the convective activity between 3D and 2D differs and discuss its possible implication to the explosion dynamics.
The paper is organized as follows. Section 2 starts with a brief description of the numerical methods employed in our 1D stellar evolution calculation as well as 2D and 3D hydrodynamics simulations. In Section 3, we present the results of the 1D stellar evolution models in Section 3.1, which is followed in order by 2D (Section 3.2) and 3D (Section 3.3) results, respectively. In Section 4, we summarize with a discussion of the possible implications. Appendices address the comparison of our 1D stellar evolution code with other reference codes and the sensitivity of our results with respect to the different parameters.
2 SETUP and NUMERICAL METHODs
2.1 1D Stellar Evolution
Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefStellar evolution models and the corresponding overshoot parameters. is the overshoot parameter during the H and He burning. is the one after the He burning. See equation (1) and text for details.
We calculate the 1D evolution of solar-metallicity massive stars with the ZAMS masses between 9 and 40 up to the onset of collapse of the iron core. The calculations are performed using an up-to-date version of the 1D stellar evolution code, HOSHI (HOngo Stellar Hydrodynamics Investigator) code333”HOSHI” is a noun in Japanese meaning a star. (e.g. Takahashi16; Takahashi18). In the code, a 300-species nuclear reaction network444The employed nuclear species are as follows: n, H, He, Li, Be, B, C, N, O, F, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, K, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Ga, Ge, As, Se, Br. The isomeric state of Al is also included. is included, the rates of which are taken from JINA REACLIB v1 (Cyburt10) except for C(O (see Takahashi16, for details).
The mass of the helium (He), carbon-oxygen (CO), and iron (Fe) cores as well as the advanced stage evolution depend on the treatment of convection. We use the Ledoux criterion for convective instability. Inside the convective region, we treat the chemical mixing by means of the MLT using the diffusion coefficients as described in Takahashi18.
In order to take into account the chemical mixing by convective overshoot, an exponentially decaying coefficient,
is included, where , , , are the diffusion coefficient at the convective boundary, the distance from the convective boundary, the overshoot parameter, and the pressure scale-height at the convective boundary, respectively (e.g., Takahashi16). We consider the following four models to see the impacts of the different overshoot parameter () on the 1D stellar evolution. First we consider the two cases ( = 0.03 or 0.01, see Table 1) during the H and He burning phases. The former and the latter values are determined based on the calibrations to early B-type stars in the Large Magellanic Cloud (LMC) (Brott11) and the main-sequence width observed for AB stars in open clusters of the Milky-Way Galaxy (Maeder89; Ekstroem12), respectively. We name the former model as model “L” after the LMC and the latter model as model ”M” after the Milky-Way Galaxy. In addition, we change the overshoot parameter from the carbon burning ( = 0 or 0.002, Table 1) to investigate the impact on the advanced stages. Note in Table 1 that the subscript “ov,c” is added to the models with .
The mass loss rate at different evolution stages is important for determining the final mass and the He and CO core masses for high mass stars. We adopt Vink01 as the mass loss rate of a main sequence star where the effective temperature ls higher than and the surface hydrogen mass fraction is higher than or equal to 0.3. The mass loss rate of Nugis00 is adopted for Wolf-Rayet (WR) stars where the effective temperature is higher than and the surface hydrogen mass fraction is lower than 0.3. When the surface temperature is lower than , we adopt Nieuwenhuijzen90 as the red supergiant mass loss rate.
2.2 Multi-D Stellar Hydrodynamics Simulations
We compute 2D and 3D models with our hydrodynamics code of 3DnSEV (3 Dimensional nuclear hydrodynamic simulation code for Stellar EVolution), that is a branch of 3DnSNe code (see Takiwaki16; Nakamura16; Sasaki17; kotake18 for the recent code development). Similar to the base code, 3DnSEV solves Newtonian hydrodynamics equations using spherical polar coordinates. A piecewise linear method with the geometrical correction of the spherical coordinates is used to reconstruct variables at the cell edge, where a modified van Leer limiter is employed to satisfy the condition of total variation diminishing (TVD) (Mignone14). The numerical flux is basically calculated by HLLC solver (Toro94). For the numerical flux of isotopes, the consistent multi-fluid advection method of Plewa99 is used. The models are computed on a spherical polar coordinate grid with an resolution of (3D) and (2D) zones. The radial grid is logarithmically spaced and that covers from the center up to the outer boundary of cm. For the polar and azimuthal angle, the grid covers all steradian. To focus on the convective activity mainly in the oxygen shell, the inner 1000 km is solved in spherical symmetry. We include self-gravity assuming a spherically symmetric (monopole) gravitational potential. Such a treatment is indispensable for reducing the computational time, whereas the non-linear coupling between the core and the surrounding shells (e.g., fuller15) s beyond the scope of this study.
We use the “Helmhotlz” equation of state (EOS: Timmes00). The neutrino cooling is taken into account (Itoh96) as a sink term in the energy equation. A nuclear reaction network of 21-isotopes (aprox21)555http://cococubed.asu.edu/code˙pages/burn˙helium.shtml (Paxton11) is implemented, where the inclusion of Fe,Fe and Cr is crucial to treat a low electron fraction in the presupernova stage. The network is as large as that of Couch15 and a little larger than 19-isotopes of bernhard16_prog. When the temperature is higher than K, the chemical composition is assumed to be in the nuclear statistical equilibrium (NSE). For numerical stability, we set an artificial upper bound in our multi-D runs, in such a way that the (absolute) sum of the local energy generation rates by thermonuclear reactions and weak interactions does not exceed the value of 100 times as large as the local neutrino cooling rate. For correctly treating the neutronization of heavy elements from Si to the iron group and the gradual shift of the nuclear abundances, one needs to use a sufficient number of isotopes (, arnett11), which is currently computationally and technically very challenging. Since the NSE region appears mainly in the Fe core, this treatment may not significantly affect our results in which we focus on convection in the O layer.
At a time of 100 s before the onset of collapse, the 1D evolution models are mapped to our multi-D hydrodynamics code and we follow the 2D and 3D evolution for the 100 s until the onset of collapse. We terminate the 2D/3D runs by a criteria that the central temperature exceeds K, because the core is dynamically collapsing at this time.
3.1 1D Stellar Evolution Models
In total 100 stellar evolution models are calculated in 1D. Four sets of models are constructed, to which different overshoot parameters are applied (Table.1). Each set consists of 25 models to cover the initial mass range of 9–40 M. These 1D models are evolved from the ZAMS stage up to the onset of collapse, which is determined using the threshold central temperature of K. Details of our 1D evolution models (e.g., comparison with reference stellar evolution codes) are given in Appendix.
The top two panels in Figure 3.1 shows the total mass (red) and masses of the He (blue), CO (green), and Fe (magenta and cyan) cores at the onset of collapse as a function of the ZAMS mass (). The top panel is assigned to show results of models L and L, while results of models M and M are shown by the second panel. In the bottom panel, results obtained for models L and M are compared. We note that models L and L result in very similar total, He-core, and CO-core masses, and differ only in the Fe-core mass. This is because the same overshoot parameter is applied until the end of the He burning phase for both cases. The same is true for the models M and M. The He core mass is defined as the largest enclosed mass having the hydrogen mass fraction less than . Similarly, the CO core mass is defined as the largest enclosed mass with the He mass fraction less than 0.1, and the Fe-core mass is defined as the largest enclosed mass with the sum of the mass fractions of elements is larger than 0.5.
The total mass at the collapse is determined by the mass loss history. Since the mass loss is relatively weak, the total mass monotonically increases with the ZAMS mass for models below – in both model L and M. The mass loss rate increases with increasing the luminosity, thus with increasing the ZAMS mass. The increasing mass loss rate explains the flat and even decreasing trends seen in 20–30 models of model L. At the same time, models in the same mass range show a stochastic trend for model M. This is caused by the bistability jump of the mass loss rate, which results in the discontinuous rate increase along the decreasing effective temperature. For more massive models above , the mass loss rate becomes so efficient to remove most of the H envelope during the He burning phase. Therefore, the total masses of these models coincide with their He core masses. This is why the total mass again shows a monotonic increase in this massive end of the ZAMS mass range. The most massive models (32, 35, and 38 models of model L and 35, 38, and 40 models of model M) finally retain only small amount of hydrogen of 0.26–0.29 in their envelopes, which will correspond to be observed as late type WN stars (Crowther07).
As an exception, model 40L (model L with ) has lost not only whole the H envelope but also most of the He layer. This is due to the even stronger WR wind mass loss during the helium and carbon burning phases. The He mass remained on the surface is 0.24 . We apply the mass loss rate of Nugis00 for the H-deficient stars. However, there is a large uncertainty in estimation of the WR wind mass loss rate. Especially, among the H-deficient stars, it has been discussed that the mass loss rate of He-deficient WC stars can be larger than the rate of Nugis00 by a factor of 10 (Yoon2017). The remained He mass can be even less if we consider more efficient WR wind mass loss. Therefore, we expect that the star will be most probably observed as a He-deficient WC star, and moreover, it will be observed as a Type Ic SN when this star explodes as a SN.
Such an overall feature is also seen in the results obtained by previous works using Kepler code (Woosley07) (e.g., Figure 4 of Ebinger18 for a concise summary). Meanwhile, the whole H envelope (and most of the He layer) of stars in the massive mass range of was stripped in our previous work (Yoshida14). Changing the mass loss rate from deJager88 to Nieuwenhuijzen90 accounts for this difference.
The He and CO core masses monotonically increase with the ZAMS mass except for the model 40L, which is affected by the strong WR wind during the He-burning phase. The mass of the helium layer, which is shown as the difference between the He-core and the CO-core masses, also increases with the ZAMS mass. As mentioned earlier, the He and CO core masses are insensitive to the overshoot parameter after the He burning, . Thus, the difference in the CO core mass between models L and L (and similarly between models M and M) is less than 0.7%. Furthermore, the difference in the He core mass is more than one order of magnitude less.
In the bottom panel of Figure 3.1 the ratios of the total, He-core, and CO-core masses between models L and M are shown. Below , the final mass is not so sensitive to the overshoot parameter for H and He core convection. On the other hand, models above 25 show a scatter within a factor of 0.3. Model L tends to show larger He and CO core masses than model M. This is simply due to the larger overshoot parameter applied during the H and He burning phases. For most of the models, the ratio of the He core masses are about a factor of 1.2–1.3, and the CO-core mass ratio is somewhat larger than the He-core mass ratio. As an exception, the He core mass ratio reaches 2.16 for the 9 models. This is due to the merging of He layer to the H envelope by the second dredge-up. Larger overshoot for model L brings about more effective convective mixing to make more massive He and CO cores. Small core mass ratios for 40 models are due to the strong mass loss occurred in the model 40L.
To select models that will show the strong convective activities in the SiO-rich layer, we utilize two measures, convection-width parameters of and . Both of them are defined based on mass fraction distributions of O and/or Si. The is basically a product of the mass fractions of O and Si weighted by the enclosed mass between cm and cm:
where is a scaling coefficient, is the mass fraction of isotope , and is the step function, which satisfies for and 0 for . Therefore, the value becomes large in a model that has a layer mainly composed of both oxygen and silicon. Such a layer would be the most preferable cite to host strong convection powered by oxygen shell-burning. This definition of has an uncertainty on what power of the local density strength of the convection depends. Hence, we also test another indicator, , in which the product of the mass fraction is weighted not by the enclosed mass but by the enclosed volume instead:
where is a scaling coefficient, is the radius in units of cm. The scaling coefficients are arbitrarily chosen. We calculate these two measures at every time steps from 120 s to 80 s before the last step of the calculations.
The result is shown in Figure 3.1, in which and are applied. We do not see clear dependencies among different treatments of overshoot. Some models in the ZAMS mass range show large () values. In the volume weighted case, the ZAMS mass range showing the models having large () values is 13–28 . From this result, we have selected 11 models, in which either or , or possibly both of them, shows a large value. Models showing the seven highest values are models 28M, 23L, 25M, 28L, 27L, 27M, and 22L. Models showing the six highest values are models 13L, 28M, 21M, 25M, 16M, 18M. Among them, models 25M and 28M show large values for both f and . The actual values of the parameters are shown in the second and third columns of Table 3.1.
For later convenience, we separate the SiO-rich layer into the “Si/O” layer and the “O/Si” layer. The “Si/O” layer has larger Si mass fraction than O mass fraction in the layer, i.e., and . Whereas the “O/Si” layer has the relation of . Then, we may classify these 11 models into two groups having different structures of the SiO-rich layer. One group has an extended O/Si layer instead of the O/Ne layer above the Si/Fe layer. The other group has a Si/O layer between the inner Si/Fe layer and the outer O/Ne layer. The former group consists of models 13L, 16M, 18M, 21M, 23L, 27L, and 28M. The top panel of Figure 3.1 shows the mass fraction distribution of model 13L as a function of radius. The radius of the outer boundary of the O/Si layer is cm. This layer was originally formed as an O/Ne layer. Neon burning has started after the core silicon burning phase, transforming neon into oxygen and silicon. In case of models 18M, 21M and 23L, a thin Si/O layer exists between the Si/Fe layer and the O/Si layer with a width of less than 3 cm.
Models in the latter group have the layered structure, in which the innermost Fe core is surrounded by the Si/Fe, Si/O, and O/Ne layers. Models 22L, 25M, 27M, and 28L comprise this group. The bottom panel of Figure 3.1 shows the mass fraction distribution of model 25M, an example of this latter group. The model has the Si/O layer between cm and cm. For other models in this group, the width of the Si/O layer is typically several times cm.
3.2 2D Stellar Hydrodynamics Simulations
In order to investigate the convective activities in a multi-dimensional space, we perform 2D hydrodynamics simulations of oxygen shell-burning. In the previous subsection, we picked up 11 models that show large convection-width parameters, and/or . Profiles of these models at 100 s before the end of the 1D calculations are taken as the initial conditions. The 2D calculations are proceeded until the central temperature reaches K, by which point the stars have started runaway collapse due to the gravitational instability.
Following bernhard16_prog, the strength of the convection is quantified by the angle-averaged Mach number,
where is the density, , , and are the radial, tangential, and azimuthal velocities, is the angle-averaged radial velocity, is the sound velocity, and is the solid angle. The maxima of evaluated at the end of the simulations are shown in the 4th column of Table 3.1, and represents the radii where are obtained. Based on the Mach number, we divide our 2D models into two groups, either showing “weak convection” or “strong convection”. The criterion of strong convection is set as , because convection with such a high Mach number potentially fosters the perturbation-aided explosion (bernhard15; bernhard16_prog). It is noted that model 21is exceptionally classified into weak convection in spite of its large Mach number. We discuss this later in this Section. The column “Layer” in Table 3.1 represents the dominant chemical composition in the convective layer, i.e., the Si/O layer or the O/Si layer. See the last part of the previous subsection for the definition of the layer.
Time evolution of convective activity
In Figure 3.2, we show the time evolution of the convective activity for representative models from weak (13L, top) and strong convection (25M, middle, and 27L, bottom). The color visualizes the angle-averaged convective Mach number (left) and the Si mass fraction (Si) (right). Note that the outer radial frame of the panels for model 27L is set to be cm, in order to show how the outer edge of the convective region keeps moving outward and reaches this radius at the end.
The model 13L has no Si/O layer. This star has an Fe core at the central region of cm, which is surrounded by the convective Si/Fe layer (2–4 cm) and the convective O/Si layer (4 cm). The Si mass fraction at Si/Fe layer is 0.5 (see the top panel of Figure 3.1). As shown in the top right panel of Figure 3.2, the Si mass fraction is small compared to that of 25M model that is shown in the middle right panel of Figure 3.2. Reflecting this structure, the convection developed in the model is weak ( ) in the inner Si/Fe layer and in the outer O/Si layer throughout the simulation (see the top left panel of Figure 3.2). However, the oxygen burning slightly enhances the Si mass fraction in the base region of the O/Si layer of 4–8 cm. Note that in Table 3.1 is estimated at the end of the simulations, and it does not refer to the peak seen at cm at 30 sec.
The features of the convection observed in model 13L are basically common to other models of weak convection (i.e., models of 13L, 16M, 18M, 21M, and 23L), in which a very weak convection develops in the O/Si layer. Besides, either no Si/O layer has been formed (13L and 16M), or only a thin Si/O layer has been formed (18M, 21M, and 23L). Thus, they do not have a strongly convective Si/O layer as well. Models 18M and 21M also have the O/Si layer on the Si/O layer but the maximum Mach number in the O/Si layer is smaller than that in the Si/O layer. The maximum Mach number in the O/Si layer of model 23L is larger than that in the Si/O layer but is still smaller than 0.1.
Readers may be confused on the models of 18M and 21M since they have Si/O layer in Table 3.1 but they are classified into the weak convection group. Actually they have a thin Si/O layer and develop marginally strong shell-convection of in the Si/O layer, but only after the last 10 s of the simulations. This is because the convection is triggered by the gravitational contraction, which amplifies the temperature at the bottom of the Si/O layer enhancing the oxygen burning rate. The shell-convection powered by the gravitational contraction has too short time to form an extended convective region, which is contrasting to the shell-convection powered by a hydrostatic burning. This is why we have selected these models as members of weak convection.
Models 22L, 25M, 27L, 27M, 28L, and 28M are categorized into models with strong convection. Strong convection powered by oxygen burning develops in the Si/O layer in these models. In these models, we pick up two models to explain typical dynamics of the convection.
Model 25M consists of the central Fe core ( cm), the Si/Fe layer (– cm), the Si/O layer (– cm), and the O/Ne layer ( cm). It is noteworthy that, despite the Si/Fe layer seems to have a homogeneous chemical composition (see the middle right panel of Figure 3.2), the outer part of – cm is actually composed of small amount of oxygen with X(O). The oxygen-free region of – cm becomes convective within a short timescale of 10 s from the start of the simulation, though the shell-convection is not extended further. At the bottom of the outer Si/O layer, the hydrostatic oxygen shell-burning takes place. In this case, the nuclear burning drives strong convection with , which is sustained for 20–110 s (see the middle left panel of Figure 4). Accordingly, convective mixing homogenizes the Si mass fraction in the region of 3–10 cm. Furthermore, oxygen burning also takes place in the oxygen-containing outer region of the Si/Fe layer. In spite of the large mean molecular weight, the heating due to the oxygen burning is strong enough to lift the silicon-rich material up into the surrounding Si/O layer. As a result, the silicon mass fraction in the Si/O layer significantly enhances. This silicon enhancement repetitively takes place at 50 and 80 s, which accompanies the enhancement of the convective Mach number as well. It seems that the repetitive mixing follows the oscillation of the outer edge of the Si/Fe convection. This will be because, having the small oxygen mass fraction, a slightly higher temperature than in the case of the Si/O layer is needed to drive the convection in the Si/Fe layer, which is explained by the temperature fluctuation during this oscillation.
In model 27L, the convective activity in the O/Si layer starts to increase at 45 s (see the bottom left panel of Figure 4). The O/Si convective region keeps expanding during the simulation time of 200 s. At the same time, the convective Mach number grows with time in almost whole region in the convective layer. As a consequence, the convective Mach number exceeds 0.15 in the wide outer region of – cm at the end of the simulation. Initially, the Si mass fraction has a decreasing distribution in the Si/O layer. The convection powered by the oxygen burning mixes material in the slightly silicon-enriched region, which is initially located below cm, into the outer slightly silicon-poor region after 70 s. The convective mixing in the 2D simulation is still not efficient enough to achieve the homogeneous chemical distribution. However, this would be due to the limitation of the calculation time. Model 28M shows similar convective properties to the model 27L.
Figure 3.2 shows the 2D distributions of the radial convection Mach number, , and the Si mass fraction taken at the last step of the simulations for models 13L, 25M and 27L. Model 13L (top panels in Figure 5) develops only weak convection of 0.01. The spherical boundary is clearly observed at cm, where the convection Mach number becomes almost zero. Inside the boundary, convection is developed in the Si/Fe layer. A more extended but even weaker convective motion is also developed in the O/Si layer above the boundary. The outer boundary of the O/Si convection may be defined at cm, but there is only a diffuse region in this case. Surrounding the inner boundary between the Si/Fe and the O/Si layers, a thin and nearly spherical band with X(Si) 0.2 exists. As a result of the weak convection in the O/Si layer, this silicon-rich material is slowly mixed into the inner region of the O/Si layer of cm.
Model 25M (middle panels in Figure 5) develops a strong convective motion in the Si/O layer ranging from cm to cm. At the end of the simulation, outflows stream into 3 directions; the northern pole direction, 45 from the polar axis, and 135 from the polar axis, and inflows are sandwiched by the outflows. These convective flows have the convection Mach number larger than 0.1. The Si mass fraction is roughly homogenized inside the convective region, having the value of 0.3–0.4, though some fluctuations are observed especially near the outer boundary.
Model 27L (bottom panels in Figure 5) has an extended O/Si layer distributed from to cm. A large-scale convective motion is developed in this layer; a broad conical outflow with the opening angle of 45 is formed in the both polar regions, and between them, a thick inflow is formed around the equatorial plane. The convective Mach number reaches 0.12. The large scale outflow mixes the silicon rich material into the O/Si layer. The silicon mass fraction in the most part of this layer is initially 0.1, while the outflow has a higher fraction of 0.16.
Typical scale of the convection
In addition to the Mach number, the dominant angular wave number in spherical harmonics also characterize the convection. That is related to the typical size of the convective flow. This quantity is important because a large-scale convective flow can amplify the explordability of core-collapse supernovae (bernhard16_prog). The power spectrum of the radial convection velocity at is calculated as
where is the spherical harmonics function of degree and order . in the table represents value, at which has a peak.
Figure 3.2 shows the power spectrum at three different times for models 13L (top), 25M (middle), and 27L (bottom). For model 13L, has a maximum at , but the spectrum is rather flat and less energetic. The radius of the highest Mach number in the O/Si layer is cm. Although the convective mixing occurs in the inner region of cm, the convection is weaker than the outer region at the last step. Large value is explained by thin width of Si/O or O/Si layer. Model 16M forms convective mixing in the inner region of the O/Si layer similarly to model 13L. Models 18M and 21M show large (see Table 3.1) because of thin Si/O layer. In models 25M and 27L, the power spectrum peaks at and 2, respectively, and the spectrum decreases with increasing above that. The power index measured for is about in the case of model 27L, and is about for model 25 M. Models 23L and 28M show a similar power spectrum to model 27L. The power spectra of models 27M and 28L show a decreasing trend for , which is not shown in model 13L.
3.3 3D Stellar Hydrodynamics Simulation
A 3D hydrodynamics simulation is conducted using the 1D model 25M as the initial condition. The hydrodynamical evolution is followed 100 s until the central temperature reaches K, at which point the Fe core is unstable enough to collapse.
Figure 3.2 shows the time evolution of the Si mass fraction distribution of model 25M. The initial distribution of the Si mass fraction is spherically symmetric (top left panel). After the start of the simulation, the convection in the Si/O layer develops from the inner region. We see that Si-enriched plumes go up into the Si/O layer (top right panel). The convective motion reaches a steady flow by 20 s. The inhomogeneous Si mass fraction distribution introduced by the convection at 30 s is shown in the middle left panel. After a while, convection activity weakens, and the convective velocity becomes small. We will discuss the mechanism of this weakening later.
At 70 s, the Si/O layer gradually contracts, triggering the strong oxygen shell-burning at the bottom of the Si/O layer. This strong burning drives strong convection and expands the Si/O layer. We see some Si enriched plumes flow from the inner region of the Si/O layer at 75 s (see red region in the middle right panel). As a result, the strong convective flow mixes with the surroundings and increases the Si mass fraction in the whole Si/O layer. The Si mass fraction in the Si/O layer slightly increases from about 0.36 at 30 s to 0.38 at 90 s (bottom left panel). The convective motion in the Si/O layer continues until the last step of the simulation. We see inhomogeneous Si mass fraction distribution at the last step (bottom right panel).
In Figure 3.3, the time evolution of and the Si mass fraction obtained from the 3D simulation is shown. As shown by the left yellow region in the top panel, the convective motion with is obtained in the Si/O layer at 20 s. The bottom panel shows that the Si mass fraction is enhanced in this layer by that time. After a while, the Si/O layer slightly expands and the convective motion weakens. The averaged Si mass fraction in the Si/O layer does not change significantly from 20 s to 70 s.
By 70 s, the Fe core contracts and silicon shell burning as well as oxygen shell burning is enhanced. Because of this, the convective Mach number becomes large not only at 3–12 cm but also at cm. Besides, the upward enhancement of the Si mass fraction takes place in the Si/O layer during 70–80 s.
From 90 s, the convective motion becomes weak again. The whole Si/O layer gradually contracts towards the core-collapse. Meanwhile, the base temperature of the Si/O layer increases, boosting the oxygen shell burning. This results in the enhancement of the convective motion at the bottom of the Si/O layer. The Si/O layer at the end of the simulation is distributed from 3.0 cm to 10.5 cm. The maximum radial convective Mach number is , which is obtained at cm.
Figure 3.3 shows the power spectrum of the radial convective velocity at cm taken at three different times. The radius, which is located in the middle of the convective Si/O layer, is determined from the radius where the maximum is obtained in the 2D simulation. Similar to the 2D case, the power spectrum in the 3D case is calculated as
The power spectrum in the 3D simulation peaks at . This is consistent with the finding in bernhard16_prog for the star. The small maximum mode means that the convective motion is dominated by a large-scale flow. The main difference between 2D and 3D power spectra is the weaker convective motion in the 3D simulation.
It is noteworthy that at cm, at which radius the 3D shows its peak, has the peak at that is much larger than the result discussed here. The large would result from small scale convective flows that are driven by the hydrodynamically enhanced oxygen-shell burning, as in the case of the 2D simulation of models 18M and 21M.
Finally, we briefly present the result of the neutrino emission at the precollapse stage of the 3D simulation. Using the same method in Yoshida16, we calculate the time evolution of the luminosity and the spectrum of neutrino emitted via pair-neutrino process for model 25M. In Figure 3.3, the emission rate and the average energy of neutrino obtained from the 1D and the 3D simulations are compared.
The overall features of neutrino spectra for the 1D and 3D simulations are in common. The neutrino emission rate and the average temperature increases with time towards the core-collapse after 30 s in both of the simulations. The emission rates of and are larger than that of and . The average energy of is slightly smaller than that of . In the 3D simulation, however, the decrease in the emission rate and the average energy of neutrinos is observed in 70–90 s.
Furthermore, we evaluate the time evolution of the event rate of by KamLAND (e.g., Gando13), assuming that the initial mass of 25 at the distance of 200 pc is at the ongoing core collapse phase, emitting neutrinos via pair-neutrino process.
KamLAND is a one-kton size liquid-scintillation type neutrino detector (e.g., Gando13). We take the neutrino oscillation into account in a simple manner: the survival probability of is set to be 0.675 and 0.024 in the normal and inverted mass ordering, respectively (Yoshida16). The livetime-to-runtime ratio and the total detection efficiency are set to be 0.903 and 0.64 (Yoshida16). Figure 3.3 shows the evolution of the neutrino event rate by KamLAND. As an overall feature, the rate increases with time independent of mass ordering. In the 3D simulation, however, the decrease in the event rate due to the oxygen shell burning at the bottom of the Si/O layer is seen at 70–90 s. The time variation in the precollapse neutrino detection thus can be used as the indicator of the multidimensional convective activity deep inside the star, although it is practically impossible to detect such a time variation by current one kilo-ton size neutrino detectors.
Next, we examine the detectability of precollapse neutrinos by Hyper-Kamiokande. Hyper-Kamiokande is a planned Water-Cherenkov type neutrino detector, which has a huge fiducial mass of 190 kton (HKPC18). If a supernova explodes in the vicinity of the earth, high quality data of supernova neutrinos with well resolved time and energy bins are expected to be obtained (e.g., Takahashi01; Mirizzi16). However, since the threshold energy of Hyper-Kamiokande is higher than that of liquid scintillation type detector, Hyper-Kamiokande is less preferable for the observations of neutrinos produced through pair-neutrino process during the precollapse stage. Yoshida16 discussed a possibility of the neutrino observations using delayed -rays from Gd in Gd-loaded Hyper-Kamiokande (Beacom04), because the energy threshold will be reduced to 1.8 MeV, a similar energy to the case of KamLAND, in this case. Considering a moderate detection efficiency of 50%, the detection rate is about 178 times as large as the rate of KamLAND. As the same value can be applied for the present case, the time variation due to the convective motion with a time scale of seconds could be observed by Hyper-Kamiokande.
We should note that we only consider pair-neutrino process for presupernova neutrinos in this study. For a few minutes before the core-collapse, the main source of will be -decays rather than the pair-neutrino process (kato17). -decays mainly take place in the innermost Fe core. We expect that the neutrino event rate by -decays also decreases for 70–90 s because the central temperature and density decrease during this period. Although the above evaluation of the neutrino event will be underestimated, the time variation of presupernova neutrino events would give us the information about dynamics in the SiO-rich layer or an collapsing Fe core.
4 Summary and Discussions
We have performed 1D, 2D, and 3D simulations of the oxygen shell-burning just before the core-collapse of massive stars. First, we have calculated the 1D evolution of solar metallicity massive stars with ZAMS mass of 9–40 . We have considered four cases of the convective overshoot parameters in hydrogen and helium burning and the following stages.
We have searched the massive stars having O and Si-rich layers located in the range of –10 cm, because the oxygen burning shell is expected to give rise to large asphericity in the mass accretion rates onto the proto-neutron star, favoring the onset of neutrino-driven explosions. For the enclosed-mass weighted convection-width parameter , the stars with 20–30 have SiO-enriched environment in 10–10 cm. For the volume weighted convection-width parameter , the mass range indicating SiO-rich layer reduces to 13–30 . Less massive models M are also selected.
We have selected eleven models showing a large SiO-rich layer from the result of the massive star evolution and performed 2D hydrodynamics simulations of the evolution for 100 s until the central temperature reaches K. We have investigated time evolution of the Mach number of convective velocity and analyzed the power spectrum of the radial convective velocity. Based upon the analysis, we have classified the eleven models into 5 models with weak convection and 6 models with strong convection.
Weak convection models are selected to have the smaller convection Mach number than 0.1 at the end of simulations for most cases. Convection in these models shows a flat power spectra, that is, convection is dominated by small-scale flows, or a peak at a large . For three models, such weak and small-scale dominated convection takes place in the O/Si layer, which was formed through convective neon burning during previous core- and shell-silicon burning phases. Meanwhile, these models have no or only thin Si/O layer. The convection of the other two models takes place in the thin Si/O layer.
Six strong convection models show the convection Mach number larger than 0.1, which results from strong oxygen shell-burning at the base of the Si/O layer. Among them, two models have a narrow Si/O layer so that the convection is dominated by small-scale flows showing . Two models with an extended convective Si/O layer of – cm develop convection with that is dominated by large-scale flows. The remaining two models with a much more extended Si/O layer of cm also show small . The Si/O layers are too thick to be mixed homogeneously within a simulation time of s, although much more efficient convective mixing than what assumed in 1D modeling takes place in the 2D simulations. The result of the 2D simulation indicates that having the extended Si/O layer will be the condition to produce strong and large-scale convective motion.
Among the four models having a thin Si/O layer, two models are categorized into “weak” convection and the others are into “strong” convection. The main difference of these models is the time of the convection activation. For the former two models, the convection is activated at the start of the core contraction. The convection of the latter models becomes effective shortly after the start of the calculations.
Strong convection is obtained in four models having a Si/O layer with the range up to cm and two models having an extended O/Si layer. These models commonly have large values. All models with large except for 23L exhibit strong convection in the 2D simulations. Our results suggest that high density in the SiO-rich layer is conducive to producing vigorous convection just before the core-collapse and that is more suitable for utilizing convection activity than .
We have performed a 3D hydrodynamics simulation for model 25M for 100 s until the central temperature reaches K. The time evolution of the convection properties is qualitatively similar between 2D and 3D simulations. The convection is dominated by large-scale flows with either in the 3D case or in the 2D case. Main difference is weaker convective motion in the 3D simulation.
Using the 3D result of the hydrodynamics, we have evaluated the time evolution of the neutrino emission through electron-positron pair annihilation. The emission rate and the average energy of neutrino evolve similarly in 1D and 3D simulations. However, in the 3D model, the neutrino emission rate shows significant variation due to the strong oxygen shell burning at 70 s. The multi-D effect of the convective burning would be observed in presupernova neutrino events by the present and future neutrino detectors such as KamLAND and Hyper-Kamiokande.
In this study, we were only able to compute a single 3D model of the star due to the high numerical cost. However, more systematic and long-term 3D simulations employing a variety of the progenitors are apparently needed (e.g., bernhard18_prog). In order to accurately deal with neutronization of heavy elements to the iron-group nuclei, we need to not only implement bigger nuclear network but also a complete set of neutrino opacities (kato17) in our multi-D models, again albeit computationally expensive. Employing the 3D progenitor models in core-collapse supernova simulations is also mandatory to see how much the precollapse inhomogenities could foster the onset of neutrino-driven explosions (Couch15; bernhard18_prog). Impact of rotation (not to mention magnetic fields!) on the burning shells during the last stage of massive stars is yet to be studied in 3D (see chatz16 for 2D evolution models with rotation). All in all, this study is only the very first step toward the more sophisticated and systematic multi-D presupernova modeling.
This study was supported in part by the Grants-in-Aid for the Scientific Research of Japan Society for the Promotion of Science (JSPS, Nos. JP26707013, JP26870823, JP16K17668, JP17H01130, JP17K14306), the Ministry of Education, Science and Culture of Japan (MEXT, Nos. 26104007, JP15H00789, JP15H01039, JP15KK0173, JP17H05205, JP17H05206 JP17H06357, JP17H06364, JP17H06365, JP24103001 JP24103006 JP26104001, JP26104007), by the Central Research Institute of Explosive Stellar Phenomena (REISEP) of Fukuoka University and the associated research projects (Nos.171042,177103), and by JICFuS as a priority issue to be tackled by using Post ‘K’ Computer. K. T. was supported by Japan Society for the Promotion of Science (JSPS) Overseas Research Fellowships.
We show some evolution properties of the massive star models. As shown in Section 2, we considered four groups of models with different overshoot parameters. The overshoot parameter during the hydrogen and helium burning affects the evolution on the Hertzsprung-Russell (HR) diagram. We compare the HR diagram of models 20L and 20M with that of the 20 star models shown in Martins13. Figure 4 shows the HR diagram of models 20L, 20M, 20 star models of Stern (Brott11), GENEC (Ekstroem12), MESA (Martins13), and Starevol (Martins13). The main difference in the HR diagram between model 20L and model 20M is the main-sequence (MS) band width, i.e., the difference of the effective temperature between ZAMS and the hydrogen burning termination. The MS band width of model 20L is almost identical to Stern model. On the other hand, model 20M is almost identical to GENEC, MESA, and Starevol. As explained in section 2, the overshoot parameters of models L and M are determined based on the calibrations to early-B type stars in the LMC similar to Stern model and the MS width observed for AB stars in open-clusters of the MW Galaxy similar to GENEC model, respectively. Except for the MS width, we do not see any large differences in the HR diagram among these models. The luminosity during the change from blue supergiant to red supergiant is almost constant in the range of –5.1 for all models. The effective temperature of all models in the red supergiant stage is commonly .
We show the HR diagram of ten different ZAMS mass models for models and in Figure 4. Stars with for model and with for model end as a red supergiant. Stars with heavier evolve to a Wolf-Rayet star.
We show the evolution of the relation between the central density and the central temperautre in Figure 4. All stars except for the ZAMS mass of 9 in models M and M end as the core-collapse of an Fe core. The 9 star in models M does not bring about Ne ignition (panel (): red line). We also confirmed no Ne ignition in the 9 star in model M until the central density becomes . These star are expected to end as an ONe white dwarf or an electron-capture SN. The off-center Ne ignition occurs for the 9 star for models L (panel (): red line) and the 11 star for models M (panel (): green line). The similar off-center neon burning is also seen in a 10 star of models M and M. The low-mass M and M models ignite silicon at an off-center region.
Whether the carbon core burning occurs convectively or radiatively depends on the stellar mass. The convective carbon-core burning occurs in the stars with for models L and L and with for models M and M. We suggest that the maximum CO-core mass for the convective carbon burning is between 4.6 and 4.9 among our models.
We calculate the massive star evolution taking account of a weak overshoot after the helium burning for models L and M. Despite the small value compared with the hydrogen and helium burning, the overshoot in the advanced stage complicatedly affects advanced evolution. Here we show Kippenhahn diagram of models 25M and 25M in Figure 4 for comparison. In these models, the carbon core burning occurs radiatively and the first carbon shell burning (C(I)) ignites at . In model 25M the convective region of the C(I) burning extends both inward and outward with the help of overshoot. The convective C layer extends in the range 1.0–2.3 . Then, the second carbon shell burning (C(II)) ignites at and the convective region extends again inward and outward. The inner convection boundary of the second carbon shell burning reaches 1.74 . This inner boundary restricts the region of the following burning. The Si core is formed through the O core burning (O(c)) and the first (O(I)) and second (O(II)) oxygen shell burning. The third oxygen shell burning (O(III)) extends the Si layer up to 1.67 but the boundary is still below the inner boundary of the second carbon shell burning. The Si core (Si(c)) and the following four silicon shell burning (Si(I)-(IV)) form an Fe core of 1.56 . The left panel of Figure 4 shows the mass fraction distribution of model 25M. There is a SiO-rich layer betwen the Si layer and the O/Ne layer but it is very thin (the width is about cm).
Model 25M indicates the evolution different from model 25M from the first carbon shell burning (C(I)). The first carbon shell burning ignites at 1.3 and the convective region extends outward to 2.9 . The convection does not extend inward in this model. Then, the second carbon shell burning (C(II)) starts at 2.9 at the coordinate. Since the second carbon burning shell has a large radius, the O core burning (O(c)) the first oxygen shell burning (O(I)) extend more effectively outward and form a larger Si core. The Si core (Si(c)) and the following three silicon shell burning (Si(I)-(III)) form an Fe core of 1.68 . The second oxygen shell burning (O(II)) makes a large Si/O layer between the top of first oxygen shell burning (1.9 ) and the bottom of the second carbon shell burning (2.9 ). The mass fraction distribution of model 25M is shown in the right panel of Figure 4. We see that the large Si/O layer in this figure.
The compactness parameter, , is considered as a quantity that correlates to the dynamics during the gravitational collapse (Oconnor11). It is defined as an inverse of the radius, inside which the mass of 2.5 is enclosed,
Several works have reported that the small compactness parameter is favorable for the supernova explosions (Nakamura15; Pejha15; Sukhbold16; Horiuchi14), although the criterion has not yet converged (Ugliano12; Suwa16; Ertl16; Burrows18). This parameter is also important to predict the neutrino emission (Horiuchi17). The tight correlation of to the CO core mass has been shown (Limongi18), and detailed study on the compactness is performed by Sukhbold14; Sukhbold18.
Figure 4 shows the compactness parameter of our models. Although the results show a large scatter, roughly increases with ZAMS mass. The scatter is reduced when the x-axis is changed from the ZAMS mass to the He-core mass (panel ()) as shown in Sukhbold18. As for the correlation between and the He core mass, the scatter is small for the models with . Besides, we do not see clear dependence on in this mass range. For more massive models, the scatter is larger, and the different overshoot parameters also give an influence on in a complicated manner.
Because of the limited number of our models, it is difficult to make a detailed comparison between our results on the compactness parameter (Figure 4) with Sukhbold18 (see their Figure 8). We do see a diversity of in , however, the range of for is almost within the same range as in Sukhbold18. Although a more systematic study needs to be done as in Sukhbold16; Sukhbold18, the overall feature (e.g., the increasing trend of the compactness parameter with the ZAMS masses) is roughly in accordance with Sukhbold18.