One-, Two-, and Three-dimensional Simulations of Oxygen Shell Burning Just Before the Core-Collapse of Massive Stars

TL;DR

This study uses 2D and 3D hydrodynamics simulations to analyze convective oxygen shell burning in massive stars before core-collapse, revealing how shell structure influences turbulence and potential supernova explosion conditions.

Contribution

It provides the first detailed 3D simulation of oxygen shell convection in massive stars, linking shell structure to turbulence and explosion likelihood.

Findings

Thick Si/O layers promote large-scale convection and turbulence.

3D models show lower turbulence velocities than 2D models.

Neutrino emission modulation could reveal presupernova structural changes.

Abstract

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 an initial mass of 9-40 $M_\odot$. Four different overshoot parameters are applied, and CO core mass trend similar to previous works is obtained in the 1D models. Selecting eleven 1D models that have a silicon and oxygen coexisting layer, we perform 2D hydrodynamics simulations of the evolution $\sim$100 s until the onset of core-collapse. We find that convection with large-scale eddies and the turbulent Mach number $\sim$0.1 is obtained in the models having a Si/O layer with a scale of 10$^8$ cm, whereas most models that have an extended O/Si layer up to a few $\times 10^9$ cm exhibit lower turbulent velocity. 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 $M_\odot$ model, which exhibits large-scale convection in the 2D models. The 3D model develops large ($\ell = 2$) convection similar to the 2D model, however, the turbulent velocity is lower. 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.

Figures

Tables

Table 1. Table 1: Stellar evolution model sets and the corresponding overshoot parameters. f ov subscript 𝑓 ov f_{\rm ov} is the overshoot parameter during the H and He burning phases. f ov , A subscript 𝑓 ov A f_{\rm ov,A} is the one during the more advanced stages, namely after the He burning phase. For the meanings of the notations “L” and “M”, see the text. Convective overshoot during the advanced stages is considered for model sets with the subscript “ A ”.

Model	$f_{\mathrm{ov}}$	$f_{\mathrm{ov},A}$
L	0.03	0
LA	0.03	0.002
M	0.01	0
MA	0.01	0.002

Table 2. Table 2: 2D Model Properties and SiO-coexistence parameters. SiO-coexistence parameters f M , SiO subscript 𝑓 𝑀 SiO f_{M,{\rm SiO}} and f V , SiO subscript 𝑓 𝑉 SiO f_{V,{\rm SiO}} are obtained from the result of 1D evolution simulations. ⟨ M ​ a 2 ⟩ max 1 / 2 subscript superscript delimited-⟨⟩ 𝑀 superscript 𝑎 2 1 2 max \langle Ma^{2}\rangle^{1/2}_{\rm max} represents the maximum convective Mach number obtained at a radius of r 𝑟 r ( ⟨ M ​ a 2 ⟩ max 1 / 2 subscript superscript delimited-⟨⟩ 𝑀 superscript 𝑎 2 1 2 max \langle Ma^{2}\rangle^{1/2}_{\rm max} ) at the end of the 2D simulations. “Layer” represents the composition of the convective region. ℓ max subscript ℓ max \ell_{\rm max} represents the ℓ ℓ \ell value at which c ℓ 2 superscript subscript 𝑐 ℓ 2 c_{\ell}^{2} has a peak (see equation ( 9 )). d c / H P subscript 𝑑 𝑐 subscript 𝐻 𝑃 d_{c}/H_{P} represents the width of convective region normalized by the local scale height. These quantities are all estimated at the last step of the simulations. See the text for more detailed definition.

Model	$f_{M,\mathrm{SiO}}$	$f_{V,\mathrm{SiO}}$	${⟨M{a}^2⟩}_{\max }^{1/2}$	$r$(${⟨M{a}^2⟩}_{\max }^{1/2}$)	Layer	${\ell }_{\max }$	$d_c/H_P$
				(108 cm)
