# Does Covid-19 Mass Testing Work? The Importance of Accounting for the
Epidemic Dynamics

Davide Ferrari  Faculty of Economics and Management, Free University of
Bolzano  Steven Stillman  Faculty of Economics and Management, Free University
of Bolzano  Mirco Tonin  Faculty of Economics and Management, Free University
of Bolzano

###### Abstract

Mass antigen testing has been proposed as a possible cost-effective tool to
contain the Covid-19 pandemic. We test the impact of a voluntary mass testing
campaign implemented in the Italian region of South Tyrol on the spread of the
virus in the following months. We do so by using an innovative empirical
approach which embeds a semi-parametric growth model - where Covid-19
transmission dynamics are allowed to vary across regions and to be impacted by
the implementation of the mass testing campaign - into a synthetic control
framework which creates an appropriate control group of other Italian regions.
We find that the mass test campaign decreased the growth rate of Covid-19 by
39% which corresponds to a reduction in the total additional cases of 18%, 30%
and 56% within ten, twenty and forty days from the intervention date,
respectively. Our results suggest that mass testing campaigns are useful
instruments for mitigating the pandemic.

Keywords: Mass testing; Covid-19; Growth model; Semi-parametric growth model;
Synthetic controls.

##  1 Introduction

A wide variety of interventions have been used in an attempt to stop the
spread of Covid-19 across the globe (i.e. lockdown, quarantines, business
closures, mobility restrictions, school closures) (Haug et al.,, 2020; Tian et
al.,, 2020). Many of these are also the main candidate policies for stopping
the spread of other large-scale epidemic outbreaks, such as Ebola. One
particularly low-cost intervention that has been tried in a few locations is
mass testing of a population in order to identify asymptomatic carriers
(Holt,, 2020; Atkeson et al.,, 2021). Theoretical work suggests that mass
testing could reduce daily infections by up to 30 percent (Bosetti et al.,,
2020), while Pavelka et al., (2021) evaluates the impact of a mass testing
campaign undertaken in Slovakia in late 2020 and finds that it temporarily
decreased the growth of Covid-19 by 70 percent. Importantly, many countries
worldwide offer open public testing indicating that frequent mass testing is a
feasible intervention (Hasell et al.,, 2020).

In this paper, we test the impact of a voluntary mass testing campaign
implemented in the Italian region of South Tyrol on the spread of the virus in
the following months. We do so by using an innovative empirical approach which
embeds a semi-parametric growth model into a synthetic control framework.
Specifically, we first use the synthetic control approach to create a control
group which is a weighted-average of other Italian regions that best follow
the dynamics of Covid-19 transmission in South Tyrol prior to the mass testing
campaign (Abadie et al., 2010a, ; Abadie and Gardeazabal, 2003a, ).111We use
this approach because it uses the data to derive the optimal control region
for the South Tyrol. Mitze et al., (2020) and Cho, (2020) use this method to
evaluate the impact of masks and other nonpharmaceutical interventions on the
spread of Covid-19, however these papers implement it in a traditional
framework which does not account for the underlying growth dynamics of an
epidemic like Covid-19. We then estimate on appropriately weighted data a
semi-parametric growth model where Covid-19 transmission dynamics (i.e. growth
rates) are allowed to vary across regions and to be impacted by the
implementation of the mass testing campaign. Importantly, this approach is in
the spirit of difference-in-difference models which compare changes in
outcomes over time in a treated location to changes in outcomes over time in
otherwise similar control locations and hence can isolate the impact of the
mass testing campaign from national level policies regarding freedom of
movement, business and school closures, hygienic measures, etc. since no other
policy changed at the same time in South Tyrol.

This approach has a number of benefits over what has so far been done in the
literature that attempts to estimate the impact of public health interventions
on the spread of Covid-19, as well as on other contagious diseases.
Specifically, previous papers use either full parametric models (Wu et al.,,
2020; Aviv-Sharon and Aharoni,, 2020; Chénangnon et al.,, 2020; Pelinovsky et
al.,, 2020) or traditional synthetic control or difference-in-difference
models (Hsiang et al.,, 2020; Tian et al.,, 2020; Mangrum and Niekamp,, 2020;
Alexander et al.,, 2020; Dave et al.,, 2021; Singh et al.,, 2021; Mitze et
al.,, 2020; Cho,, 2020). While fully parametric models based on biological
growth theory have proved very useful for assessing the dynamics of disease
outbreaks, their overly rigid a-priori parametric assumptions do not allow
them to fit the typical variation in short- and medium-run dynamics seen in
the spread of Covid-19. This is a potential concern with the findings of
Pavelka et al., (2021), especially since the testing campaign in Slovakia was
a national level campaign and other policies were being implemented at the
same time to reduce transmission rates. On the other hand, traditional
difference-in-difference and synthetic control models, by assuming that
treatment and control regions would follow parallel trends in a counterfactual
reality without any intervention, are potentially biased because growth
dynamics in different regions will depend on the prior contagion rate.

Although there is a rich literature on estimating model-free growth curves in
other fields, the use of non-parametric or semi-parametric approaches has not
been sufficiently explored in the context of health policy
analysis.222Exceptions include Chénangnon et al., (2020) which fits a flexible
growth model curve to reported positive Covid-19 cases and Lee et al., (2020)
which specifies dependence of the parameters in the Richards growth model on
covariates to assess country-level risk factors due to Covid-19 within a
hierarchical Bayesian framework. However, neither of these papers estimate the
impact of a particular policy intervention. We show that the estimated impact
of the mass test campaign in South Tyrol is sensitive to how one models
Covid-19 transmission dynamics and flexible models better fit the data and
generate more robust estimates of the impact of the mass testing campaign.
Furthermore, our approach has the added benefit that it generates estimates of
the disease proliferation rate both with and without intervention, thus
gaining insight on the transmission pathways.

Overall, we find that the mass test campaign in the South Tyrol decreased the
growth rate of Covid-19 by 39% (95% confidence interval: 29-49%). This
corresponds to a reduction in the total additional cases of 14%, 18%, 30% and
56% within seven, ten, twenty and forty days from the intervention date,
respectively. Importantly, this large impact was achieved even though the
campaign was entirely voluntary with no incentives to participate.333In
Slovakia, individuals who did not participate were required to quarantine for
ten days, hence it was not entirely a voluntary test. Our results are in line
with the predictions made in Bosetti et al., (2020) based on an
epidemiological model. We find a smaller impact than Pavelka et al., (2021)
reports for the mass testing campaign in Slovakia; however, the Slovakian
intervention featured a multiple round campaign in addition to concurrent
interventions, including a one-week lockdown.444In line with our findings, a
recent working paper by Kahanec et al., (2021) uses a different approach to
analyze the impact of the Slovakia campaign and finds that infections were
reduced by about 25-30% in the two weeks after individuals in some districts
were tested a second time.

##  2 The Mass Test

The population of South Tyrol was invited to take part in a mass testing
campaign in late November 2020 using rapid antigen tests, which involve a
nasal and throat swab. Authorities set up around 300 testing centers, where
professional health care workers carried out the tests, with the support of
volunteers from the civil protection agency, the voluntary fire services and
other organizations for handling the logistics and the administration. All
residents were invited to participate, with the exception of children below
the age of five, people with Covid-19 symptoms, those on sick leave, those who
had tested positive and isolated in the last three months, and those who had
recently tested positive or were in quarantine or self-isolating. People with
a prior appointment for a PCR test, those regularly tested for work reasons,
and individuals in social care were also not tested.

Testing centers generally operated from 8am to 6pm from Friday, 20 November to
Sunday, 22 November. During this period, people could show up at any of the
centers throughout the region. In some municipalities, it was possible to
register online and some published suggested centers and time slots based on
the address of residence. It was also possible to be tested at some pharmacies
and GPs in the period 18 to 25 November. People only needed a valid ID and a
European Health Insurance card. They filled in a form with an email address,
where they would receive, generally within a day, an encrypted file with the
outcome, and a mobile number, where they would receive an SMS with the code to
open the file. In case of a negative result, people were advised to continue
following prevention measures like social distancing and mask wearing. In case
of a positive result, people had to isolate for 10 days if asymptomatic and
contact their doctor if they developed symptoms.

Participation in the mass testing was voluntary and encouraged by a massive
communication campaign, providing information (with material available also in
Albanian, Arabic, English, French and Urdu, as well as in simple language for
kids), as well as endorsements by public figures. The goal was to identify
asymptomatic cases in the population and hence reduce virus transmission. The
headline of the campaign was “Together against coronavirus”, using appeals
like “Let’s break the infection wave together and pave the way towards a
gradual return to normality!”. In the end, 72 percent of eligible residents
volunteered to be tested (Stillman and Tonin,, 2021). In comparison, 83
percent of the eligible population in Slovakia decided to be tested when the
alternative was to quarantine for 10 days.

##  3 Methods

In this section, we first discuss how we create the optimal control group for
evaluating the mass testing campaign in the South Tyrol. We then discuss how
we implement the semi-parametric growth model, how it compares to parametric
models and how we imbed a difference-in-difference type framework in this
model. Finally, we discuss statistical inference.

###  3.1 Defining the control group

As there is no obvious a priori control group, in order to evaluate the effect
of the screening intervention in South Tyrol, we constructed an optimal
control group (Synthetic South Tyrol) using the synthetic control methodology
of Abadie et al., 2010b . This method constructs a weighted combination of
data from the control group to approximate the behavior of the intervention
group in terms of pre-intervention characteristics.

Suppose we observe new cases $Y_{it}$ in regions $i=1,\dots,N+1$ at times
$t=1,\dots,T$ and assume intervention occurs at time $T_{0}+1$ so that
$1,\dots,T_{0}$ are pre-intervention periods. Without loss of generality the
first region ($i=1$) represents South Tyrol, which is exposed to the
intervention, so we have $N$ remaining donor regions contributing to the
synthetic control (in our case, the other 20 Italian regions or autonomous
provinces in the donor pool). Let $\mathbf{w}=(w_{1},\dots,w_{N})^{\top}$ be a
$(N\times 1)$ vector of non-negative weights such that $\sum_{i=1}^{N}w_{i}=1$
and $w_{i}\geq 0$ for all $i=1,\dots,N$, where $w_{i}$ denotes the weight of
region $i$ in the synthetic South Tyrol. Let
$\mathbf{u}_{i}=(x_{i1},\dots,x_{iT_{0}},z_{i1},\dots,z_{iT_{0}})^{\top}$ be a
$(2T_{0}\times 1)$ vector containing pre-treatment values for cumulative
number of cases ($x_{it}$) and number of tests ($z_{it}$) for region $i$.
Thus, the vector $\mathbf{u}_{1}$ contains pre-treatment values for South
Tyrol and the matrix $(2T_{0}\times N)$ matrix
$\mathbf{U}_{0}=(\mathbf{u}_{2},\dots,\mathbf{u}_{N+1})$ with columns
$\mathbf{u}_{2},\dots,\mathbf{u}_{N+1}$ contains pretreatment values for all
the remaining donor regions.

The synthetic control method selects the vector of weights
$\widehat{}\mathbf{w}$ that determines the best control region by solving the
following quadratic optimization problem:

$\text{min}_{\mathbf{w}}\ \
(\mathbf{u}_{1}-\mathbf{U}_{0}\mathbf{w})^{\top}\mathbf{V}(\mathbf{u}_{1}-\mathbf{U}_{0}\mathbf{w}),\
\ \text{subject to }\sum_{i=1}^{N}w_{i}=1,$ (1)

where $\mathbf{V}$ is a $2T_{0}\times 2T_{0}$ symmetric and positive
semidefinite matrix to allow different weights to the variables in
$\mathbf{U}_{0}$ depending on their predictive power on the outcome. To stress
dependence on $\mathbf{V}$, we use the notation
$\widehat{}\mathbf{w}=\widehat{}\mathbf{w}(\mathbf{V})$ for the solution to
(1). To select $\mathbf{V}$, we minimize the mean squared prediction error
(MSPE) of the outcome variable over pre-intervention periods as described
inAbadie and Gardeazabal, 2003b  and Abadie et al., 2010b . Specifically, if
$\mathbf{Y}_{1}$ is the $(T_{0}\times 1)$ vector with the values of the
outcome variable (new cases) for South Tyrol and $\mathbf{Y}_{0}$ is the
$(T_{0}\times N)$ analogous matrix for the control units, then $\mathbf{V}$ is
selected by minimizing the
$\text{MSPE}(\mathbf{V})=\|\mathbf{Y}_{1}-\mathbf{Y}_{0}\widehat{}\mathbf{w}(\mathbf{V})\|^{2}$
over all positive definite and diagonal matrices $\mathbf{V}$.

The above procedure is implemented in the function synth of the R package
Synth which uses the Nelder-Mead and BFGS algorithms as default options, and
then picks the solution with the lowest MSPE (Abadie et al.,, 2011). We then
use the weights generated by this procedure when estimating the model
described in the next section.

###  3.2 Semi-parametric growth model

To analyze the dynamics of the Covid-19 epidemic over time, we develop a
flexible semi-parametric growth curve approach. Let $x_{t}$ denote the
cumulative size of the detected infected population at time $t$. In classic
parametric growth curve analysis, the dynamics of $x_{t}$ is represented by
the derivative of $x_{t}$, which is often assumed to have the form

$\dfrac{\partial x_{t}}{\partial t}=\rho\times x^{p}\times g(x_{t}),$ (2)

where $\rho$ is the intrinsic growth rate determining the time scale of the
epidemic process, $p\in[0,1]$ is the “deceleration of growth" parameter
capturing different growth profiles and $g$ is a smooth non-increasing
function of $x_{t}$ possibly depending on other model parameters.

A number of well-known growth models can be recovered from Equation (2)
depending on the choice of $g$; although the literature is too vast to be
covered here, we provide a few examples. A basic model is the exponential
growth (EG) model corresponding to $g(x_{t})=1$. The EG model assumes that an
epidemic continues to grow following the same process as in the past with the
growth path completely specified by $p$: constant incidence ($p=0$), sub-
exponential growth ($0<p<1$) and exponential growth ($p=1$); see, e.g. Wu et
al., (2020). Although this may be useful for representing the early stages of
an epidemic, it is an upper bound scenario since outbreaks often slow down and
reach saturation capacity after initial exponential growth. One popular
extension is the generalized logistic growth (GLG) model
$g(x_{t})=x_{t}^{p}(1-x_{t}/k)$, where $k$ represents the total size of the
epidemic, i.e. the asymptotic number of infections over the whole epidemics
(Wu et al.,, 2020). A slightly more flexible model is the generalized Richards
growth (GRG) model $g(x_{t})=x_{t}^{p}[1-(x_{t}/k)^{a}]$, where the parameter
$a$ allows for deviations from the S-shaped dynamics of the classical GLG (Wu
et al.,, 2020; Chowell,, 2017).

Despite the wide range of parametric growth models available, none of these
rigid apriori specifications do a good job at fitting the rich dynamics
observed in real Covid-19 data, which potentially leads to biased inferences.
Instead, we estimate a more flexible model, referred to as a semi-parametric
growth (SPG) model, that includes additional flexibility in two dimensions.

We first add flexibility by estimating the function $g$ in the following semi-
parametric specification:

$g(x_{t})=\exp\left\\{h(x_{t})\right\\}=\exp\left\\{\sum_{j=1}^{q}\eta_{j}b_{j}(x_{t})\right\\},$
(3)

where the $\eta_{i}$s are unknown coefficients and the $b_{j}(x_{t})$ are
given basis functions, such basis spline functions. The basis expansion
$h(x_{t})=\sum_{i=1}^{q}\eta_{i}b_{j}(x_{t})\\}$ allows us to approximate an
arbitrary smooth function, provided that $q$ is large enough. Overall, this
specification requires fewer assumptions about the data and fits the data
better in situations where the true outbreak trajectory is hard to specify in
advance, as is the case for the Italian Covid-19 data.

We also allow the growth rate parameter $\rho$ to depend on other explanatory
variables. Since our main inferential interest is to assess the impact of the
policy treatment on pre- and post-treatment growth differences, we assume the
growth rate follows a log-linear difference-in-differences type of
specification:

$\log\\{\rho_{it}(\theta)\\}=\alpha+\beta\ \texttt{Int}_{t}+\gamma\
\texttt{Reg}_{i}+\delta\
\texttt{Int}_{t}\times\texttt{Reg}_{i}+\boldsymbol{\xi}^{\top}\mathbf{z}_{it},$
(4)

where $\texttt{Int}_{t}$ is a dummy variable for the policy intervention
taking value 0 up to the intervention date and 1 afterwards,
$\texttt{Reg}_{i}$ is a dummy variable taking value 1 for the geographical
region exposed to treatment (South Tyrol) and 0 otherwise, and
$\mathbf{z}_{it}$ is a $q\times 1$ vector of additional controls (such as the
difference in the number of diagnostic tests compared to the previous day,
weekly seasonality effects, etc). The overall vector describing the growth
parameters is denoted by
$\theta=(\alpha,\beta,\gamma,\delta,\boldsymbol{\xi}^{\top})^{\top}$.

Our main inferential interest is in the coefficient $\delta$ for the
interaction between the intervention and the geographical area. This parameter
measures the effect of the policy intervention and is calculated in a similar
way as in a traditional difference-in-differences model by comparing the
change over time in the outcome variable for the South Tyrol compared to the
change over time for the control region (Synthetic South Tyrol). We also
assess the impact of policy intervention by examining its impact on the
transmission growth rate. From equation (4), the relative change in the growth
rate can be computed as
$\Delta\rho_{i}(\theta)=[\rho_{i1}(\theta)-\rho_{i0}(\theta)]/\rho_{i0}(\theta)=\exp\left\\{\beta+\delta\
\texttt{Reg}_{i}\right\\}-1$. Therefore, the relative change in the treatment
and control groups are $\Delta\rho_{1}=\exp\\{\beta+\delta\\}-1$ and
$\Delta\rho_{0}=\exp\\{\beta\\}-1$, respectively.

###  3.3 Statistical inference

Let $\\{Y_{it}\\}$ denote new COVID-19 cases observed in region $i$ at time
$t$. We follow Chénangnon et al., (2020) and estimate the model parameters by
assigning to each $Y_{it}$ an appropriate statistical distribution with mean
$E[Y_{it}]=\partial x_{it}/\partial t=\rho_{it}(\theta)x_{it}^{p}g(x_{it})$,
where $x_{it}$ is the cumulative number of cases observed in region $i$ at
time $t$ and $\rho_{it}(\theta)$ is a region- and time- specific growth rate
parameter. For $Y_{it}$ we consider both the Poisson and negative binomial
(NB) distributions, which are appropriate models for count data. Using a log-
link function to relate the expected new cases to other predictors, the
working statistical model for inference can be written as

$\displaystyle\log\left\\{E(Y_{it})\right\\}$
$\displaystyle=\log\left\\{\partial x_{it}/\partial
t\right\\}=\log\\{\rho_{it}(\theta)\\}+p\log(x_{it})+h(x_{it}),$ (5)

$i=0,1$ and $t=1,\dots,T$, where $\rho_{it}(\theta)$ is the time- and region-
specific growth rate defined in Equation 4.

Equation 5 represents a generalized additive model (GAM); e.g., see Hastie and
Tibshirani, (1990); Wood, (2017) for an introduction on GAMs. Estimates for
the growth rate parameters
$\widehat{\theta}=(\widehat{\alpha},\widehat{\beta},\widehat{\gamma},\widehat{\delta},{\widehat{}\boldsymbol{\xi}}^{\top})^{\top}$,
growth deceleration $\widehat{p}$ and for the smooth function
$\widehat{h}(\cdot)$ are obtained by running a penalized likelihood estimator
with penalty depending on a smoothing parameter for the non-parametric part of
the model. We use thin-plate radial basis spline functions for the terms
$b_{i}(x_{t})$ and smoothing parameter tuned via maximum likelihood
estimation.

Estimates are obtained using the function gam in the R package mgcv (Wood and
Wood,, 2015). See Wood, (2017) and references therein for details on the
estimation procedure and implementation. Approximate variances and covariances
for the estimated parameters are extracted using the function vcov.gam in the
R package mgcv. An estimate of $\widehat{\Delta\rho_{i}}$, the relative growth
change due to intervention as described in Section 3.2, is obtained by
plugging-in the parameter estimates $\widehat{\beta}$ and $\widehat{\delta}$
into the expression $\Delta\rho_{i}(\theta)$. Since the parameter estimates
are asymptotically normal, the standard errors for $\widehat{\Delta}\rho_{i}$
can be derived using the Delta method (e.g., see Van der Vaart, (2000)),
obtaining $SE(\widehat{\Delta\rho_{0}})=\exp\\{\widehat{\beta}\\}\times
SE(\widehat{\beta})$ and
$SE(\widehat{\Delta\rho_{1}})=\exp\\{\widehat{\beta}+\widehat{\delta}\\}\times\sqrt{SE(\widehat{\beta})^{2}+SE(\widehat{\delta})^{2}+\widehat{Cov}(\widehat{\beta},\widehat{\delta})}$,
where $SE(\widehat{\beta})$ and $SE(\widehat{\delta})$ are the standard errors
for $\widehat{\beta}$ and $\widehat{\delta}$ and
$\widehat{Cov}(\widehat{\beta},\widehat{\delta})$ is an estimate of the
covariance between $\widehat{\beta}$ and $\widehat{\delta}$.

##  4 Data and Results

###  4.1 Data

The analysis in this paper focuses on the second COVID-19 outbreak wave that
occurred in Italy in late 2020. Daily data on cumulative cases and number of
tests from 1 September 2020 to 31 December 2020 on 19 Italian regions and 2
autonomous provinces were obtained from the online repository of the Italian
Civil Protection Department https://github.com/pcm-dpc/COVID-19. The
autonomous province of Bozen/Bolzano - South Tyrol is the intervention group
of our study as it was the only Italian territory that implemented mass
screening in this time period. Residents of South Tyrol were invited to take a
COVID-19 antigen rapid tests from Friday, 20 November through Sunday, 22
November; see Section 2 for details. 361,781 out of 500,607 (72.3 percent)
eligible residents volunteered to take a test.

###  4.2 Constructing the control group

The first step of our estimation process is to construct the control group of
regions against which we will evaluate the impact of the mass testing in South
Tyrol. Figure 1 shows the estimated synthetic control weights for the 20 donor
regions using the approach described in Section 2. The regions Valle d’Aosta,
Friuli Venezia Giulia and Veneto are by far the major contributors in the
synthetic control with percent weights equal to 71.0%, 20.0% and 7.5%,
respectively. Valle d’Aosta is, like South Tyrol, a small mountain region in
the North of Italy, but is not contiguous with South Tyrol. To assess the
goodness-of-fit of this selection, we computed the R-squared type statistic
$R^{2}=1-\|\mathbf{Y}_{1}-\mathbf{Y}_{0}\widehat{}\mathbf{w}\|^{2}/\|\mathbf{Y}_{1}-\bar{\mathbf{Y}}\|^{2}=0.871$,
with $\bar{\mathbf{Y}}_{1}$ denoting the arithmetic average of the elements of
vector $\mathbf{Y}_{1}$. Moreover, we found a Pearson correlation coefficient
between $\mathbf{Y}_{1}$ and $\mathbf{Y}_{0}\widehat{}\mathbf{w}$ equal to
$0.936$. Both the $R^{2}$ statistic and correlation coefficient show a very
good match between South Tyrol and the synthetic control region in the pre-
intervention period.

###  4.3 Results

Table 1 shows the estimates for the semiparametric growth models (SPGs)
described in Section 3.3 using Poisson (Pois) and negative binomial (NBin)
response functions. Each model includes an indicator for being in the
treatment region (Reg), an indicator for being in the time-period after the
mass-testing intervention (Int), and interaction between the two, and an
intrinsic growth parameter ($p$). We also estimate a second specification of
each model that includes additional (unreported) controls (+Ctrl) that could
impact short term (daily) movements of the outbreak trajectory, specifically
i) difference in the number of new tests carried out, ii) day of the week
fixed effects and iii) interactions between the number of new tests and day of
the week fixed effects.

Thus, overall we estimate four SPG models: Pois, Pois + Ctrl, NBin, and NBin +
Ctrl. For each SPG model, we report chi-squared statistics for the non-
parametric component of our model corresponding to null hypothesis
$H_{0}:h(x)=0$. For comparison purposes, Table 1 also shows in the lower part
analogous setups for pure exponential growth models obtained by setting
$h(x_{it})=0$ in the semi-parametric models. For all the considered models, we
compute the adjusted R-squared and explained deviance statistics as well as
the Akaike information criterion (AIC) and Bayesian information criterion
(BIC) for model selection.