Low-$Ma$
13LA	0.27–0.28	0.95–0.96	0.018	11.6	O/Si	12	6.22
16MA	0.24–0.24	0.90–0.91	0.015	3.9	O/Si	4	3.20
18MA	0.57–0.58	0.91–0.91	0.131	3.1	Si/O	14	1.06
21MA	0.47–0.47	0.91–0.95	0.134	3.0	Si/O	8	4.42
23LA	0.75–0.80	0.78–0.80	0.069	11.5	O/Si	4	5.20
High-$Ma$
22L	0.57–0.61	0.77–0.82	0.108	9.4	Si/O	2	2.50
25M	0.75–0.79	0.91–0.94	0.160	5.8	Si/O	3	3.65
27LA	0.59–0.66	0.76–0.76	0.179	45.0	O/Si	2	4.56
27M	0.58–0.65	0.37–0.40	0.134	4.7	Si/O	10	2.44
28LA	0.60–0.68	0.37–0.42	0.117	5.3	Si/O	8	1.81
28M	0.83–0.90	0.90–0.95	0.369	14.6	O/Si	2	4.08

Equations

$$D_{\mathrm{cv}}^{\mathrm{ov}}=D_{\mathrm{cv},0}\exp \left(-2\frac{\Delta r}{f_{\mathrm{ov}}{H}_{P0}}\right),$$

$${\partial }_t\rho +\nabla \cdot \left(\rho \mathbf{v}\right)$$

$${\partial }_t\left(\rho \mathbf{v}\right)+\nabla \cdot \left(\rho \mathbf{v}\mathbf{v}+p\mathbf{\delta }\right)$$

$${\partial }_t{e}_{\mathrm{tot}}+\nabla \cdot \left[(e_{\mathrm{tot}}+p)\mathbf{v}\right]$$

$${\partial }_t\rho X_i+\nabla \cdot \left(\rho X_i\mathbf{v}\right)$$

$$f_{M,\mathrm{SiO}}$$

$$f_{V,\mathrm{SiO}}$$

$$⟨M{a}^2{⟩}^{1/2}(r)={\left[\frac{\int \rho \{(v_r-⟨v_r⟩{)}^2+v_{\theta }^2+v_{\varphi }^2\}d\Omega }{\int \rho c_s^{2}d\Omega }\right]}^{1/2},$$

$$c_{\ell }^2={\left|\int (v_r-⟨v_r⟩)Y_{\ell 0}^{\ast }(\theta )d\Omega \right|}^2,$$

$$c_{\ell }^2=\sum _{m=-\ell }^{\ell }{\left|\int (v_r-⟨v_r⟩)Y_{\ell m}^{\ast }(\theta ,\varphi )d\Omega \right|}^2.$$

$${\xi }_{2.5}=\frac{2.5}{R(M_r=2.5M_{\odot })/1000\mathrm{km}}$$

Full text

stars: massive – supernovae:general – convection – hydrodynamics

††journal: ApJ††software: HOSHI (Takahashi et al., 2016, 2018), 3DnSNe (Takiwaki et al., 2016; Nakamura et al., 2016; Kotake et al., 2018)

1 Introduction

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 et al. (2015); Janka et al. (2016); Patat (2017) for reviews). To facilitate the neutrino-driven mechanism of CCSNe (Bethe & Wilson, 1985), multi-D hydrodynamics instabilities such as neutrino-driven convection and the standing accretion shock instability (Blondin et al., 2003) 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., Vartanyan et al. (2019); O’Connor & Couch (2018); Müller et al. (2017); Roberts et al. (2016); Nakamura et al. (2016); Melson et al. (2015b); Summa et al. (2016); Lentz et al. (2015); Takiwaki et al. (2014); Hanke et al. (2013) 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 typically observed value (e.g., $\sim$1051erg, Tanaka et al., 2009). 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 (2017) for a review), general relativity (GR, e.g., Müller et al. (2012); Ott et al. (2013); Kuroda et al. (2012, 2016)), rapid rotation (e.g., Marek & Janka (2009); Suwa et al. (2010); Takiwaki et al. (2016); Summa et al. (2018); Harada et al. (2018)) and/or magnetic fields (e.g., Obergaulinger et al. (2006); Mösta et al. (2014); Guilet & Müller (2015); Masada et al. (2015); Obergaulinger & Aloy (2017)), and sophistication in the neutrino opacities (Melson et al., 2015a; Bollig et al., 2017; Burrows et al., 2018; Kotake et al., 2018) and in the neutrino transport schemes (e.g., Sumiyoshi & Yamada (2012); Richers et al. (2017); Nagakura et al. (2018); Just et al. (2018)). In this work, we focus on the first item listed in the above list.