Based on both AIC and BIC model selection criteria, the best fitting model is
NBin+Ctrl (negative binomial response with additional controls), which has an
adjusted R-squared of 88.8% and percent of explained deviance of 93.9%.
Overall, all the SPG models fit the data well with adjusted R-squared always
exceeding 83%. The standard EG models perform worse in terms of the goodness-
of-fit metrics and for all the combinations of response distribution and
predictors compared to the SPG models. The appropriateness of the semi-
parametric models is also confirmed by the statistically significant chi-
squared statistics for the non-parametric function $h$ in all considered
cases.

Figure 2 plots the estimated number of new cases $\widehat{y}_{it}$ over time
along with the actual observed data for both South Tyrol and the synthetic
control region. The plots show the results from the Poisson and negative
binomial models with and without the included additional short-term controls.
Regardless of whether additional controls are included, the fitted models show
a sizable and sudden decrease in the estimated new cases after intervention
date. The models with additional covariates fit the short- and long-term data
trajectories for both control and treatment regions remarkably well.

In each of the SPG models, we find that the difference-in-differences
interaction Int$\times$Reg, which indicates the impact of the mass testing
campaign on Covid-19 growth in South Tyrol, is negative and highly
significant. The estimated coefficient from our preferred SPG model NBin+Ctrl
is -0.512 which implies a decrease of growth rate of $39\%$ (with 95%
confidence interval 29-49%). These results are robust to both the response
function used and whether controls are included with the estimated growth rate
in these alternative specifications ranging from 38% to 51%. If instead one
uses standard EG models, the estimated impacts are generally smaller, are less
robust to model choice, and are less precisely estimated.

We next run a placebo test to gauge the possibility that our estimated impact
of the mass testing campaign occurs by chance. We do this by matching each of
the other 20 regions of Italy to synthetic versions of themselves and then
estimating the SPG model (NBin + Ctrl) assuming there was an intervention at
the same time as the real intervention in the South Tyrol. Figure 3 compares
the estimated impact in South Tyrol with the distribution of the estimates for
the other regions. Further, we examine the distribution of the estimated
interaction effect $\widehat{\delta}$ for across placebo models. The estimated
effect for the model with South Tyrol as treatment region is
$\widehat{\delta}=-0.512$ while a 95% confidence interval based on the
estimates for the other regions is $(-0.312,0.160)$. The estimate for South
Tyrol clearly stands out, with none of the other regions having large changes
in transmission rates around the time of the intervention.

Finally, to evaluate the overall impact of the mass testing campaign, we
compute the number of cases in South Tyrol over time by summing the predicted
number of cases each day after the intervention assuming that the post-
intervention transmission growth rate remained constant. Based on the
estimates from our preferred model (SPG NBin+Ctrl), there were 1,743, 2,312,
5,754 and 10,533 fewer cumulative Covid-19 cases in South Tyrol compared to a
scenario of unmitigated growth 7, 10, 20 and 40 days after the intervention,
respectively. This corresponds to an overall reduction in cases of 14%, 18%,
30% and 56%).

Our findings are in line with the predictions by Bosetti et al., (2020) which
uses a SEIR dynamic ordinary differential equations model to predict outbreak
dynamics under various contagion scenarios. When they consider the scenario
where 75% of the population is tested – which is close to the 72% observed in
South Tyrol – the predicted reduction in prevalence in the population within
10 days of mass testing is between 10% and 30% where we estimate a reduction
of 14% under the same time horizon.

| Semi-parametric growth (SPG) models:
---|---
| Pois | Pois + Ctrl | NBin | NBin + Ctrl
| Parameter estimate (SE)
Int | -0.322 (0.034)* | -0.038 (0.035) | -0.277 (0.278) | -0.078 (0.201)
Reg | 0.072 (0.011)* | 0.085 (0.011)* | 0.011 (0.063) | 0.014 (0.047)
Int $\times$ Reg | -0.754 (0.018)* | -0.569 (0.021)* | -0.745 (0.118)* | -0.512 (0.094)*
$p$ | 0.189 (0.026)* | 0.234 (0.027)* | 0.084 (0.047) | 0.106 ( 0.039)*
|  $\chi^{2}$ statistic (approx df)*
$h(x)$ | 3413 (9)* | 2669 (9)* | 267.7 (9)* | 339.8 (9)*
| Estimated growth change $\widehat{\Delta\rho}\times 100\%$ (SE)
South Tyrol | -49.47 (0.47) | -38.38 (1.09) | -52.01 (4.82) | -39.22 (5.12)
Synth Control | 7.43 (1.16) | 8.85 (1.21) | 1.13 (6.36) | 1.44 (4.78)
| Goodness-of-fit metrics
Adjusted-$R^{2}$ (%) | 83.9 | 92.2 | 82.9 | 88.8
Dev. Explained (%) | 88.9 | 94.5 | 88.1 | 93.9
AIC | 7209.8 | 4436.1 | 2637.9 | 2507.7
BIC | 7258.7 | 4530.3 | 2689.3 | 2604.8
| Exponential growth (EG) models:
| Pois | Pois + Ctrl | NBin | NBin + Ctrl
| Parameter estimate (SE)
Int | -0.830 (0.013)* | -0.704 (0.013)* | -0.561 (0.128)* | -0.426 (0.104)*
Reg | 0.063 (0.011)* | 0.072 (0.011)* | -0.001 (0.088) | 0.004 (0.070)
Int $\times$ Reg | -0.643 (0.017)* | -0.413 (0.018)* | -0.586 (0.151)* | -0.178 (0.126)
$p$ | 0.748 (0.005)* | 0.701 (0.005)* | 0.653 (0.025)* | 0.609 (0.023)*
| Estimated growth change $\widehat{\Delta\rho}\times 100\%$ (SE)
South Tyrol | -44.05 (0.73) | -28.87 (1.06) | -44.42 (6.79) | -15.96 (9.15)
Synth Control | 6.46 (1.14) | 7.51 (1.19) | -0.15 (8.81) | 0.39 (7.05)
| Goodness-of-fit metrics
Adjusted-$R^{2}$ (%) | 76.5 | 85.5 | 74.5 | 56.8
Dev. Explained (%) | 82.2 | 89.3 | 72.5 | 83.8
AIC | 10560.7 | 7013.6 | 2830.4 | 2726.5
BIC | 10578.2 | 7076.4 | 2851.3 | 2792.8
Table 1: Fitted semi-parametric growth (SPG) and exponential growth (EG) models. using Poisson (Pois) and negative binomial (NBin) response with additional control variables (+Ctrl) and without additional controls. For the parametric terms estimates are reported with standard errors in parenthesis. For the SPG models the $\chi^{2}$ statistics corresponding to null hypothesis “$H_{0}:h(x)=0$" are given with the approximate degrees of freedom (df). Significant results (p-value<0.01) are marked by “*". For all the models we report estimated percent relative growth change $\widehat{\Delta\rho}\times 100\%$ for South Tyrol and synthetic control groups as well as goodness-of-fit statistics. Figure 1: Synthetic control weights used to construct the synthetic control region. Pois | Pois + Ctrl
---|---
|
NBin | NBin +Ctrl
|
Figure 2: Fitted semi-parametric growth models against actual data on new
cases for South Tyrol and synthetic control.The first and second row
correspond, respectively, to Poisson (Pois) and negative binomial (NBin)
models for new cases. The first column corresponds to models without short-
term control factors, while the second column shows models with the additional
controls included (+Ctrl). Figure 3: Differences between new cases in
treatment and synthetic control regions over time, obtained by running the
synthetic control method for each of the 21 Italian regions. The solid line
corresponds to South Tyrol, the dashed line represents point-wise average
difference for the remaining regions, the shaded area represents points within
$1.96$ standard deviations from the point-wise average difference.

##  5 Conclusions

In this paper, we employ an innovative empirical approach to study the impact
of mass testing on Covid-19. In particular, we combine a synthetic control
methodology to define a control group and a semi-parametric growth model to
model the epidemiological developments. We show that using a semi-parametric
approach rather than more standard parametric exponential growth models makes
a difference in the evaluation of the impact of the mass testing that took
place in the Italian region of South-Tyrol in November 2020 on Covid-19
growth.

In our preferred specification, we find that mass testing reduced the growth
rate of Covid-19 by 39%. This suggests that mass testing can be an useful tool
to contain the pandemic. Importantly, the mass testing campaign in South Tyrol
was entirely voluntary and still managed to attract 72 percent of eligible
residents to participate. This likely occurred because it was underlined that
participating would help reduce the spread of Covid-19 and it could
potentially lead to fewer restrictions in other dimensions. In many countries,
public testing is openly available, suggesting that supply constraints are
less and less binding. We find that a more organized approach to testing can
make a big difference in reducing the spread of Covid-19.

## References

  * (1) Abadie, A., Diamond, A., and Hainmueller, J. (2010a).  Synthetic control methods for comparative case studies: Estimating the effect of california’s tobacco control program.  Journal of the American Statistical Association, 105(490):493–505.
  * (2) Abadie, A., Diamond, A., and Hainmueller, J. (2010b).  Synthetic control methods for comparative case studies: Estimating the effect of california’s tobacco control program.  Journal of the American Statistical Association, 105(490):493–505.
  * Abadie et al., (2011) Abadie, A., Diamond, A., and Hainmueller, J. (2011).  Synth: An r package for synthetic control methods in comparative case studies.  Journal of Statistical Software, Articles, 42(13):1–17.
  * (4) Abadie, A. and Gardeazabal, J. (2003a).  The economic costs of conflict: A case study of the basque country.  American economic review, 93(1):113–132.
  * (5) Abadie, A. and Gardeazabal, J. (2003b).  The economic costs of conflict: A case study of the basque country.  American Economic Review, 93(1):113–132.
  * Alexander et al., (2020) Alexander, A., Martin, H., and Lackner, M. (2020).  Mass gatherings contributed to early covid-19 spread: Evidence from us sports.  Technical report, Working paper.
  * Atkeson et al., (2021) Atkeson, A., Droste, M. C., Mina, M. J., and Stock, J. (2021).  Economic benefits of covid-19 screening tests with a vaccine rollout.  medRxiv.
  * Aviv-Sharon and Aharoni, (2020) Aviv-Sharon, E. and Aharoni, A. (2020).  Generalized logistic growth modeling of the covid-19 pandemic in asia.  Infectious Disease Modelling, 5:502–509.
  * Bosetti et al., (2020) Bosetti, P., Tran Kiem, C., Yazdanpanah, Y., Fontanet, A., Lina, B., Colizza, V., and Cauchemez, S. (2020).  Impact of mass testing during an epidemic rebound of sars-cov-2: A modelling study.  medRxiv.
  * Chénangnon et al., (2020) Chénangnon, F. T., Lokonon, B. E., and Kakaï, R. G. (2020).  On the use of growth models to understand epidemic outbreaks with application to covid-19 data.  PLoS One, 15(10):e0240578.
  * Cho, (2020) Cho, S.-W. S. (2020).  Quantifying the impact of nonpharmaceutical interventions during the COVID-19 outbreak: The case of Sweden.  The Econometrics Journal, 23(3):323–344.
  * Chowell, (2017) Chowell, G. (2017).  Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts.  Infectious Disease Modelling, 2(3):379–398.
  * Dave et al., (2021) Dave, D., Friedson, A. I., Matsuzawa, K., and Sabia, J. J. (2021).  When do shelter-in-place orders fight covid-19 best? policy heterogeneity across states and adoption time.  Economic Inquiry, 59(1):29–52.
  * Hasell et al., (2020) Hasell, J., Mathieu, E., Beltekian, D., Macdonald, B., Giattino, C., Ortiz-Ospina, E., Roser, M., and Ritchie, H. (2020).  A cross-country database of covid-19 testing.  Scientific Data, 7(1):345.
  * Hastie and Tibshirani, (1990) Hastie, T. J. and Tibshirani, R. J. (1990).  Generalized additive models, volume 43.  CRC press.
  * Haug et al., (2020) Haug, N., Geyrhofer, L., Londei, A., Dervic, E., Desvars-Larrive, A., Loreto, V., Pinior, B., Thurner, S., and Klimek, P. (2020).  Ranking the effectiveness of worldwide covid-19 government interventions.  Nature Human Behaviour, 4(12):1303–1312.
  * Holt, (2020) Holt, E. (2020).  Slovakia to test all adults for sars-cov-2.  The Lancet, 396(10260):1386–1387.
  * Hsiang et al., (2020) Hsiang, S., Allen, D., Annan-Phan, S., Bell, K., Bolliger, I., Chong, T., Druckenmiller, H., Huang, L. Y., Hultgren, A., Krasovich, E., et al. (2020).  The effect of large-scale anti-contagion policies on the covid-19 pandemic.  Nature, 584(7820):262–267.
  * Kahanec et al., (2021) Kahanec, M., Lafférs, L., and Schmidpeter, B. (2021).  The impact of mass antigen testing for covid-19 on the prevalence of the disease.  GLO Discussion Paper 775, Essen.
  * Lee et al., (2020) Lee, S. Y., Lei, B., and Mallick, B. K. (2020).  Estimation of covid-19 spread curves integrating global data and borrowing information.  arXiv preprint arXiv:2005.00662.
  * Mangrum and Niekamp, (2020) Mangrum, D. and Niekamp, P. (2020).  Jue insight: College student travel contributed to local covid-19 spread.  Journal of Urban Economics, page 103311.
  * Mitze et al., (2020) Mitze, T., Kosfeld, R., Rode, J., and Wälde, K. (2020).  Face masks considerably reduce covid-19 cases in germany.  Proceedings of the National Academy of Sciences, 117(51):32293–32301.
  * Pavelka et al., (2021) Pavelka, M., Van-Zandvoort, K., Abbott, S., Sherratt, K., Majdan, M., Analýz, I. Z., Jarčuška, P., Krajčí, M., Flasche, S., and Funk, S. (2021).  The impact of population-wide rapid antigen testing on sars-cov-2 prevalence in slovakia.  Science.
  * Pelinovsky et al., (2020) Pelinovsky, E., Kurkin, A., Kurkina, O., Kokoulina, M., and Epifanova, A. (2020).  Logistic equation and covid-19.  Chaos, Solitons & Fractals, 140:110241.
  * Singh et al., (2021) Singh, S., Shaikh, M., Hauck, K., and Miraldo, M. (2021).  Impacts of introducing and lifting nonpharmaceutical interventions on covid-19 daily growth rate and compliance in the united states.  Proceedings of the National Academy of Sciences, 118(12).
  * Stillman and Tonin, (2021) Stillman, S. and Tonin, M. (2021).  Communities and testing for covid-19.
  * Tian et al., (2020) Tian, H., Liu, Y., Li, Y., Wu, C.-H., Chen, B., Kraemer, M. U., Li, B., Cai, J., Xu, B., Yang, Q., et al. (2020).  An investigation of transmission control measures during the first 50 days of the covid-19 epidemic in china.  Science, 368(6491):638–642.
  * Van der Vaart, (2000) Van der Vaart, A. W. (2000).  Asymptotic statistics, volume 3.  Cambridge university press.
  * Wood and Wood, (2015) Wood, S. and Wood, M. S. (2015).  Package ‘mgcv’.  R package version, 1:29.
  * Wood, (2017) Wood, S. N. (2017).  Generalized additive models: an introduction with R.  CRC press.
  * Wu et al., (2020) Wu, K., Darcet, D., Wang, Q., and Sornette, D. (2020).  Generalized logistic growth modeling of the covid-19 outbreak in 29 provinces in china and in the rest of the world.  arXiv preprint arXiv:2003.05681.


# Hydrodynamic Simulations of Electron-capture Supernovae: Progenitor and
Dimension Dependence

Shuai Zha,1 Evan P. O’Connor,2 Sean M. Couch,3,4,5 Shing-Chi Leung,6 and
Ken’ichi Nomoto,7
1Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai 200240,
People’s Republic of China
2The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-106 91 Stockholm, Sweden
3Department of Physics and Astronomy, Michigan State University, East Lansing,
Michigan 48824, USA
4Department of Computational Mathematics, Science, and Engineering, Michigan
State University, East Lansing, Michigan 48824, USA
5Facility for Rare Isotope Beams, Michigan State University, East Lansing, MI
48824, USA
6TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA
91125, USA
7Kavli Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study,
The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
<EMAIL_ADDRESS>
(Accepted 09 April 2022. Received 20 March 2022; in original form 23 December
2021)

###### Abstract

We present neutrino-transport hydrodynamic simulations of electron-capture
supernovae (ECSNe) in FLASH with new two-dimensional (2D) collapsing
progenitor models. These progenitor models feature the 2D modelling of oxygen-
flame propagation until the onset of core collapse. We perform axisymmetric
simulations with 6 progenitor models that, at the time of collapse, span a
range of propagating flame front radii. For comparison, we also perform a
simulation with the same setup using the canonical, spherically-symmetrical
progenitor model n8.8. We found that the variations in the progenitor models
inherited from simulations of stellar evolution and flame propagation do not
significantly alter the global properties of the neutrino-driven ECSN
explosion, such as the explosion energy ($\sim 1.36$–$1.48\times 10^{50}$ erg)
and the mass ($\sim 0.017$–$0.018M_{\odot}$) and composition of the ejecta.
Due to aspherical perturbations induced by the 2D flame, the ejecta contains a
small amount ($\lesssim 1.8\times 10^{-3}~{}M_{\odot}$) of low-$Y_{e}$
($0.35<Y_{e}<0.4$) component. The baryonic mass of the protoneutron star is
$\sim 1.34~{}M_{\odot}$ ($\sim 1.357~{}M_{\odot}$) with the new (n8.8)
progenitor models when simulations end at $\sim 400$ ms and the discrepancy is
due to updated weak-interaction rates in the progenitor evolutionary
simulations. Our results reflect the nature of ECSN progenitors containing a
strongly degenerate ONeMg core and suggest a standardized ECSN explosion
initialized by ONeMg core collapse. Moreover, we carry out a rudimentary
three-dimensional simulation and find that the explosion properties are fairly
compatible with the 2D counterpart. Our paper facilitates a more thorough
understanding of ECSN explosions following the ONeMg core collapse, though
more three-dimensional simulations are still needed.

###### keywords:

(stars:) supernovae: general, hydrodynamics, stars: neutron

††pubyear: 2021††pagerange: Hydrodynamic Simulations of Electron-capture
Supernovae: Progenitor and Dimension Dependence–References

##  1 Introduction

Electron-capture supernovae (ECSNe) have been proposed as the complementary
channel for producing neutron stars (NS; see Miyaji et al., 1980; Nomoto,
1984), other than the canonical channel, i.e. core-collapse supernova (CCSN;
see the recent review by Burrows & Vartanyan, 2021). This additional channel
can explain some bi-modality in the properties of observed NSs, such as the
spin-orbit period (Knigge et al., 2011) and mass distribution (Schwab et al.,
2010; Farrow et al., 2019). Moreover, NSs born in ECSNe are believed to obtain
a small kick velocity, $\lesssim 100$ km/s (Podsiadlowski et al., 2004;
Gessner & Janka, 2018; Stockinger et al., 2020), so the ECSN channel can be
rather important for making compact binaries with a NS component (Kruckow et
al., 2018).

Progenitors of ECSNe have a zero-age-main-sequence (ZAMS) mass of
$\sim$8–10$~{}M_{\odot}$, lying just below that of CCSN progenitors (Nomoto &
Hashimoto, 1988; Woosley et al., 2002). Therefore, ECSNe can be abundant due
to the bottom-heavy initial mass function of stars (Salpeter, 1955) and thus
contribute to a significant portion of supernova-like transients. Depending on
the mass loss during the last evolutionary stage of their progenitors, ECSNe
may lead to normal to low-luminosity type IIp supernovae (Botticella et al.,
2009; Tominaga et al., 2013; Kozyreva et al., 2021; Hiramatsu et al., 2021),
Type IIn supernovae (Moriya et al., 2014), fast evolving luminous transients
(Tolstov et al., 2019), and intermediate-luminosity red transients (Cai et
al., 2021).

Unlike a CCSN, an ECSN progenitor contains a degenerate oxygen-neon-magnesium
(ONeMg) core at its final evolutionary stage (Nomoto & Leung, 2017). An oxygen
flame develops due to the heating by electron captures on 24Mg and 20Ne, and
during its outward propagation as a deflagration wave, electron captures on
the nuclear-statistical-equilibrium (NSE) ash behind the flame front lead to
the collapse of the ONeMg core. However, from the beginning of ECSN studies,
it has been debated in theory that the oxygen flame may as well lead to a
(partially) thermonuclear explosion of ONeMg core, similar to that in a Type
Ia supernova (Nomoto & Kondo, 1991; Timmes & Woosley, 1992; Jones et al.,
2016; Leung & Nomoto, 2019; Zha et al., 2019; Kirsebom et al., 2019).

The collapse-explosion puzzle of ECSNe cannot be solved solely from the
theoretical side, due to the subtleties involved in simulations, such as the
weak-reaction rates (Suzuki et al., 2019; Strömberg et al., 2021; Langanke et
al., 2021), convective URCA process (Paczyński, 1973; Denissenkov et al.,
2015; Schwab et al., 2017), convective mixing process driven by oxygen burning
(Nomoto, 1987; Zha et al., 2019; Takahashi et al., 2019), and physics of the
oxygen flame (Leung et al., 2020; Schwab et al., 2020). Each scenario, i.e.
collapse or thermonuclear explosion, needs to be examined carefully with
sophisticated simulations, which make solid predictions to confront the
existing and upcoming observations.

For the collapse scenario of ECSN, hydrodynamic simulations of the collapse
and explosion phases (Kitaura et al., 2006; Janka et al., 2008; Fischer et
al., 2010; Radice et al., 2017; Stockinger et al., 2020) have been using the
spherically symmetrical progenitor model n8.8 provided in Nomoto (1984, 1987).
This canonical ECSN progenitor model has led to substantial knowledge on the
explosion properties of ECSNe. However, important updates have been made to
the weak-interaction rates (Toki et al., 2013; Martínez-Pinedo et al., 2014;
Suzuki et al., 2016; Kirsebom et al., 2019; Suzuki et al., 2019) and the
numerical schemes for the flame propagation, e.g. more than 1D modelling
(Jones et al., 2016; Leung et al., 2020). In this work, we perform a suite of
2D axisymmetric neutrino-transport simulations of ECSNe with the updated
progenitor models possessing a collapsing ONeMg core provided in Zha et al.
(2019). These progenitor models feature the updated weak-interaction rates and
2D simulations for the propagation of the oxygen flame. We focus on how the
uncertainties in the progenitor evolution transforms into the explosion
properties of ECSNe. Moreover, we carry out a rudimentary 3D octant simulation
to examine the 3D effects. Our results show that unlike CCSNe whose
progenitors are more massive, the 3D ECSN simulation exhibits only small
differences compared to its 2D axisymmetric counterpart.

Our paper is organized as follows. In Section 2 we describe the ECSN
progenitor models and methods used in our numerical simulations. Results of
progenitor and dimension dependence are presented in Section 3 and conclusions
are given in Section 4.

##  2 Methods

###  2.1 Progenitor models

We use the collapsing ECSN progenitor models of Zha et al. (2019). Here we
briefly describe their evolution prior to collapse for the sake of
completeness, and refer readers to Zha et al. (2019) and Leung et al. (2020)
for details.

A solar-metallicity star with a ZAMS mass of $8.4~{}M_{\odot}$ was evolved
with the MESA code (Paxton et al., 2011; Paxton et al., 2013, 2015, 2018,
2019) from the main-sequence phase to the formation of a degenerate ONeMg
core. The simulation of the whole star was terminated when its helium (He)
layer becomes thin ($\simeq 0.2M_{\odot}$). Afterwards the hydrogen (H)
envelope and He layer were removed and the ONeMg core was evolved with a
constant mass accretion rate, $10^{-6}$ or $10^{-7}$ $M_{\odot}$ yr-1, to
mimic the H-He thermal pulses. During the growth of the ONeMg core, two
extreme criteria, i.e. Ledoux and Schwarzschild (Miyaji & Nomoto, 1987;
Kippenhahn et al., 2012), were used to explore the uncertain semiconvection
driven by the heating due to electron captures on 24Mg and 20Ne. The MESA
calculation was terminated at the ignition of oxygen burning, when its energy
release exceeds that of thermal neutrino losses. The ONeMg core has a central
density of $\sim 10^{10}$ g cm-3 at this moment.