Couch & Ott (2013) were the first to demonstrate that the inhomogeneities seeded by convective shell burning fosters the onset of a neutrino-driven explosion (see also Fernández et al. (2014); Couch & Ott (2015); Müller & Janka (2015); Burrows et al. (2018)). This is because the infalling perturbation that could be amplified in the supersonic accretion (Takahashi & Yamada, 2014; Nagakura et al., 2013, 2019) enhances turbulence behind the postshock material leading to the reduction of the critical neutrino luminosity for shock revival (e.g., Müller & Janka (2015); Abdikamalov et al. (2016)). In these studies, the non-spherical structures in the burning shells, although physically motivated, were treated in a parametric manner due to the paucity of the multi-D stellar evolution models covering the lifespan of massive stars up to the iron core-collapse. Currently one-dimensional (1D) stellar evolution calculations are the only way to accomplish this (Woosley et al., 2002; Woosley & Heger, 2007; Sukhbold et al., 2018), where the errors introduced from the omission of multi-D effects are absorbed into the free parameters of MLT, namely the mixing length theory, (e.g., Kippenhahn et al. (2012)).

The truly multi-D hydrodynamics stellar evolution calculations have been done over several turnover timescales of convection (limited by the affordable computational resources) in selected burning shells (e.g., Meakin & Arnett (2007); Viallet et al. (2013); Campbell et al. (2016); Cristini et al. (2017); Cristini et al. (2019) for different burning shells, and see Arnett & Meakin (2016) 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 (Arnett, 1994; Bazan & Arnett, 1994; Bazán & Arnett, 1998; Asida & Arnett, 2000; Kuhlen et al., 2003; Meakin & Arnett, 2006, 2007; Arnett & Meakin, 2011; Chatzopoulos et al., 2014, 2016; Jones et al., 2017).

More recently, ground-breaking attempts to evolve convective shells in 3D prior to the onset of collapse have been first reported by Couch et al. (2015) for silicon shell burning in a 15 $M_{\odot}$ star and by Müller et al. (2016) for oxygen shell burning in an 18 $M_{\odot}$ star (and also in 11.8, 12, and 12.5 $M_{\odot}$ stars by Müller et al. (2018)). Couch et al. (2015) obtained earlier onset of a neutrino-driven explosion for the 3D progenitor model of the $15M_{\odot}$ star compared to that in the corresponding 1D progenitor model. By performing 3D GR simulations with more advanced neutrino transport scheme, Müller et al. (2017) obtained a neutrino-driven explosion with the seed perturbations. In comparison, this 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 for precollapse seed perturbation growth, Collins et al. (2018) 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 $\sim$16 and 26 $M_{\odot}$ are most likely to exhibit large-scale convective overturn with high convective Mach numbers, leading to the most favorable condition for perturbation-aided explosions. In fact the 3D progenitor model of the 18 $M_{\odot}$ star (Müller et al., 2016) is 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 zero-age main sequence (ZAMS) masses between 9 and 40 $M_{\odot}$ with the HOSHI code developed by Takahashi et al. (2016, 2018). Based on the 1D results, we select $\sim$ ten 1D progenitors that have extended and enriched O and Si layers, presumably leading to vigorous convection. At a time of $\sim 100$ s before the onset of collapse, the 1D evolution models are mapped to multi-D hydrodynamics code (a branch of 3DnSNe, e.g., Takiwaki et al. (2016); Nakamura et al. (2016); Kotake et al. (2018)). We perform axisymmetric (2D) simulations for the selected progenitors having an extended O and Si-rich layer and investigate the features of their convective motion, especially the convective-eddy scale and the turbulent Mach number. 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 features 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

In this section, we briefly summarize the numerical setups of our stellar evolution calculations in 1D (Section 3.1), 2D (Section 3.2), and 3D (Section 3.3).

2.1 1D Stellar Evolution