The oxygen burning will eventually turn into nuclear runaway and initiation of
a deflagration wave that propagates the oxygen-flame front. On the other hand,
oxygen burning can drive convection to allow energy being transported away
from the burning region and thus defer the nuclear runaway to a higher central
density. Two approaches were used to follow the evolution after the oxygen
ignition. Firstly, the nuclear runaway was assumed to happen immediately after
the oxygen ignition. Secondly, to accommodate the effect of the convection
driven by oxygen burning, the last MESA profile was modified with a fully
mixed central region of $\sim 0.14~{}M_{\odot}$, while the central density was
parameterized with a larger value between $10^{10}$ to $10^{10.2}$ g cm-3.

The last MESA profile was then mapped into a 2D cylindrical grid with an
initial flame sitting at the spot of oxygen burning. A 2D hydrodynamic code
with the capability of flame-capturing (Leung et al., 2015, 2020) was used to
follow the propagation of the oxygen flame and electron captures on the NSE
ash behind it. The simulations were terminated when the ONeMg core reaches a
central density of $3\times 10^{10}$ g cm-3, signaling its collapse.

Figure 1: Profiles of density, temperature, electron fraction, and transverse velocity ($v_{\theta}$) for the progenitor models s1, m1 and l1. The red line in each panel marks the locus of the oxygen-flame front in the progenitor models. Note that the models are only in the first quadrant, and the variables except density are mirrored to other quadrants. Model | $M_{\rm core}$ | $R_{\rm core}$ | $v_{\rm infall}$ | $M_{\rm NSE}$ | $R_{\rm flame,\max}$ | $R_{\rm flame,\min}$ | $\langle R_{\rm flame}\rangle$ | $E_{\rm kin,\theta}$ | Nomenclature in | Mixing scheme
---|---|---|---|---|---|---|---|---|---|---
| ($M_{\odot}$) | (km) | ($10^{8}$ cm s-1) | ($M_{\odot}$) | (km) | (km) | (km) | ($10^{46}$erg) | Zha et al. (2019) |
n8.8 | 1.377 | 1090 | 0.8 | 0.090 | 110 | 110 | 110 | N.A. | N.A. | N.A.
s1 | 1.361 | 1190 | 1.6 | 0.270 | 219 | 158 | 198 | 7.09 | 6-1020-046-30 | L_no_mix
s2 | 1.361 | 1180 | 1.6 | 0.250 | 208 | 166 | 191 | 3.65 | 7-1010-046-00 | L_no_mix
m1 | 1.359 | 1210 | 2.4 | 0.539 | 338 | 212 | 308 | 41.7 | 6-0999-049-00 | S_$\rho$_mix
m2 | 1.358 | 1180 | 2.1 | 0.498 | 319 | 205 | 294 | 32.6 | 7-0997-049-00 | S_$\rho$_mix
l1 | 1.358 | 1200 | 3.3 | 0.790 | 460 | 254 | 363 | 117 | 6-1000-049-00 | S_$\rho$_mix
l2 | 1.359 | 1210 | 3.0 | 0.709 | 426 | 209 | 341 | 108 | 6-0999-049-30-LM | Ledoux mix o-burn
Table 1: Properties of the progenitor models used in this work. $M_{\rm core}$
and $R_{\rm core}$ are the mass and radius of the ONeMg core. $v_{\rm infall}$
is the maximum infall velocity in the ONeMg core at collapse. $M_{\rm NSE}$ is
the mass of the NSE ash enclosed by the flame front. $R_{\rm flame,\max}$,
$R_{\rm flame,\min}$ and $\langle R_{\rm flame}\rangle$ are the maximum,
minimum and angular-averaged radius of the flame front. $E_{\rm kin,\theta}$
is the total transverse kinetic energy of the progenitor models. The
corresponding nomenclature and mixing scheme used in Section 3.2 of Zha et al.
(2019) are given for clarity. “L_no_mix" and “S_$\rho$_mix" refer to MESA
simulations using the Ledoux and Schwarzschild criteria for the convective
instability, respectively. They both assumes nuclear runaway occurs
immediately after the oxygen ignition. “Ledoux mix o-burn" uses the Ledoux
criterion and considers the convection driven by the oxygen burning.

Zha et al. (2019) found that the configuration of the ONeMg core at the oxygen
runaway determines the subsequent flame propagation. In general, a higher
central density and a more mixed central region favor electron captures,
accelerate core contraction, and result in a shorter distance of the flame
propagation until the onset of core collapse. Here we further investigate how
the propagated distance of the oxygen flame will affect properties of ECSN
explosions. In Table 1 we list the progenitor models of Zha et al. (2019)
employed in this work. The model names indicate the size of the region
enclosed by the flame front, with s, m and l corresponding to small, medium
and large. Our nomenclature for the models is different from that in Section
3.2 of Zha et al. (2019), which is also given in Table 1 for clarity. We
choose 2 models for each group of flame size and the mean radius of the flame
front in the model ‘2’ is smaller by $\sim 10-20$ km than that in the model
‘1’. We diversify the progenitor models with respect to the parameters used in
the stellar evolutionary calculations (Zha et al., 2019). For example, in
groups ‘s’ and ‘m’, the model ‘1’ uses a mass accretion rate of
$10^{-6}~{}M_{\odot}~{}{\rm yr}^{-1}$ and the model ‘2’ uses a mass accretion
rate of $10^{-7}~{}M_{\odot}~{}{\rm yr}^{-1}$. These models cover the extremes
of propagation of the oxygen flame and reflect the uncertainties in the
progenitor evolutionary simulations. We note that the mass of the NSE ash
enclosed by the flame front ranges from $18\%$ to $58\%$ of the ONeMg core.
The models used here are publicly available at
https://doi.org/10.5281/zenodo.5748457. We also perform simulations with the
progenitor model n8.8 (Nomoto, 1984, 1987) for comparison. One can see from
Table 1 that n8.8 has the smallest mass of NSE ash and radius of the flame
front.

The 2D profiles of density, temperature, electron fraction ($Y_{e}$), and
transverse velocity ($v_{\theta}$) are rendered in Fig. 1 for progenitor
models s1, m1 and l1. The red line in each plot marks the locus of the oxygen-
flame front in the progenitor models.

###  2.2 Neutrino-transport simulations

We utilize the FLASH code (Fryxell et al., 2000) for the neutrino-transport
hydrodynamic simulations that follow the collapse of the ONeMg core and the
subsequent explosion. The equation of state (EOS) of the stellar matter is
described by the SFHo EOS (Steiner et al., 2013) at densities above $4\times
10^{6}$ g cm-3 and the Helmholtz EOS (Timmes & Swesty, 2000) at densities
below $2\times 10^{6}$ g cm-3, with a linear interpolation for smooth
transition in between (Witt et al., 2021). Neutrino transport with full
velocity dependence is simulated with a two-moment scheme and the analytic
“M1" closure (Shibata et al., 2011; Cardall et al., 2013; O’Connor, 2015;
O’Connor & Couch, 2018a). Neutrino-matter interactions are described by NuLib
(O’Connor, 2015), with 12 energy groups from 0 to 300 MeV spaced
logarithmically. We use the same baseline prescription for neutrino-matter
interactions as that of (O’Connor et al., 2018). We note explicitly that these
FLASH simulations, unlike the FLASH simulations in (O’Connor et al., 2018),
include neutrino-electron inelastic scattering.

We use a block structured adaptive-mesh-refinement grid for our simulations.
In the 2D simulations, the cylindrical grid extends to $1.2\times 10^{4}$ km
and includes the first quadrant, assuming reflective symmetry at the
equatorial plane as in the previous simulations of the deflagration phase (Zha
et al., 2019). Our fiducial resolution is $\sim$250 m at $\lesssim 40$ km to
resolve the protoneutron star (PNS) and maintains an effective angular
resolution of 0.6∘ at $\gtrsim 40$ km. In the 3D octant simulation, we use a
coarser resolution, $\sim$500 m at the center and $1.5^{\circ}$ outside $40$
km. To disentangle the effects of resolution and dimension, we run a 2D
simulation (m1-2d-lr) with a comparable resolution as the 3D run (m1-3d-lr).
The timestep is smaller in 3D (0.43 $\mu$s) than that in 2D (0.65 $\mu$s) to
maintain numerical stability of the neutrino transport. The results of
integrated quantities such as the explosion energy and mass of ejecta are
multiplied by 2 in 2D and 8 in 3D according to the imposed symmetries. In
addition, we perform two 1D simulations with the fiducial resolution using
n8.8 (n8.8-1d) and spherically-averaged m1 (m1-1d) as the progenitor models.

In the simulations of the deflagration phase the computational domain only
included the ONeMg core. To follow the evolution of shock after its breaking
out of the core surface, we attach a convective envelope above the core with
$X_{\rm H}=0.6$ and $X_{\rm He}=0.4$, where $X_{i}$ is the mass fraction of
element $i$. The envelope has a mass of $\lesssim 10^{-5}~{}M_{\odot}$ and
should not affect the results of the ECSN explosion, such as the explosion
energy and the properties of the ejecta (Janka et al., 2008).

##  3 Results

In this section, we first present a systematic investigation of the progenitor
dependence with 2D axisymmetric simulations. Then we compare the results of
the 2D and 3D octant simulations with the same progenitor and show the impact
of the 3rd dimension. We further compare the properties of the ejecta in
different simulations and discuss the implications for nucleosynthesis of
ECSNe. Finally we discuss the caveat that our neutrino-transport simulations
lacks a nuclear network to simulate the flame and burning of oxygen and neon.

###  3.1 Axisymmetric simulations and progenitor dependence

Figure 2: Mean shock radius ($\langle R_{\rm shock}\rangle$, upper panel, thick lines), mass of ejecta ($M_{\rm ej}$, upper panel, thin lines), and diagnostic explosion energy ($E_{\rm exp}$, lower panel) as a function of the postbounce time in the axisymmetric simulations. Model | $t_{\rm b}$ | $E_{\rm exp}$ | $M_{\rm ej}$ | $M_{\rm PNS}$
---|---|---|---|---
| (ms) | $(10^{51}\rm erg)$ | $(10^{-2}M_{\odot})$ | $(M_{\odot})$
n8.8-1d | 62.2 | 0.111 | 1.39 | 1.362
m1-1d | 63.8 | 0.109 | 1.35 | 1.344
n8.8 | 62.8 | 0.146 | 1.74 | 1.357
s1 | 68.5 | 0.141 | 1.72 | 1.341
s2 | 69.7 | 0.148 | 1.79 | 1.340
m1 | 67.2 | 0.143 | 1.73 | 1.339
m2 | 66.5 | 0.139 | 1.69 | 1.339
l1 | 75.8 | 0.136 | 1.69 | 1.340
l2 | 71.1 | 0.136 | 1.70 | 1.340
m1-2d-lr | 67.2 | 0.188 | 2.25 | 1.333
m1-3d-lr | 67.1 | 0.177 | 2.06 | 1.334
Table 2: Results of the hydrodynamic simulations. n8.8-1d and m1-1d are the 1D
simulations with the progenitor models n8.8 and spherically-averaged m1,
respectively. s1 to l2 are the 2D axisymmetric simulations with the
corresponding progenitor models. These are with the fiducial resolution.
m1-2d-lr and m1-3d-lr are the low-resolution simulations with the progenitor
model m1 in 2D and 3D, respectively. $t_{\rm b}$ is time of bounce. $E_{\rm
exp}$, $M_{\rm ej}$ and $M_{\rm PNS}$ are the diagnostic explosion energy,
mass of ejecta and baryonic mass of protoneutron star at the end of
simulations, which is 400 ms after bounce except for m1-3d-lr and m1-2d-lr
which are at 300 ms after bounce).

Due to the high initial central density ($3\times 10^{10}$ g cm-3) of ECSN
progenitors, the collapse of ONeMg core lasts for only $\sim$60-70 ms until
the formation of the protoneutron star (PNS) and the bounce shock wave. The
bounce time, $t_{\rm b}$, as listed in Table 2, differs by 10 ms at most for
2D progenitors and shows a non-monotonic trend with the size of flame front
due to the competition between electron captures and thermal effects. Model
n8.8 has a central density of $4.43\times 10^{10}$ g cm-3 and the smallest
$t_{\rm b}$ of 62.8 ms. As plotted in the upper panel of Fig. 2, the bounce
shock reaches a radius of $100$ km in $\sim 20$ ms, and then gradually expands
until breaking out of the core surface at $\sim 450$ km and $\sim 80$ ms
postbounce. Due to the steep density gradient at the core surface, the shock
gets accelerated there and reaches 1000 km in another 20 ms.

In the lower panel of Fig. 2 we plot the instantaneous diagnostic explosion
energy, defined as

$E_{\rm exp}=\int_{e_{\rm tot}>0,v_{r}>0}dV\rho e_{\rm tot},$ (1)

where $\rho$ and $v_{r}$ are the matter density and radial velocity of a fluid
element. $e_{\rm tot}$ is the specific total energy

$e_{\rm tot}=e-e_{\rm cold}+\Phi+0.5|\vec{v}|^{2},$ (2)

where $e$, $\Phi$ and $0.5|\vec{v}|^{2}$ are the specific internal,
gravitational and kinetic energy. $e_{\rm cold}$ is the the specific internal
energy of matter with the same density and $Y_{e}$ of the fluid element, but
with temperature set to be the minimum value of the SFHo EOS table (0.01 MeV).
Once the shock breaks out of the core surface, $E_{\rm exp}$ quickly becomes
nonzero due to neutrino heating and gradually increases to $\sim 1.3\times
10^{50}$ erg in the following $\sim 200$ ms. Afterwards $E_{\rm exp}$
increases slowly to $\sim 1.4\times 10^{50}$ erg at the end of simulations and
differs by $0.12\times 10^{50}$ erg ($\sim$9$\%$) at most among all progenitor
models (cf. Table 2).

The mass of the ejecta, defined as

$M_{\rm ej}=\int_{e_{\rm tot}>0,v_{r}>0}dV\rho,$ (3)

is $\sim$0.017-0.018$~{}M_{\odot}$ at the end of the simulations (thin lines
in the upper panel of Fig. 2). We note that due to the shock breakout of the
outer boundary of our computational domain at $\sim 150$ ms postbounce, a
small amount of matter with mass$\sim 2\times 10^{-4}~{}M_{\odot}$ and $E_{\rm
exp}\sim 3\times 10^{48}$ erg has left the grid by the end of the simulations.
This accounts for $\sim 1\%$ of $M_{\rm ej}$ and $\sim 2\%$ of $E_{\rm exp}$
which we add back in Fig. 2 and Table 2. As shown in Fig. 2 and Table 2, the
explosion properties show a modest dependence on the ECSN progenitor, while 2D
simulations have a larger explosion energy and ejecta mass by 25%-30%
comparing to their 1D counterparts (n8.8-1d and m1-1d). In Section 3.3 we will
further show the detailed composition of the ejecta and the implications for
ECSN nucleosynthesis.

Figure 3: Colormaps depicting the magnitude of angular averaged transverse
velocity $\langle|v_{\theta}|\rangle$ as a function of radius and the
postbounce time in model n8.8 (left) and m1 (right). Black lines indicate the
radii that enclose specific masses. The mass resolution in the left (right)
plot is 0.1 $M_{\odot}$, 0.01 $M_{\odot}$ and 0.001 $M_{\odot}$ below $M=$1.3
(1.3), 1.36 (1.34) and 1.377 (1.359) $M_{\odot}$, which are indicated by the
blue dashed lines. The white, red and yellow lines mark the mean shock radius,
mean gain radius and PNS surface; the radius of the PNS surface is defined as
the mean radius with a density of $10^{11}~{}{\rm g~{}cm}^{-1}$.

Recent CCSN simulations found that precollapse aspherical perturbations in the
convective oxygen and silicon shells (Fields & Couch, 2020, 2021; Müller et
al., 2017; Vartanyan et al., 2022) are amplified when crossing the standing
accretion shock and induce violent turbulent convection which can enhance the
neutrino energy deposition (see e.g. Couch & Ott, 2013; Müller et al., 2017;
O’Connor & Couch, 2018b; Vartanyan et al., 2018). Our progenitor models also
possess such perturbations inside and around the oxygen-flame front (see
bottom right corner of Fig. 1). To see how this affects the postbounce
turbulence, in Fig. 3 we plot the angular-averaged transverse velocity
$\langle|v_{\theta}|\rangle$ as a function of radius and the postbounce time
in models n8.8 and m1. As the bounce shock is launched the postshock region
gains a $\langle|v_{\theta}|\rangle$ up to $\sim 3\times 10^{8}$ cm s-1 in
both models. This is the so-called prompt convection due to the negative
entropy gradient behind the shock (Bruenn & Mezzacappa, 1994). This prompt
convection spans from $\sim 10$ km to the shock and lasts for $\sim 20$ ms.
Model m1 shows a slightly stronger turbulence than model n8.8 and this is
compared quantitatively later. After this episode, convection develops in the
neutrino-heated postshock region and below the PNS surface, similar to that in
the simulations of massive progenitors. Here, we define the PNS surface as the
mean radius with a density of $10^{11}~{}{\rm g~{}cm}^{-3}$. The magnitudes of
$\langle|v_{\theta}|\rangle\sim 1.5\times 10^{8}$ cm s-1 are quantitatively
similar in the two models for both convective regions. We note a similar
magnitude of $\langle|v_{\theta}|\rangle$ just below the fast expanding shock
wave (white lines in Fig. 3) especially for model n8.8. The transverse kinetic
energy of this fluid element is small ($\sim 1\%$) comparing to that in the
gain region due to its small mass ($\sim 10^{-4}~{}M_{\odot}$).

Figure 4: Left panel: Total transverse kinetic energy of the matter inside the
shock as a function of the postbounce time in the 0.05 s after bounce. Right
panel: Same as the left panel, but for the matter in the gain region from 0.05
s to 0.4 s after bounce.

It is more instructive to look at the difference in the total kinetic energy
of the transverse motion ($E_{\rm kin,\theta}$) in turbulent regions, which
serves as a proxy for the overall turbulence. Because in the early postbounce
phase the prompt convection spans from $\sim 10$ km to the shock, we use
$E_{\rm kin,total,\theta}$ for matter inside the shock as a measure of
turbulence in the 0.05 s after bounce (left panel of Fig. 4). It manifests the
impact of precollapse asphericities due to the oxygen flame on the prompt
convection in the postshock region. Generally, a larger flame front leads to
stronger turbulence in the postshock region when the aspherical perturbations
pass through the shock. While models with a small flame front have the same
strength of prompt convection as model n8.8, models with a large flame front
acquire a larger turbulent energy by a factor of 2. In the left panel of Fig.
4, the slightly larger $E_{\rm kin,\theta}$ seen in the model n8.8 compared
with models s1 and s2 suggests that the turbulent power seeded by the
numerical perturbation of our cylindrical grid is at least comparable to that
seeded by asphericities of the flame with the smallest propagated radius.

At later times when convection develops rigorously in the PNS surface region
($t_{\rm pb}\gtrsim$0.05 ms), one can separate $E_{\rm kin,\theta}$ in the
gain region and the PNS convective region. In the right panel of Fig. 4 we
plot $E_{\rm kin,gain,\theta}$ for the matter in the gain region for $t_{\rm
pb}=$0.05-0.4 s. All simulations with 2D progenitors still possess a larger
$E_{\rm kin,\theta}$ than that with model n8.8 but there is no clear relation
between $E_{\rm kin,gain,\theta}$ and initial flame size. This is likely due
to a joint effect of neutrino heating and turbulence. After $t_{\rm pb}\approx
0.25$ s, $E_{\rm kin,\theta}$ becomes similar in all simulations as the
neutrino heating becomes weak.

Figure 5: Same as Fig. 2, but for the mass (top panel), radius (middle panel)
and turbulent kinetic energy (bottom panel) of the protoneutron star as a
function of the postbounce time in axisymmetric simulations. Note that at
early postbounce time $E_{\rm kin,PNS,\theta}$ includes contribution of prompt
convection.

In Figure 5, we plot the mass, radius, and transverse kinetic energy of the
PNS as a function of the postbounce time. When the simulations end, the PNS
baryonic mass differs by at most $\sim 0.002~{}M_{\odot}$ among 2D progenitor
models, which reflects the same difference in the ONeMg core mass (cf. Table
1). This means that non-rotating ECSNe produce a NS with a generic baryonic
mass of $\sim 1.34~{}M_{\odot}$, irrespective of the uncertainty in the flame
physics. The corresponding gravitational mass is $\sim 1.23~{}M_{\odot}$ for a
cold and $\beta$-equilibrated neutron star with the same baryonic mass and the
SFHo EOS. The updated weak interaction rates used in Zha et al. (2019) lead to
a lower $Y_{e}$ and a smaller mass of the ONeMg core than that of n8.8. The
latter results in a PNS baryonic (gravitational) mass of 1.357 (1.245)
$M_{\odot}$111We note that Radice et al. (2017) determined a final PNS
baryonic mass of $1.30~{}M_{\odot}$ with the n8.8 progenitor. This is due to
the misinterpretion of “radius” in the data file which results in an ONeMg
core mass of $1.32~{}M_{\odot}$ (Radice, 2022). The “radius” in n8.8 is given
at the inner edge of each mass zone.. The PNS radius has almost no difference
after $\sim$200 ms postbounce and the transverse kinetic energy shows a small
difference among all simulations. The strength of PNS convection, as measured
by $E_{\rm kin,PNS,\theta}$, reaches a maximum value of 4-5$\times 10^{48}$
erg at 0.2 s postbounce, which is similar to that in low-mass CCSN models but
about an order of magnitude smaller than more massive progenitors (Nagakura et
al., 2020; Eggenberger Andersen et al., 2021) due to the absence of PNS
accretion. Note that vigorous PNS convection develops from $\sim 0.05$ s after
bounce and the early postbounce $E_{\rm kin,PNS,\theta}$ is due to the prompt
convection.

###  3.2 Three-dimensional simulation

To explore the potential impact of the imposed axisymmetry on the core-
collapse phase of ECSNe in our 2D simulations, we performed a 3D simulation in
an octant with m1 as the progenitor model. The 3D simulation has a coarser
resolution (see Section 2.2) than the 2D runs in Section 3.1, so we
additionally performed a 2D simulation (m1-2d-lr) with a similar resolution to
disentangle the effects of resolution and dimension.

Figure 6: Same as Fig. 3, but for the magnitude of angular averaged transverse
velocity $\langle|v_{\theta}|\rangle$ as a function of radius and the
postbounce time in model m1-3d-lr. The mass resolution is same as the right
plot in Fig. 3 for model m1, namely 0.1 $M_{\odot}$, 0.01 $M_{\odot}$ and
0.001 $M_{\odot}$ below $M=$1.3, 1.34 and 1.359 $M_{\odot}$, which are
indicated by the blue dashed lines. The white, red and yellow lines mark the
mean shock radius, mean gain radius and PNS surface.

Fig. 6 is the 3D counterpart of the right panel in Fig. 3. It is obvious that
the explosion dynamics and structure of turbulent convection are similar in
the 2D and 3D simulations. Compared to the 2D results, in 3D
$\langle|v_{\theta}|\rangle$ looks smoother in the gain region during 0.1-0.2
s after bounce and the PNS convective region is more extended. The magnitude
of $\langle|v_{\theta}|\rangle$ is comparable between 2D and 3D simulations.

Figure 7: Same as Fig. 2, but for the postbounce evolution of the mean shock
radius, mass of ejecta and explosion energy in models m1 (solid), m1-2d-lr
(dash-dotted) and m1-3d-lr (dotted). Note that the discrepancy is small
between 2D and 3D simulations with a similar resolution. Figure 8: Same as
Fig. 5, but for the postbounce evolution of the mass, radius and transverse
kinetic energy of PNS in models m1 (solid), m1-2d-lr (dash-dotted) and
m1-2d-3d (dotted). Note that at early postbounce time $E_{\rm kin,PNS,\theta}$
includes contribution of prompt convection. The difference between m1-3d-lr
and m1-2d-lr at $t_{\rm pb}=0.05-0.15$ s is due to the more extended
convective regions in 3D (see Fig. 6) which makes it more difficult to
separate the PNS and postshock convection.

As listed in Table 2, when the 3D simulation terminates at 300 ms after
bounce, more mass (by $\sim$20%) is ejected with a larger explosion energy (by
$\sim 30\%$) than the fiducial 2D simulation with the same progenitor.
However, when compared to the low-resolution 2D counterpart (m1-2d-lr), the
numerical results are more compatible, see Fig. 7 for the explosion dynamics
and Fig. 8 for PNS properties. This suggests that the differences seen here
are _not_ due to dimension but result from different resolutions.

Figure 9: Mass (top panel), neutrino heating rate (middle panel) and
transverse kinetic energy (bottom panel) in the gain region as a function of
the postbounce time for model m1 (solid), m1-2d-lr (dash-dotted) and m1-3d-lr
(dotted). Note that for the 3D simulation, the transverse kinetic energy
includes both $\theta$\- and $\phi$-direction motions.

To have a closer look at the difference between 2D and 3D simulations, in Fig.
9 we compare the mass, neutrino heating rate and transverse kinetic energy
$E_{\rm kin,gain,\theta}$ in the gain region for the model m1, m1-2d-lr and
m1-3d-lr. Note that for the 3D simulation, $E_{\rm kin,gain,\theta}$ includes
both $\theta$\- and $\phi$-direction contributions. There is a visible
difference between simulations with the fiducial resolution and low resolution
(2D and 3D), especially after $t_{\rm pb}\approx$0.05 s. The larger mass and
heating rate in the gain region are due to a less compact PNS in low-
resolution runs (see Fig. 8). The difference in $E_{\rm kin,gain,\theta}$ can
result from both different heating rates and angular resolutions
($0.6^{\circ}$ vs. $1.5^{\circ}$).

For the low-resolution simulations m1-2d-lr and m1-3d-lr, the large
discrepancy of $E_{\rm kin,gain,\theta}$ in the early postbounce phase is
likely due to a different manifestation of the violent prompt convection
inside the shock. Afterwards, quantitative agreement is found between the low-
resolution runs in the episode when the explosion energy increases quickly
($\sim 0.1-0.3$ s after bounce). Our results indicate a minor impact of 3D
over 2D simulations for the collapse and explosion phases of ECSN, comparing
to the 25%-30% difference between 2D and 1D simulations (see also Radice et
al., 2017).

###  3.3 Properties of ejecta

Figure 10: Slices of entropy (upper panel) and $Y_{e}$ (lower panel) at 300 ms
postbounce in the 6 axisymmetric simulations. Figure 11: Similar as Fig. 10,
but for the slices of entropy (upper panel) and $Y_{e}$ (lower panel) at 300
ms postbounce in the model m1-2d-lr and m1-3d-lr. Note that for m1-3d-lr, the
$r$-axis stands for the diagonal in the $xy$-plane.

It has been suggested that ECSNe can be a site for $r$-process nucleosynthesis
due to the early shock revival and prompt explosion (Hillebrandt et al., 1984;
Wanajo et al., 2003; Ning et al., 2007). However, recent ECSN simulations with
accurate modelling of neutrino transport (Kitaura et al., 2006; Janka et al.,
2008) have shown that it is unlikely to eject material with a low $Y_{e}$
($\ll 0.4$). A more detailed study showed that ECSN can produce abundant light
trans-Fe elements from Zn to Zr (Wanajo et al., 2018). Here we present the
properties of the ejecta in our simulations and the implications for the ECSN
nucleosynthesis.

In Fig. 10, we show the 2D profiles of entropy and $Y_{e}$ at 300 ms after
bounce in the 6 axisymmetric simulations with our fiducial resolution. The
ejecta contains three components, which are the neutrino-driven PNS wind,
neutrino-heated and early-shocked, from small to large radius. The early-
shocked ejecta has a $Y_{e}$ of $\sim$0.5 and entropy from high ($\sim
35~{}k_{\rm B}~{}{\rm nuc}^{-1}$) to low ($\sim 10~{}k_{\rm B}~{}{\rm
nuc}^{-1}$). The neutrino-heated ejecta has a $Y_{e}$ of $\sim$0.35-0.5 and an
entropy of $\sim 10$–$20~{}k_{\rm B}~{}{\rm nuc}^{-1}$. This component has
inherited the aspherical structure from the neutrino-driven turbulent
convection. The later ejecta in the neutrino-driven PNS wind is proton-rich
with a $Y_{e}$ slightly above 0.5 ($\lesssim 0.52$) and an entropy $\gtrsim
20~{}k_{\rm B}~{}{\rm nuc}^{-1}$. There is a modest and subtle dependence of
the morphology of ejecta on the employed progenitor model, especially for the
neutrino-heated component. In Fig. 11, we show the slices of entropy and
$Y_{e}$ for the low-resolution runs at 300 ms after bounce. Note that for the
3D model, the $r$-axis stands for the diagonal in the $xy$-plane.

Figure 12: Mass distribution of ejecta as a function of the electron fraction
($Y_{e}$) and entropy per baryon ($s$) at 300 ms postbounce in models m1-2d-lr
and m1-3d-lr.

To diagnose the impact of 3D, in Fig. 12 we compare the mass distribution of
ejecta in the $Y_{e}$-entropy space at 300 ms after bounce in m1-2d-lr (left
panel) and m1-3d-lr (right panel). The lack of high-entropy ejecta with
$Y_{e}\approx 0.5$ (vertical strike at $Y_{e}\approx 0.5$) in m1-3d-lr is due
to matter that has left the computational grid, with a small domain in 3D
($8\times 10^{3}$ km) than in 2D ($1.2\times 10^{4}$ km). Otherwise the mass
distribution looks very similar between 2D and 3D simulations and is
consistent with the slice view of the fiducial 2D runs (Fig. 10).

Figure 13: Mass distribution of ejecta as a function of the electron fraction ($Y_{e}$) in all our multidimensional ECSN simulations. The width of binning in $Y_{e}$ is 0.005. Model | $M_{\rm ej}~{}(10^{-3}M_{\odot})$
---|---
| $Y_{e}\leq 0.40$ | $0.4<Y_{e}\leq 0.45$ | $0.45<Y_{e}\leq 0.5$ | $Y_{e}>0.5$
n8.8-1d | 0 | 0 | 11.68 | 1.77
m1-1d | 0 | 0 | 11.59 | 1.46
n8.8 | 0.53 | 5.03 | 8.32 | 2.61
s1 | 0.59 | 4.18 | 9.11 | 2.40
s2 | 0.47 | 4.47 | 9.71 | 2.32
m1 | 0.76 | 4.11 | 9.34 | 2.25
m2 | 0.49 | 4.91 | 9.10 | 1.51
l1 | 0.19 | 5.19 | 9.07 | 1.54
l2 | 1.82 | 5.52 | 7.27 | 1.52
m1-2d-lr | 1.55 | 7.78 | 11.32 | 1.16
m1-3d-lr | 0.92 | 7.75 | 9.88 | 1.87
Table 3: Detailed mass distribution of ejecta as a function of $Y_{e}$ at 300
ms after bounce in all our simulations (cf. Table 2).

In Fig. 13, we plot the mass distribution of ejecta as a function of $Y_{e}$
at 300 ms after bounce for all the simulations. The left panel in Fig. 13
provides an overview for the progenitor dependence of ejecta. Progenitor
models with a large $R_{\rm flame}$ (groups m and l) result in an extended
tail at $Y_{e}\lesssim 0.4$ down to $\sim$0.35. The total mass in this
low-$Y_{e}$ tail depends also on the early shock expansion after bounce so
that it shows a non-monotonic trend with respect to $R_{\rm flame}$. For model
l2 which ejects the most low-$Y_{e}$ material, ejecta with $0.35\leq
Y_{e}<0.4$ has a mass of $\sim 1.8\times 10^{-3}~{}M_{\odot}$ and an entropy
of 10-15 $k_{\rm B}~{}{\rm nuc}^{-1}$, belonging to the neutrino-heated
ejecta. This may contribute to some neutron-rich nucleosynthesis and deserves
further calculation with careful post-process calculation. Matter with
$0.4<Y_{e}<0.52$ comprises the bulk ($\sim 99\%$ for l1 to $\sim 91\%$ for l2)
of ejecta. The 3D result (m1-3d-lr) is consistent with that of m1-2d-lr (right
panel of Fig. 13) except for a little less low-$Y_{e}$ component. The low-
resolution simulations (both in 2D and 3D) yield a larger mass of ejecta with
$0.4\lesssim Y_{e}\lesssim 0.45$, in agreement with their more powerful
neutrino heating. This again shows the minor impact of 3D on ECSN explosions.

Overall, our results are fairly compatible with Wanajo et al. (2011); Wanajo
et al. (2018) though two differences are found. One is the low-$Y_{e}$ ejecta
($Y_{e}<0.4$) due to the 2D progenitor models of Zha et al. (2019) employed in
this work, which have an aspherical flame front and larger $R_{\rm flame}$
compared to model n8.8. The other is that the neutrino-driven wind has a
$Y_{e}$ of $\lesssim 0.52$ in this work while its maximum $Y_{e}$ reaches
$\sim$0.55 in Wanajo et al. (2018). This may be due to the different neutrino-
transport schemes used in this work and Wanajo et al. (2018) that deserves
further exploration.

###  3.4 Discussion

The degenerate ONeMg core of ECSN progenitors possess an outwardly propagating
oxygen flame, which is not captured in our neutrino-transport simulations.
According to Jones et al. (2016) and Leung et al. (2020), the flame speed,
either laminar or turbulent, is less than $\sim 0.1c_{\rm s}\simeq
10^{-3}c=300~{}$km s-1, where $c_{\rm s}$ ($c$) is the speed of sound (light).
Therefore, before the shock breaks out of the core surface which takes $\sim
0.15$ s from the start of simulations, the flame front can propagate outwards
for $\sim 45$ km at most. The propagated distance is small compared with both
the radius of the initial flame front $R_{\rm flame}$ and the difference in
$R_{\rm flame}$ among different progenitor models employed in this work. This
suggests that the ignorance of the propagating oxygen flame should not affect
the properties of ECSN explosion too much.

On the other hand, the deflagration wave might transit to a detonation wave
with a supersonic propagation speed, similar to that in a Type Ia supernova
(Khokhlov, 1991). If such a transition would happen in an ECSN, it might form
a shock that could eject the outer layer. However, the propagation of the
oxygen flame through $R_{\rm flame}$ as simulated in Zha et al. (2019) does
not reach the condition of the deflagration-detonation transition (DDT)
adopted in Leung & Nomoto (2018) and Leung et al. (2020). This is because the
nuclear fuel is oxygen which is more difficult to ignite compared with carbon
and also the density at the oxygen flame front is higher and less turbulent
compared with the Type Ia supernova model, although the condition of the DDT
is highly uncertain (Röpke, 2017). We also note that when including the flame
propagation in a spherically-symmetric ECSN simulation, Takahashi et al.
(2019) found an enhancement of the flame speed to $\sim 1000$ km s-1 due to
neutrino-electron scattering. This new driven mechanism for flame propagation
has not been confirmed elsewhere, and should be further checked in multi-
dimensional neutrino-transport simulations with a proper flame-capturing
scheme.

Our simulations have treated the matter outside the flame front as iron for
densities above $4\times 10^{6}$g cm-3 and ignored nuclear burning at lower
densities. Things are flashed into NSE composition when crossing the
transition density of EOS. The binding energy difference between 56Ni and O-Ne
of 0.1 $M_{\odot}$ is $\sim 1.4\times 10^{50}$ erg, which is comparable to the
ECSN explosion energy found here. Nevertheless, the explosion energy shows a
small difference ($\sim 0.1\times 10^{50}$ erg) among progenitor models with a
range of mass of the NSE ash (0.09 to 0.79 $M_{\odot}$, cf. Table 1). It
suggests that as most ($\sim 98.5\%$) of the ONeMg core is quickly accreted by
the PNS and remains bounded, the detailed composition has a minor effect on
the explosion energy and dynamics. We leave a proper treatment for materials
outside the flame front (e.g., Janka et al., 2008) for a future study.

##  4 Conclusions

In this paper, we have presented 2D axisymmetric ECSN simulations using FLASH
with progenitor model n8.8 of Nomoto (1984, 1987) and six 2D progenitor models
of Zha et al. (2019). We found a weak dependence of the properties of ECSN
explosions on progenitor models. The diagnostic explosion energy at the
termination of simulations (400 ms after bounce) is $\sim 1.36$–$1.48\times
10^{50}$ erg with $\sim 0.017$–$0.018M_{\odot}$ being ejected, in the range of
results reported previously in Janka et al. (2008); Fischer et al. (2010);
Tauris et al. (2017). The PNS has a baryonic mass of 1.34 $M_{\odot}$ with the
new progenitor models, slightly lighter than that with n8.8 ($\sim
1.36M_{\odot}$). This reflects the difference in the ONeMg core mass due to
updated weak-interaction rates. In addition, the extended flame front leads to
more violent convective motions behind the shock, and results in a low-$Y_{e}$
(0.35-0.4) tail in the mass distribution of the ejecta. This low-$Y_{e}$
ejecta has a maximum mass of $\sim 1.8\times 10^{-3}~{}M_{\odot}$ and may
contribute to some weak neutron-rich nucleosynthesis, which deserves further
investigation with post-process calculations.

With the 3D octant simulation and its 2D counterpart with a similar
resolution, we demonstrated that 3D has a small impact on the dynamical
processes in ECSN, while multi-dimensional simulations result in 25%-30%
difference in explosion properties comparing to 1D simulations. The major
phenomenon observed in the 3D run is the smoother and more extended turbulent
structure inside convective zones than that in 2D simulations. We attribute
the small impact of 3D (comparing to the case of CCSN simulations) to the
subdominant role of turbulence in the ECSN explosion. Future studies can
investigate the effect of resolution in 3D ECSN simulations, perhaps in $4\pi$
geometry and with a true 3D progenitor model (Jones et al., 2016).

Finally, we expect that our neglect of the flame propagation after the onset
of the ONeMg core collapse will not significantly alter our conclusions.
Nevertheless, rigorous modelling of the nuclear burning during the collapse
and explosion is needed to yield more accurate numerical results for the
explosion energy, the mass and composition of the ejecta, and the PNS mass,
which are the necessary inputs for the predictions of electromagnetic signals
and nucleosynthesis of ECSN.

## Acknowledgements

This work is supported by the Swedish Research Council (Project No.
2020-00452). The simulations were enabled by resources provided by the Swedish
National Infrastructure for Computing (SNIC) at PDC and NSC partially funded
by the Swedish Research Council through grant agreement No. 2016-07213. S.C.L.
acknowledges support by NASA grants HST-946AR-15021.001-A and 80NSSC18K1017.
K.N. is supported by the World Premier International Research Center
Initiative (WPI), MEXT, Japan, and the Japan Society for the Promotion of
Science (JSPS) KAKENHI grant JP17K05382, JP20K04024, and JP21H04499.

S.M.C. is supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, Early Career Research Program under Award Number
DE-SC0015904. This material is based upon work supported by the U.S.
Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research and Office of Nuclear Physics, Scientific Discovery through
Advanced Computing (SciDAC) program under Award Number DE- SC0017955. This
work was supported by the Exascale Computing Project (17-SC-20-SC), a
collaborative effort of the U.S. Department of Energy Office of Science and
the National Nuclear Security Administration.

## Data Availability

The ECSN progenitor models used in this work are publicly available at
https://doi.org/10.5281/zenodo.5748457. Results of the hydrodynamic
simulations will be shared on reasonable request to the corresponding author.

## References

  * Botticella et al. (2009) Botticella M. T., et al., 2009, MNRAS, 398, 1041
  * Bruenn & Mezzacappa (1994) Bruenn S. W., Mezzacappa A., 1994, ApJ, 433, L45
  * Burrows & Vartanyan (2021) Burrows A., Vartanyan D., 2021, Nature, 589, 29
  * Cai et al. (2021) Cai Y. Z., et al., 2021, arXiv e-prints, p. arXiv:2108.05087
  * Cardall et al. (2013) Cardall C. Y., Endeve E., Mezzacappa A., 2013, Phys. Rev. D, 87, 103004
  * Couch & Ott (2013) Couch S. M., Ott C. D., 2013, ApJ, 778, L7
  * Denissenkov et al. (2015) Denissenkov P. A., Truran J. W., Herwig F., Jones S., Paxton B., Nomoto K., Suzuki T., Toki H., 2015, MNRAS, 447, 2696
  * Eggenberger Andersen et al. (2021) Eggenberger Andersen O., Zha S., da Silva Schneider A., Betranhandy A., Couch S. M., O’Connor E. P., 2021, arXiv e-prints, p. arXiv:2106.09734
  * Farrow et al. (2019) Farrow N., Zhu X.-J., Thrane E., 2019, ApJ, 876, 18
  * Fields & Couch (2020) Fields C. E., Couch S. M., 2020, ApJ, 901, 33
  * Fields & Couch (2021) Fields C. E., Couch S. M., 2021, ApJ, 921, 28
  * Fischer et al. (2010) Fischer T., Whitehouse S. C., Mezzacappa A., Thielemann F. K., Liebendörfer M., 2010, A&A, 517, A80
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Gessner & Janka (2018) Gessner A., Janka H.-T., 2018, ApJ, 865, 61
  * Hillebrandt et al. (1984) Hillebrandt W., Nomoto K., Wolff R. G., 1984, A&A, 133, 175
  * Hiramatsu et al. (2021) Hiramatsu D., et al., 2021, Nature Astronomy, 5, 903
  * Janka et al. (2008) Janka H. T., Müller B., Kitaura F. S., Buras R., 2008, A&A, 485, 199
  * Jones et al. (2016) Jones S., Röpke F. K., Pakmor R., Seitenzahl I. R., Ohlmann S. T., Edelmann P. V. F., 2016, A&A, 593, A72
  * Khokhlov (1991) Khokhlov A. M., 1991, A&A, 246, 383
  * Kippenhahn et al. (2012) Kippenhahn R., Weigert A., Weiss A., 2012, Stellar Structure and Evolution, doi:10.1007/978-3-642-30304-3.
  * Kirsebom et al. (2019) Kirsebom O. S., et al., 2019, Phys. Rev. Lett., 123, 262701
  * Kitaura et al. (2006) Kitaura F. S., Janka H.-T., Hillebrandt W., 2006, A&A, 450, 345
  * Knigge et al. (2011) Knigge C., Coe M. J., Podsiadlowski P., 2011, Nature, 479, 372
  * Kozyreva et al. (2021) Kozyreva A., Baklanov P., Jones S., Stockinger G., Janka H.-T., 2021, MNRAS, 503, 797
  * Kruckow et al. (2018) Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G., 2018, MNRAS, 481, 1908
  * Langanke et al. (2021) Langanke K., Martínez-Pinedo G., Zegers R. G. T., 2021, Reports on Progress in Physics, 84, 066301
  * Leung & Nomoto (2018) Leung S.-C., Nomoto K., 2018, ApJ, 861, 143
  * Leung & Nomoto (2019) Leung S.-C., Nomoto K., 2019, Publications of the Astronomical Society of Australia, 36, e006
  * Leung et al. (2015) Leung S. C., Chu M. C., Lin L. M., 2015, MNRAS, 454, 1238
  * Leung et al. (2020) Leung S.-C., Nomoto K., Suzuki T., 2020, ApJ, 889, 34
  * Martínez-Pinedo et al. (2014) Martínez-Pinedo G., Lam Y. H., Langanke K., Zegers R. G. T., Sullivan C., 2014, Phys. Rev. C, 89, 045806
  * Miyaji & Nomoto (1987) Miyaji S., Nomoto K., 1987, ApJ, 318, 307
  * Miyaji et al. (1980) Miyaji S., Nomoto K., Yokoi K., Sugimoto D., 1980, Publications of the Astronomical Society of Japan, 32, 303
  * Moriya et al. (2014) Moriya T. J., Tominaga N., Langer N., Nomoto K., Blinnikov S. I., Sorokina E. I., 2014, A&A, 569, A57
  * Müller et al. (2017) Müller B., Melson T., Heger A., Janka H.-T., 2017, MNRAS, 472, 491
  * Nagakura et al. (2020) Nagakura H., Burrows A., Radice D., Vartanyan D., 2020, MNRAS, 492, 5764
  * Ning et al. (2007) Ning H., Qian Y. Z., Meyer B. S., 2007, ApJ, 667, L159
  * Nomoto (1984) Nomoto K., 1984, ApJ, 277, 791
  * Nomoto (1987) Nomoto K., 1987, ApJ, 322, 206
  * Nomoto & Hashimoto (1988) Nomoto K., Hashimoto M., 1988, Phys. Rep., 163, 13
  * Nomoto & Kondo (1991) Nomoto K., Kondo Y., 1991, ApJ, 367, L19
  * Nomoto & Leung (2017) Nomoto K., Leung S.-C., 2017, Handbook of Supernovae, eds. Alsabti, A.W. & Murdin, P. (Springer), 2017. p. p. 483, doi:10.1007/978-3-319-21846-5_118
  * O’Connor (2015) O’Connor E., 2015, ApJS, 219, 24
  * O’Connor & Couch (2018a) O’Connor E. P., Couch S. M., 2018a, ApJ, 854, 63
  * O’Connor & Couch (2018b) O’Connor E. P., Couch S. M., 2018b, ApJ, 865, 81
  * O’Connor et al. (2018) O’Connor E., et al., 2018, Journal of Physics G Nuclear Physics, 45, 104001
  * Paczyński (1973) Paczyński B., 1973, Astrophys. Lett., 15, 147
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Podsiadlowski et al. (2004) Podsiadlowski P., Langer N., Poelarends A. J. T., Rappaport S., Heger A., Pfahl E., 2004, ApJ, 612, 1044
  * Radice (2022) Radice D., 2022, private communication
  * Radice et al. (2017) Radice D., Burrows A., Vartanyan D., Skinner M. A., Dolence J. C., 2017, ApJ, 850, 43
  * Röpke (2017) Röpke F. K., 2017, Combustion in Thermonuclear Supernova Explosions. p. 1185, doi:10.1007/978-3-319-21846-5_58
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Schwab et al. (2010) Schwab J., Podsiadlowski P., Rappaport S., 2010, ApJ, 719, 722
  * Schwab et al. (2017) Schwab J., Martínez-Rodríguez H., Piro A. L., Badenes C., 2017, ApJ, 851, 105
  * Schwab et al. (2020) Schwab J., Farmer R., Timmes F. X., 2020, ApJ, 891, 5
  * Shibata et al. (2011) Shibata M., Kiuchi K., Sekiguchi Y., Suwa Y., 2011, Progress of Theoretical Physics, 125, 1255
  * Steiner et al. (2013) Steiner A. W., Hempel M., Fischer T., 2013, ApJ, 774, 17
  * Stockinger et al. (2020) Stockinger G., et al., 2020, MNRAS, 496, 2039
  * Strömberg et al. (2021) Strömberg D. F., Martínez-Pinedo G., Nowacki F., 2021, arXiv e-prints, p. arXiv:2104.02614
  * Suzuki et al. (2016) Suzuki T., Toki H., Nomoto K., 2016, ApJ, 817, 163
  * Suzuki et al. (2019) Suzuki T., Zha S., Leung S.-C., Nomoto K., 2019, ApJ, 881, 64
  * Takahashi et al. (2019) Takahashi K., Sumiyoshi K., Yamada S., Umeda H., Yoshida T., 2019, ApJ, 871, 153
  * Tauris et al. (2017) Tauris T. M., et al., 2017, ApJ, 846, 170
  * Timmes & Swesty (2000) Timmes F. X., Swesty F. D., 2000, ApJS, 126, 501
  * Timmes & Woosley (1992) Timmes F. X., Woosley S. E., 1992, ApJ, 396, 649
  * Toki et al. (2013) Toki H., Suzuki T., Nomoto K., Jones S., Hirschi R., 2013, Phys. Rev. C, 88, 015806
  * Tolstov et al. (2019) Tolstov A., Nomoto K., Sorokina E., Blinnikov S., Tominaga N., Taniguchi Y., 2019, The Astrophysical Journal, 881, 35
  * Tominaga et al. (2013) Tominaga N., Blinnikov S. I., Nomoto K., 2013, ApJ, 771, L12
  * Vartanyan et al. (2018) Vartanyan D., Burrows A., Radice D., Skinner M. A., Dolence J., 2018, MNRAS, 477, 3091
  * Vartanyan et al. (2022) Vartanyan D., Coleman M. S. B., Burrows A., 2022, MNRAS, 510, 4689
  * Wanajo et al. (2003) Wanajo S., Tamamura M., Itoh N., Nomoto K., Ishimaru Y., Beers T. C., Nozawa S., 2003, ApJ, 593, 968
  * Wanajo et al. (2011) Wanajo S., Janka H.-T., Müller B., 2011, ApJ, 726, L15
  * Wanajo et al. (2018) Wanajo S., Müller B., Janka H.-T., Heger A., 2018, ApJ, 852, 40
  * Witt et al. (2021) Witt M., et al., 2021, ApJ, 921, 19
  * Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
  * Zha et al. (2019) Zha S., Leung S.-C., Suzuki T., Nomoto K., 2019, ApJ, 886, 22