# Sequential Inference of Hospitalization Electronic Health Records Using
Probabilistic Models

Alan D. Kaplan and Priyadip Ray  Lawrence Livermore National Laboratory
Livermore, USA  John D. Greene and Vincent X. Liu  Kaiser Permanente Division
of Research
Oakland, USA

###### Abstract

In the dynamic hospital setting, decision support can be a valuable tool for
improving patient outcomes. Data-driven inference of future outcomes is
challenging in this dynamic setting, where long sequences such as laboratory
tests and medications are updated frequently. This is due in part to
heterogeneity of data types and mixed-sequence types contained in variable
length sequences. In this work we design a probabilistic unsupervised model
for multiple arbitrary-length sequences contained in hospitalization
Electronic Health Record (EHR) data. The model uses a latent variable
structure and captures complex relationships between medications, diagnoses,
laboratory tests, neurological assessments, and medications. It can be trained
on original data, without requiring any lossy transformations or time binning.
Inference algorithms are derived that use partial data to infer properties of
the complete sequences, including their length and presence of specific
values. We train this model on data from subjects receiving medical care in
the Kaiser Permanente Northern California integrated healthcare delivery
system. The results are evaluated against held-out data for predicting the
length of sequences and presence of Intensive Care Unit (ICU) in
hospitalization bed sequences. Our method outperforms a baseline approach,
showing that in these experiments the trained model captures information in
the sequences that is informative of their future values.

###### Index Terms:

Electronic Health Records, Machine Learning, Inference, Probabilistic
Modeling, Latent Variable Methods

LLNL-TR-849343

##  I Introduction

Many different machine learning approaches have been explored in the space of
EHR modeling and prediction (see, e.g. [1, 2]). Methods have been developed
for use in disease prediction [3, 4, 5], length of stay prediction [6, 7],
biomarker discovery [8, 9], treatment planning [10, 11], and subgroup analysis
[12, 13, 14, 15]. These approaches also vary with the type of input data used,
which can include summary statistics, images, time series, and sequences.

For this work, we construct a probabilistic model of sequences contained in
episodic EHR data collected during hospitalization. Probabilistic methods
allow us to compute distributions of quantities of interest, which is useful
in an uncertain clinical environment. In addition, one model can be trained
over all of the data and used for multiple inference tasks by computing
conditional probabilities of interest. Modeling of multiple sequences, such as
medications and laboratory tests, has several applications, including
prognosis of patient health trajectories, resource allocation, treatment
planning, and outcome prediction. A single model that can perform multiple
tasks is beneficial in a dynamic clinical environment, where both the nature
of the known data and desired predictive outcomes may be changing often.

Our goal for this work is to develop a flexible approach that can predict
properties of future values of a patient’s EHR given partial sequences of
their EHR, using a single probabilistic model. This requires the construction
of the underlying probability model and derivation of inference algorithms
that use this model.

We use the formulation of mixture models, which have been previously
investigated in the context of EHR data [16, 17, 18, 15]. Both probabilistic
and neural network based approaches have been used that incorporate sequential
dynamics of EHR data [19, 20, 21, 22, 15, 23, 24]. In addition, combinations
of probabilistic subgrouping methods and dynamic models have been developed
[25, 26, 27, 28]. In contrast with this work, these methods do not incorporate
sequences with simultaneous observations (such as multiple medications
administered at one timepoint). Also, these prior works do not derive general
inference methods that can be used to estimate arbitrary future values of EHR
sequences. Prior work that developed mixture models for EHR data showed that
learned subgroup phenotypes can be used to gain insight from sequences in EHR
data [15]. The model in the present work utilizes a mixture model at the top
level, but differs in its latent variable structure. In addition, in this work
we focus on inference, by deriving expressions for distributions of future
values and perform validation of prediction.

The model developed in this work is defined by a full joint probability
distribution over the data components and can fit the data without requiring
any lossy transformations. Inference procedures are derived for estimating
final sequence lengths at discharge and future values using conditional
likelihoods from the original model. Once trained, we evaluate the performance
of the method on data from Kaiser Permanente Northern California (KPNC). We
perform these tasks at both the individual and population-levels and show that
they improve performance beyond baseline population statistic approaches.

##  II Data

KPNC is a highly integrated healthcare delivery system with 21 medical centers
caring for an overall population of 4 million members. For this work, we use a
dataset consisting of 244,248 inpatient hospitalization visits with a
suspected or confirmed infection and sepsis diagnosis, drawn from KPNC medical
centers between 2009 and 2013 [29].

We consider the following data elements:

  * •

Age: $\phi_{a}$

  * •

Sex: $\phi_{s}$

  * •

Death: $\phi_{d}$

  * •

Beds: $\boldsymbol{\beta}=\\{\beta_{1},\ldots,\beta_{|\boldsymbol{\beta}|}\\}$

  * •

Admission Diagnoses:
$\boldsymbol{\alpha}=\\{\alpha_{1},\ldots,\alpha_{|\boldsymbol{\alpha}|}\\}$

  * •

Discharge Diagnoses:
$\boldsymbol{\delta}=\\{\delta_{1},\ldots,\delta_{|\boldsymbol{\delta}|}\\}$

  * •

Laboratory Tests:
$\boldsymbol{\lambda}=\\{\boldsymbol{\lambda}_{1},\ldots,\boldsymbol{\lambda}_{|\boldsymbol{\lambda}|}\\}=\\{\\{\lambda_{1,1},\ldots,\lambda_{1,|\boldsymbol{\lambda}_{1}|}\\},\ldots,\\{\lambda_{|\boldsymbol{\lambda}|,1},\ldots,\lambda_{|\boldsymbol{\lambda}|,|\boldsymbol{\lambda}_{|\boldsymbol{\lambda}|}|}\\}\\}$

  * •

Neurological Tests:
$\boldsymbol{\nu}=\\{\boldsymbol{\nu}_{1},\ldots,\boldsymbol{\nu}_{|\boldsymbol{\nu}|}\\}=\\{\\{\nu_{1,1},\ldots,\nu_{1,|\boldsymbol{\nu}_{1}|}\\},\ldots,\\{\nu_{|\boldsymbol{\nu}|,1},\ldots,\nu_{|\boldsymbol{\nu}|,|\boldsymbol{\nu}_{|\boldsymbol{\nu}|}|}\\}\\}$

  * •

Medications:
$\boldsymbol{\mu}=\\{\boldsymbol{\mu}_{1},\ldots,\boldsymbol{\mu}_{|\boldsymbol{\mu}|}\\}=\\{\\{\mu_{1,1},\ldots,\mu_{1,|\boldsymbol{\mu}_{1}|}\\},\ldots,\\{\mu_{|\boldsymbol{\mu}|,1},\ldots,\mu_{|\boldsymbol{\mu}|,|\boldsymbol{\mu}_{|\boldsymbol{\mu}|}|}\\}\\}$

Bold symbols indicate vectors of values. In cases where we refer to any of
these elements, we use the symbol $\boldsymbol{x}$. The number of values in
$\boldsymbol{x}$ is $|\boldsymbol{x}|$. Admission Diagnoses and Discharge
Diagnoses are collections of ICD-9 codes that are given apon admission and
discharge, respectively. Sequential information that are collected throughout
the each episode include Beds, Laboratory Tests, Neurological Tests, and
Medications. These sequences are of varying length and may contain multiple
observations of varying quantity at each timepoint. For these sequences, the
time associated with a value is denoted $t_{x,i}$. For example, the timestamp
(in hours) of the collection of laboratory tests $\boldsymbol{\lambda_{i}}$ is
$t_{\lambda,i}$ and the timestamp for the medications $\boldsymbol{\mu_{i}}$
is $t_{\mu,i}$. The collection of data for one episode is:
$\boldsymbol{y}=\left(\phi_{a},\phi_{s},\phi_{d},\boldsymbol{\beta},\boldsymbol{\alpha},\boldsymbol{\delta},\boldsymbol{\lambda},\boldsymbol{\nu},\boldsymbol{\mu}\right)$.

##  III Methods

###  III-A Model structure

In this section we describe the probabilistic model that computes the
likelihood of an episode, $f\left(\boldsymbol{y}\right)$. Figure 1 shows the
structure of the model. Colored circles represents single values from the data
elements and arrows indicate probabilistic dependencies.

Figure 1: This figure shows a graphical representation of the model. Arrows
indicate dependencies incorporated the model. Each colored circle represents a
single value of data.

The model utilizes a mixture model formulation with a latent state $z$ that
characterizes associations between data elements. The data elements of the
episode are independent conditioned on the latent variable, so that

$f\left(\boldsymbol{y}|Z=z\right)=f\left(\phi_{a}|Z=z\right)f\left(\phi_{s}|Z=z\right)\cdots
f\left(\boldsymbol{\mu}|Z=z\right),$ (1)

and the distribution of the model is formed by summing over the latent
variable,

$f\left(\boldsymbol{y}\right)=\sum_{z}p_{z}f\left(\boldsymbol{y}|Z=z\right).$

Each product term in Eq 1 is specified with distribution families that match
the corresponding data types. For Age, we use a quantized and truncated
Gaussian distribution,

$f\left(\phi_{a}|Z=z\right)=\frac{1}{C\sqrt{2\pi\sigma_{z}^{2}}}e^{\frac{\left(m_{z}-\phi_{a}\right)^{2}}{2\sigma_{z}^{2}}},$

where the normalizing value ensures that the distribution sums to 1,
$C=\sum_{\phi_{a}}\frac{1}{\sqrt{2\pi\sigma_{z}^{2}}}e^{\frac{\left(m_{z}-\phi_{a}\right)^{2}}{2\sigma_{z}^{2}}}$.

For Sex and Death, we we use Bernoulli distributions,

$f\left(\phi_{s}|Z=z\right)=p_{s,z}^{1-\phi_{s}}\left(1-p_{s,z}\right)^{\phi_{s}},$
$f\left(\phi_{d}|Z=z\right)=p_{d,z}^{1-\phi_{d}}\left(1-p_{d,z}\right)^{\phi_{d}}.$

For Admission Diagnoses, Discharge Diagnoses, Laboratory Tests, and
Neurological Tests, we use products of categorical (i.e. sometimes referred to
as “bag-of-words”) models. The length of each of these sequences for the
entire episode (e.g. $|\boldsymbol{\alpha}|$ for Admission Diagnoses and
$\sum_{i=1}^{|\boldsymbol{\lambda}|}|\boldsymbol{\lambda}_{i}|$ for Laboratory
Tests) is modeled using a Poisson distribution. For these sequences we have,

$f\left(\boldsymbol{\alpha}|Z=z\right)=\frac{l_{\alpha,z}^{|\boldsymbol{\alpha}|}}{|\boldsymbol{\alpha}|!}e^{-l_{\alpha,z}}\prod_{i=1}^{|\boldsymbol{\alpha}|}p_{\alpha,z,\alpha_{i}},$
$f\left(\boldsymbol{\delta}|Z=z\right)=\frac{l_{\delta,z}^{|\boldsymbol{\delta}|}}{|\boldsymbol{\delta}|!}e^{-l_{\delta,z}}\prod_{i=1}^{|\boldsymbol{\delta}|}p_{\delta,z,\delta_{i}},$
$f\left(\boldsymbol{\lambda}|Z=z\right)=\frac{l_{\lambda,z}^{\sum_{i=1}^{|\boldsymbol{\lambda}|}|\boldsymbol{\lambda}_{i}|}}{\sum_{i=1}^{|\boldsymbol{\lambda}|}|\boldsymbol{\lambda}_{i}|!}e^{-l_{\lambda,z}}\prod_{i=1}^{|\boldsymbol{\lambda}|}\prod_{j=1}^{|\boldsymbol{\lambda}_{i}|}p_{\lambda,z,\lambda_{i,j}},$
$f\left(\boldsymbol{\nu}|Z=z\right)=\frac{l_{\nu,z}^{\sum_{i=1}^{|\boldsymbol{\nu}|}|\boldsymbol{\nu}_{i}|}}{\sum_{i=1}^{|\boldsymbol{\nu}|}|\boldsymbol{\nu}_{i}|!}e^{-l_{\nu,z}}\prod_{i=1}^{|\boldsymbol{\nu}|}\prod_{j=1}^{|\boldsymbol{\nu}_{i}|}p_{\nu,z,\nu_{i,j}},$

where $p_{x,z,i}$ is the probability of item $i$ in sequence $x$ for latent
state $z$.

For the Beds sequence we use a Markov Chain to model transitions between
timepoints, so that the order of observations is significant. A Poisson
distribution characters the length of the sequence, leading to,

$f\left(\boldsymbol{\beta}|Z=z\right)=\frac{l_{\beta,z}^{|\boldsymbol{\beta}|}}{|\boldsymbol{\beta}|!}e^{-l_{\beta,z}}p_{\beta,z,\beta_{1}}\prod_{i=1}^{|\boldsymbol{\beta}|-1}q_{\beta,z,\beta_{i},\beta_{i+1}},$

where $p_{\beta,z,i}$ is the probability that the first item is $i$
conditioned on $Z=z$, and $q_{\beta,z,i,j}$ is the probability of
transitioning from item $i$ to $j$ conditioned on $Z=z$.

We also capture the order and transitions between timepoints in the
Medications sequence. However, unlike the Beds sequence, the Medications can
contain multiple observations at any single timepoint. We use a sequence of
latent states, $\boldsymbol{S}$, where each of these states carries a
conditional distribution over a collection of medications of arbitrary length.
In this way we are able to model the ordered sequence and express multiple
simultaneous observations using a Hidden Markov Model with variable length
observations for each state, $s$, characterized by Poisson distributions. The
distribution is,

$f\left(\boldsymbol{\mu}|Z=z\right)=\frac{l_{\mu,z}^{|\boldsymbol{\mu}|}}{|\boldsymbol{\mu}|!}e^{-l_{\mu,z}}\sum_{\boldsymbol{S}}\left(p_{\mu,z,s_{1}}\prod_{i=1}^{|\boldsymbol{\mu}|-1}q_{\mu,z,s_{i},s_{i+1}}\right)\\\
\left(\prod_{i=1}^{|\boldsymbol{\mu}|}\frac{l_{\mu,z,s_{i}}^{|\boldsymbol{\mu}_{i}|}}{|\boldsymbol{\mu}_{i}|!}e^{-l_{\mu,z,s_{i}}}\prod_{j=1}^{|\boldsymbol{\mu}_{i}|}p_{\mu,z,s_{i},\mu_{i,j}}\right),$

where $p_{\mu,z,i}$ is the probability that state $i$ is the initial state,
$q_{\mu,z,i,j}$ is the transition probability from state $i$ to state $j$, and
$p_{\mu,z,i,j}$ is the probability of medication $j$ conditioned on states $i$
and $z$.

###  III-B Estimation

The estimation of the model parameters follows the standard Expectation
Maximization (EM) procedure (see e.g. [30]). The hyperparameters of the model
are the number of top-level latent states $|Z|$ and the number of HMM states
for the Medications sequence, $\mathcal{C}_{S}$. In order to choose these 2
values, we compute the Bayesian Information Criterion (BIC), which penalizes
the model fit by a function of the number of parameters:
$BIC\left(d\right)=d\ln N-2\ln f\left(\boldsymbol{y}\right)$, where $d$ is the
total parameter count and $N$ is the number of episodes. We perform a 2D grid
search over $|Z|$ and $\mathcal{C}_{S}$ and compute the BIC to select these
hyperparameters.

###  III-C Inference

Given a trained model and input data, we can derive inference algorithms to
compute the likelihood of an event of interest. This can include estimation of
the final sequence length (e.g., length of Beds or Laboratory Tests) or
prediction of specific future observations.

Of particular interest, is the inference of future events given a partial set
of EHR data. The collection of data for an episode up until time $t$ is

$\boldsymbol{y}_{[1\ldots
t]}=\left(\phi_{a},\phi_{s},\phi_{d},\boldsymbol{\beta}_{[1\ldots
t]},\boldsymbol{\alpha},\boldsymbol{\lambda}_{[1\ldots
t]},\boldsymbol{\nu}_{[1\ldots t]},\boldsymbol{\mu}_{[1\ldots t]}\right),$

where $\boldsymbol{x}_{[1\ldots
t]}=\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{w}\\}$ for the largest $w$
such that $t_{x,w}<t$ for component $x$. We do not include Discharge Diagnoses
or Death as inputs, since they are unknown until the end of the episode.

The central quantity that is needed to compute is the conditional distribution
of the latent variable $Z$ given the partial data:
$f\left(Z=z|\boldsymbol{y}_{[1\ldots t]}\right)$. This distribution is the
partial sequence likelihood weighted by cumulative Poisson probabilities,

$f\left(Z=z|\boldsymbol{y}_{[1\ldots
t]}\right)=p_{z}f\left(\boldsymbol{y}_{[1\ldots
t]}|Z=z\right)\prod_{x}e^{-l_{x,z}}\sum_{k=t}^{\infty}\frac{l_{x,z}^{k}}{k!}.$

In this expression, we evaluate the likelihood of the partial data for each
$z$ and multiply by the remaining cumulative probability of future sequence
lengths. This accounts for the lengths of the input sequences that have
different Poisson parameters for each latent state, and can be approximated
easily as the values in the Poisson distribution approach zero.

Once this distribution over $Z$ is computed using the partial input sequence
$\boldsymbol{y}_{[1\ldots t]}$, we can infer quantities of interest. To infer
the length of sequence $x$ we compute the probability of the sequence length
given the partial sequence data,

$f\left(k_{x}|\boldsymbol{y}_{[1\ldots
t]}\right)=\sum_{z}f\left(Z=z|\boldsymbol{y}_{[1\ldots
t]}\right)\frac{l_{x,z}^{k}}{k!}e^{-l_{x,z}}.$

To infer future values, the latent state probabilities are used as
coefficients to mix over component distributions. In general, inferring the
distribution of component $\boldsymbol{x}$ can be accomplished by,

$f\left(\boldsymbol{x}|\boldsymbol{y}_{[1\ldots
t]}\right)=\sum_{z}f\left(Z=z|\boldsymbol{y}_{[1\ldots
t]}\right)f\left(\boldsymbol{x}|Z=z\right).$

####  III-C1 Individualized Sequence Length Prediction

Individualized point predictions for sequence lengths may be computed by
finding the mode of the resulting distributions,

$\hat{k}_{x,t}=\arg\max_{k}f\left(k|\boldsymbol{y}_{[1\ldots t]}\right).$

####  III-C2 Population Level Inference

Considering a sequence $\boldsymbol{x}$ and using data up until time $t$, the
expected value of the length of this sequence for subject $i$ is

$E_{x,i,t}=\sum_{k=0}^{\infty}kf\left(k|\boldsymbol{y}_{i,[1\ldots
t]}\right).$

The total population level sequence length estimate is computed by summing
across episodes,

$E_{x,t}=\sum_{i}E_{x,i,t}.$

This is the expected value of the total length across all subjects.

We are also interested in inferring the presence of a specific item in a
sequence. In this work, we will infer the presence of future ICU in the Bed
sequence. In order to do this, we first infer the Bed sequence distribution
using episode data up until time $t$,

$f\left(\boldsymbol{\beta}\right)=\sum_{z}f\left(Z=z|\boldsymbol{y}_{[1\ldots
t]}\right)f\left(\boldsymbol{\beta}|Z=z\right).$

Under this model, the likelihood that ICU will not occur from time $t$ to time
$t+s$ is

$\sum_{z}\frac{l_{\beta,z}^{t+s}}{\left(t+s\right)!}\prod_{t^{\prime}=1}^{s}q_{t+t^{\prime},z},$

where $q_{t,z}$ is the probability that the sequence at time $t$ is not ICU
given $Z=z$,

$q_{t,z}=\frac{l_{\beta,z}^{t}}{t!}e^{-l_{\beta,z}}p_{\beta,z,\beta_{t-1}}\sum_{i\neq
ICU}q_{\beta,z,\beta_{t-1},i}.$

Under all future times, we can compute,

$g=\sum_{z}\sum_{s=1}^{\infty}\frac{l_{\beta,z}^{t+s}}{\left(t+s\right)!}\prod_{t^{\prime}=1}^{s}q_{t+t^{\prime},z},$

which is the probability that ICU will not occur in the future. Our
probability of interest is the ICU will occur, which is $1-g$. This can be
viewed as a Bernoulli trial in which the probability of success (ICU) changes
at every step. The population inference is computed by summing this
probability over the episodes.

##  IV Results

From our dataset of 244,248 episodes, we randomly selected 80% for training.
The two hyperparameters, $|Z|$ and $\mathcal{C}_{S}$ were selected by
performing a linear search and computing the BIC on the training set as
described in Section III. BIC values were computed on a grid consisting of
$|Z|$ taking values 1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100; and
$\mathcal{C}_{S}$ taking values 1, 5, 10, 15. This resulted in 55 trained
models. The model order $|Z|=30$ and $\mathcal{C}_{S}=5$ attained the minimum
BIC value over these models.

###  IV-A Individualized Inference of Sequence Length

Estimated sequence lengths are computed at every available timepoint for each
test episode. Our method is compared to two baseline approaches. The first,
“Constant”, uses the training population average sequence lengths, which is
2.68 for Beds, 197 for Laboratory Tests, and 15.2 for Medications per episode.
The second, “Hold Last”, outputs the current sequence length at every
timepoint (i.e., estimates that the final sequence length is equal to the
present sequence length).

Given the true value for the length of component $x$, $k_{x}^{\star}$, we
quantify the Absolute Error (AE),

$e_{x,t}=|\hat{k}_{x,t}-k_{x}^{\star}|,$

where $\hat{k}_{x,t}$ is the prediction at time $t$ and $k_{x}^{\star}$ is the
final true value.

The results in Table I show the AE per episode for the three methods. Values
in parenthesis are the percentage improvement over the “Constant” method
performance.

Method | Beds | Laboratory Tests | Medications
---|---|---|---
Constant | 0.93 (0) | 129 (0) | 6.9 (0)
Hold Last | 0.72 (23) | 85 (34) | 6.2 (10)
Model | 0.45 (52) | 73 (43) | 4.3 (38)
TABLE I: Absolute error of sequence length prediction results. Values are the
average absolute error per episode and values in parenthesis are percentage
improvement over the Constant baseline approach.

###  IV-B Population-Level Inference

For population-level inference we compute the average absolute error,

$e=100\frac{|y_{pred}\left(t\right)-y\left(t\right)|}{\sum_{t}v\left(t\right)},$

where $y_{pred}\left(t\right)$ is the prediction at time $t$,
$y\left(t\right)$ is the true value at time $t$, and $v\left(t\right)$ is the
number of episodes available at time $t$. The values $y_{pred}$,
$y\left(t\right)$, and $v\left(t\right)$ are accumulated over all episodes in
the test set.

####  IV-B1 Sequence Length Inference

For this problem, we compute the population level inference for sequence
length as given in Section III. We compare our method with the baseline
population rate for each sequence, as determined from the training data. Table
II shows the results.

Method | Beds | Laboratory Tests | Medications
---|---|---|---
Baseline | 3.27 | 364.06 | 22.33
Model | 0.64 | 214.84 | 8.43
TABLE II: Average absolute error for the baseline and model approaches in
inferring the length of each sequence.

####  IV-B2 ICU Presence Inference

For ICU inference, we compare our approach to the population baseline rate of
ICU admittance in the training set. For reference, we also compute the error
for the two edge cases of strictly predicting that ICU will and will not
exist. We call these methods “All” and “None”, respectively. Table III shows
the averaged absolute error as defined above in this Section.

Method | Average AE
---|---
All | 60.86
None | 39.14
Baseline | 15.40
Model | 12.32
TABLE III: Average absolute error for ICU inference.

##  V Discussion and Conclusions

In this work, we defined a model for episodic EHR data containing mixed
sequences and static information. The model is a mixture over probability
distributions tailored to each data type, including collections and sequences.

Expressions were derived that enable inference of future values given partial
input sequence data. These are based on inferring the underlying latent
variable of the mixture model. Inference depends not only on the values, but
also on the lengths of the input sequences.

We trained the model on data from KPNC. Sequence length prediction and
presence of future ICU in the Beds sequence was performed. We find that our
approach outperforms the population average baseline, indicating that the
model is capturing individualized information and is capable of generalizing
beyond the training set.

In Table I we see that the model outperforms the baseline methods for all
three of the sequences. The Hold Last approach produces lower AE than the
Constant approach. This is because there is substantial spread in the sequence
lengths across individuals and the population average does not provide an
accurate estimate. These results show that our method is leveraging
individualized information contained within the sequences to predict sequence
length.

Table II shows that at the population level, our model also outperforms the
baseline. In this approach, we are estimating the total sequence length across
all subjects the test set. The sequence length of the Laboratory Tests were
more difficult to predict than the Beds or Medications. Since these sequences
are very long, it may be that the dynamics are more difficult to learn
compared to other sequences. Prediction of sequence length may be a useful
task for resource planning. Being able to predict utilization of resources
using individualized information rather than population averages can lead to
more accurate estimates of short-term future resource needs.

Prediction of future ICU presence (Table III) shows that our model produces
lower error than the baseline methods. The Baseline method is significantly
better than the edge cases of assuming all patients will be in the ICU (All)
and that no patients will be in the ICU (None). This indicates that the
average rate of ICU presence may be a relatively stable value. There was an
decrease in error, however, when using our model, showing that the model is
capturing patterns that are informative for predicting ICU presence. This
problem could also be a test case of potential use in resource planning, where
more specific information may be needed in addition to the total sequence
length.

Although the inference algorithms are computationally fast to compute, the
training is expensive. Being able to update the model given new training data
may be an important feature to develop for practical adoption of this method.
Various techniques exist for online training of latent variable models, and
the exploration of those methods for this problem is left as future work.

## References

  * [1] J. Wu, J. Roy, and W. F. Stewart, “Prediction modeling using EHR data: challenges, strategies, and a comparison of machine learning approaches,” _Med. Care_ , vol. 48, no. 6 Suppl, pp. S106–13, Jun. 2010.
  * [2] B. Shickel, P. J. Tighe, A. Bihorac, and P. Rashidi, “Deep EHR: A survey of recent advances in deep learning techniques for electronic health record (EHR) analysis,” _IEEE J Biomed Health Inform_ , vol. 22, no. 5, pp. 1589–1604, Sep. 2018.
  * [3] M. Chen, Y. Hao, K. Hwang, L. Wang, and L. Wang, “Disease prediction by machine learning over big data from healthcare communities,” _IEEE Access_ , vol. 5, pp. 8869–8879, 2017.
  * [4] P. S. Kohli and S. Arora, “Application of machine learning in disease prediction,” in _2018 4th International Conference on Computing Communication and Automation (ICCCA)_. ieeexplore.ieee.org, Dec. 2018, pp. 1–4.
  * [5] S. Uddin, A. Khan, M. E. Hossain, and M. A. Moni, “Comparing different supervised machine learning algorithms for disease prediction,” _BMC Med. Inform. Decis. Mak._ , vol. 19, no. 1, p. 281, Dec. 2019.
  * [6] D. H. Gustafson, “Length of stay: prediction and explanation,” _Health Serv. Res._ , vol. 3, no. 1, pp. 12–34, 1968.
  * [7] A. Awad, M. Bader-El-Den, and J. McNicholas, “Patient length of stay and mortality prediction: A survey,” _Health Serv. Manage. Res._ , vol. 30, no. 2, pp. 105–120, May 2017.
  * [8] V. A. Huynh-Thu, Y. Saeys, L. Wehenkel, and P. Geurts, “Statistical interpretation of machine learning-based feature importance scores for biomarker discovery,” _Bioinformatics_ , vol. 28, no. 13, pp. 1766–1774, Jul. 2012.
  * [9] D. Ledesma, S. Symes, and S. Richards, “Advancements within modern machine learning methodology: Impacts and prospects in biomarker discovery,” _Curr. Med. Chem._ , vol. 28, no. 32, pp. 6512–6531, 2021.
  * [10] A. Nicolae, G. Morton, H. Chung, A. Loblaw, S. Jain, D. Mitchell, L. Lu, J. Helou, M. Al-Hanaqta, E. Heath, and A. Ravi, “Evaluation of a Machine-Learning algorithm for treatment planning in prostate Low-Dose-Rate brachytherapy,” _Int. J. Radiat. Oncol. Biol. Phys._ , vol. 97, no. 4, pp. 822–829, Mar. 2017.
  * [11] G. Valdes, C. B. Simone, 2nd, J. Chen, A. Lin, S. S. Yom, A. J. Pattison, C. M. Carpenter, and T. D. Solberg, “Clinical decision support of radiotherapy treatment planning: A data-driven machine learning strategy for patient-specific dosimetric decision making,” _Radiother. Oncol._ , vol. 125, no. 3, pp. 392–397, Dec. 2017.
  * [12] K. C. Vranas, J. K. Jopling, T. E. Sweeney, M. C. Ramsey, A. S. Milstein, C. G. Slatore, G. J. Escobar, and V. X. Liu, “Identifying distinct subgroups of ICU patients: A machine learning approach,” _Crit. Care Med._ , vol. 45, no. 10, pp. 1607–1615, Oct. 2017.
  * [13] I. Pinal-Fernandez and A. L. Mammen, “On using machine learning algorithms to define clinically meaningful patient subgroups,” _Ann. Rheum. Dis._ , vol. 79, no. 10, p. e128, Oct. 2020.
  * [14] Y. Wang, Y. Zhao, T. M. Therneau, E. J. Atkinson, A. P. Tafti, N. Zhang, S. Amin, A. H. Limper, S. Khosla, and H. Liu, “Unsupervised machine learning for the discovery of latent disease clusters and patient subgroups using electronic health records,” _J. Biomed. Inform._ , vol. 102, p. 103364, Feb. 2020.
  * [15] A. D. Kaplan, J. D. Greene, V. X. Liu, and P. Ray, “Unsupervised probabilistic models for sequential electronic health records,” _J. Biomed. Inform._ , vol. 134, p. 104163, Oct. 2022.
  * [16] G. J. McLachlan and D. Peel, _Finite Mixture Models_. John Wiley & Sons, Mar. 2004\.
  * [17] A. Najjar, C. Gagné, and D. Reinharz, “Two-Step heterogeneous finite mixture model clustering for mining healthcare databases,” in _2015 IEEE International Conference on Data Mining_. ieeexplore.ieee.org, Nov. 2015, pp. 931–936.
  * [18] L. C. Cheung, Q. Pan, N. Hyun, M. Schiffman, B. Fetterman, P. E. Castle, T. Lorey, and H. A. Katki, “Mixture models for undiagnosed prevalent disease and interval-censored incident disease: applications to a cohort assembled from electronic health records,” _Stat. Med._ , vol. 36, no. 22, pp. 3583–3595, Sep. 2017.
  * [19] F. Stella and Y. Amer, “Continuous time bayesian network classifiers,” _J. Biomed. Inform._ , vol. 45, no. 6, pp. 1108–1119, Dec. 2012.
  * [20] Y.-Y. Liu, S. Li, F. Li, L. Song, and J. M. Rehg, “Efficient learning of Continuous-Time hidden markov models for disease progression,” _Adv. Neural Inf. Process. Syst._ , vol. 28, pp. 3599–3607, 2015.
  * [21] J. Futoma, S. Hariharan, and K. Heller, “Learning to detect sepsis with a multitask gaussian process RNN classifier,” in _International Conference on Machine Learning_. proceedings.mlr.press, 2017, pp. 1174–1182.
  * [22] A. M. Alaa and M. van der Schaar, “Bayesian inference of individualized treatment effects using multi-task gaussian processes,” _Adv. Neural Inf. Process. Syst._ , vol. 30, 2017.
  * [23] E. Choi, M. T. Bahadori, A. Schuetz, W. F. Stewart, and J. Sun, “Doctor AI: Predicting clinical events via recurrent neural networks,” _JMLR Workshop Conf. Proc._ , vol. 56, pp. 301–318, Aug. 2016.
  * [24] B. Jin, C. Che, Z. Liu, S. Zhang, X. Yin, and X. Wei, “Predicting the risk of heart failure with EHR sequential data modeling,” _IEEE Access_ , vol. 6, pp. 9256–9261, 2018.
  * [25] A. M. Alaa, J. Yoon, S. Hu, and M. van der Schaar, “Personalized risk scoring for critical care prognosis using mixtures of gaussian processes,” _IEEE Trans. Biomed. Eng._ , vol. 65, no. 1, pp. 207–218, Jan. 2018.
  * [26] S. Cui, E. C. Yoo, D. Li, K. Laudanski, and B. E. Engelhardt, “Hierarchical gaussian processes and mixtures of experts to model COVID-19 patient trajectories,” _Pac. Symp. Biocomput._ , vol. 27, pp. 266–277, 2022\.
  * [27] I. C. McDowell, D. Manandhar, C. M. Vockley, A. K. Schmid, T. E. Reddy, and B. E. Engelhardt, “Clustering gene expression time series data using an infinite gaussian process mixture model,” _PLoS Comput. Biol._ , vol. 14, no. 1, p. e1005896, Jan. 2018.
  * [28] A. D. Kaplan, Q. Cheng, K. A. Mohan, L. D. Nelson, S. Jain, H. Levin, A. Torres-Espin, A. Chou, J. R. Huie, A. R. Ferguson, M. Mccrea, J. Giacino, S. Sundaram, A. J. Markowitz, and G. T. Manley, “Mixture model framework for traumatic brain injury prognosis using heterogeneous clinical and outcome data,” _IEEE J Biomed Health Inform_ , vol. PP, Jul. 2021\.
  * [29] V. X. Liu, V. Fielding-Singh, J. D. Greene, J. M. Baker, T. J. Iwashyna, J. Bhattacharya, and G. J. Escobar, “The timing of early antibiotics and hospital mortality in sepsis,” _Am. J. Respir. Crit. Care Med._ , vol. 196, no. 7, pp. 856–863, Oct. 2017.
  * [30] T. K. Moon, “The expectation-maximization algorithm,” _IEEE Signal Process. Mag._ , vol. 13, no. 6, pp. 47–60, Nov. 1996.


# Light curves of the explosion of ONe WD+CO WD merger remnant and type Icn
supernovae

Chengyuan Wu Yunnan Observatories, Chinese Academy of Sciences, Kunming
650216, China  International Centre of Supernovae, Yunnan Key Laboratory,
Kunming, 650216, China<EMAIL_ADDRESS>Shuai Zha  Yunnan
Observatories, Chinese Academy of Sciences, Kunming 650216, China
International Centre of Supernovae, Yunnan Key Laboratory, Kunming, 650216,
China  Yongzhi Cai  Yunnan Observatories, Chinese Academy of Sciences, Kunming
650216, China  International Centre of Supernovae, Yunnan Key Laboratory,
Kunming, 650216, China  Zhengyang Zhang Yunnan Observatories, Chinese Academy
of Sciences, Kunming 650216, China  International Centre of Supernovae, Yunnan
Key Laboratory, Kunming, 650216, China  University of Chinese Academy of
Sciences, Beijing 100049, China  Yi Yang  Physics Department and Tsinghua
Center for Astrophysics, Tsinghua University, Beijing, 100084, People’s
Republic of China  Danfeng Xiang  Physics Department and Tsinghua Center for
Astrophysics, Tsinghua University, Beijing, 100084, People’s Republic of China
Weili Lin  Physics Department and Tsinghua Center for Astrophysics, Tsinghua
University, Beijing, 100084, People’s Republic of China  Xiaofeng Wang
Physics Department and Tsinghua Center for Astrophysics, Tsinghua University,
Beijing, 100084, People’s Republic of China  Bo Wang  Yunnan Observatories,
Chinese Academy of Sciences, Kunming 650216, China  International Centre of
Supernovae, Yunnan Key Laboratory, Kunming, 650216, China

###### Abstract

Type Icn supernovae (SNe Icn) are a newly detected rare subtype of interacting
stripped-envelope supernovae which show narrow P-Cygni lines of highly ionized
carbon, oxygen, and neon in their early spectra due to the interactions of the
SNe ejecta with dense hydrogen- and helium-deficient circumstellar material
(CSM). It has been suggested that SNe Icn may have multiple progenitor
channels, such as the explosion of carbon-rich Wolf–Rayet stars, or the
explosion of stripped-envelope SNe which undergo binary interactions. Among
the SNe Icn, SN 2019jc shows unique properties, and previous work inferred
that it may stem from the ultra-stripped supernova, but other possibilities
still exist. In this work, we aim to simulate the light curves from the
explosions of oxygen-neon and carbon-oxygen double white dwarf (WD) merger
remnants, and to further investigate whether the corresponding explosions can
appear as some particular SNe Icn. We generate the light curves from the
explosive remnants and analyse the influence of different parameters on the
light curves, such as the ejecta mass, explosion energy, mass of ${}^{\rm
56}{\rm Ni}$ and CSM properties. Comparing our results with some SNe Icn, we
found that the light curves from the explosions of double WD merger remnants
can explain the observable properties of SN 2019jc, which inferred that this
special SN Icn may have a different progenitor. Our results indicated that
double WD merger may be an alternative model in producing at least one of the
SNe Icn.

White dwarf stars (1799); Stellar mergers (2157); Supernovae (1668); Core-
collapse supernovae (304)

##  1 Introduction

It is widely believed that massive stars with initial masses greater than
$8-12{M}_{\odot}$ will end their lives exploding as core collapse supernovae
(CCSNe), which can be classified into different types from their spectral
properties (e.g., Filippenko 1997; Gal-Yam 2017). CCSNe that show hydrogen (H)
lines in their spectra are classified as Type II SNe, while for those that
show complete absence of H features in their spectra, or show evidence of H
only at early times are called H-poor or stripped-envelope SNe (SESNe),
including types Ib, Ic, and IIb SNe. Some CCSNe show evidence for interactions
between the SNe ejecta and circumstellar material (CSM) are named as
interacting SNe. The interacting SNe that exhibit bright and narrow Balmer
lines of H in their spectra are classified as Type IIn SNe, while the events
display narrow lines of helium (He) in their spectra without H, indicating the
presence of a H-poor but He-rich CSM around the exploding stars, are named as
Type Ibn SNe. Recently, an even more extreme subtype of interacting SESNe,
named as Type Icn supernovae (SN Icn) has been discovered (e.g., Fraser et al.
2021; Gal-Yam et al. 2022). The early-time spectra of SNe Icn are blue with
narrow P-Cygni lines of highly ionized carbon (C), oxygen (O) and neon (Ne),
instead of He. These SNe show rapid photometric evolution with a high peak
brightness. After the maximum light, as the SN ejecta becomes optically thin,
their spectra appear more similar to normal SNe Ibc or SNe Ibn but still
remain different (e.g. Gal-Yam et al. 2022; Perley et al. 2022; Davis et al.
2023).

Up to now, there are only five SNe Icn have been reported, i.e., SN 2019hgp,
SN 2019jc, SN 2021csp, SN 2021ckj and SN 2022ann (e.g., Fraser et al. 2021;
Perley et al. 2022; Gal-Yam et al. 2022; Pellegrino et al. 2022; Nagao et al.
2023; Davis et al. 2023). The origins of SNe Icn are still under debate. Gal-
Yam et al. (2022) proposed a Wolf-Rayet (WR) star to be the progenitor of SN
2019hgp based on its observational properties, and suggested a possibility
that the differences between Type Ibn and Type Icn SNe arise from their
different types of WR star progenitors, i.e., He- and nitrogen (N)-rich WN
stars for Type Ibn SNe and C-rich WC stars for Type Icn SNe. Perley et al.
(2022) suggested that a WR star directly collapse to a black hole can launch a
subrelativistic jet, and the jet interacts with the dense CSM may explain the
observational properties of SN 2021csp. Nagao et al. (2023) suggested that SN
2021csp and SN 2021ckj may have aspherical ejecta and spherical CSM, and the
different observational properties between these two SNe arise from the
viewing angle effects (see also Fraser et al. 2021), i.e., the SN 2021ckj is
observed from a direction close to the jet pole, while SN 2021csp is observed
from an off-axis direction. Davis et al. (2023) found that the CSM of SN
2022ann only has a velocity of about $800\,{\rm{km}}/{\rm s}$, which is
inconsistent with the typical velocity of WR star winds. Moreover, the host
galaxy of SN 2022ann is a low-mass dwarf galaxy with low star formation rate
and low averaged metallicity, suggesting that the progenitor of SN 2022ann has
low metallicity. This implies that it is difficult to fully stripped the H/He
envelope from its progenitor through stellar wind (massive WR star may not be
the progenitor of SN 2022ann). Instead, they suggested a binary-stripped
progenitor for SN 2022ann, rather than a single massive WR progenitor.
Pellegrino et al. (2022) analyzed a sample of four SNe Icn, including SN
2019hgp, SN 2019jc, SN 2021ckj, and SN 2021csp. They found that SN 2019jc
shows different features comparing with other SNe Icn in the following
aspects: first, this SN is unique in having a light curve that is both faster
evolving and fainter than those of other SNe Icn; secondly, SN 2019jc is
located at the outer edge of its host galaxy, where the star formation rate
density is low. Based on a circumstellar interaction model, Pellegrino et al.
(2022) fitted the bolometric light curves of these four SNe Icn, and found
that SN 2019jc has low explosion energy and low ejecta mass, suggesting a low
mass, ultra-stripped SN as its progenitor. However, for other SNe Icn, the WR
star progenitors or binary systems may better explain their observational
properties.

Recently, Wu et al. (2023) investigated the evolution of the post-merger
remnants from the coalescence of ONe and CO WD pairs. They found that the
merger remnants could evolve to H/He-deficient giant-like stars due to the
shell carbon burning. The merger remnants with masses greater than
$1.95{M}_{\odot}$ can explode as electron-capture SNe (ECSNe), and could
produce CO-rich CSM resulted from the stellar wind during their giant phases.
The interactions between the SNe ejecta and CSM may produce narrow C, O and Ne
emission lines similar to SNe Icn, making the double WD merger is able to be a
potential channel in producing SNe Icn. Even though some of the SNe Icn
observed at early times show He emission lines (for example, the spectrum of
SN 2019jc at maximum light shows He II $\lambda 4686$ feature), suggesting
that the CSM surrounding SNe Icn may not fully devoid of He. Note that there
could exist a He surface layer on the CO WD, which causes the merger remnant
involving at least one CO WD could contain amount of He (e.g., Dan et al.
2014; Schwab 2021; Wu et al. 2023). The subsequent evolution of such merger
remnant may produce CSM that is dominated by C, O and Ne but not completely
devoid of He through stellar wind, which means that it is possible to detect
He emission lines from the early-time spectra of some SNe Icn, if these SNe
are stemmed from the double WD merger channel. Since the optical light curves
from the explosions of the massive double WD merger remnants have not been
well studied, in this work, we are going to simulate the light curves from the
explosions of merger remnants, and then to compare them with some SNe Icn.

The rest of the article is organized as follows: in Sect. 2, we introduce the
method in generating the light curves of exploding double WD merger remnants.
In Sect. 3, we analyse how different sets of parameters in influencing the
light curves. In Sect. 4, we compare our simulated light curves to some SNe
Icn, and discuss whether our model can explain SNe Icn. Finally, the summary
are provided in Sect. 5.

##  2 Model setups

As investigated by previous works, ONe+CO WD merger remnants with masses
greater than $1.95{M}_{\odot}$ will end their lives as ECSNe (e.g., Wu et al.
2023). The collapse of the ONe core releases a large amount of gravitational
potential energy, which can be transformed to thermal and kinetic energy to
power a bright transient. In order to calculate the light curves from the
explosions of double WD merger remnants, we take the stellar structures at
their final evolutionary stages (just prior to the core collapse) from Wu et
al. (2023) as our initial models. Then we evolve the models from a few seconds
after the central explosions triggered by core collapse to a time just before
the outgoing shocks reach the stellar surface by employing the Modules for
Experiments in Stellar Astrophysics ($\tt{MESA}$, version 12778; e.g., Paxton
et al. 2011; Paxton et al. 2013; Paxton et al. 2015; Paxton et al. 2018;
Paxton et al. 2019). Subsequently, we use radiation hydrodynamics code
$\tt{STELLA}$ (e.g., Blinnikov & Bartunov 1993; Blinnikov & Panov 1996;
Blinnikov et al. 1998; Blinnikov et al. 2000; Blinnikov et al. 2006; Blinnikov
& Tolstov 2011) to evolve the models from the shock breakout to several tens
of days to TELLAgenerate the light curves.

In detail, our approach is as follows. First, we consider the most massive
merger remnant described in Wu et al. (2023) as our initial model (i.e.,
$2.25{M}_{\odot}$, resulted from an $1.20{M}_{\odot}$ ONe WD merge with an
$1.05{M}_{\odot}$ CO WD). The central core temperature of this model achieved
${\rm{log}}({T}_{\rm c}/{\rm K})=9.3$ when core collapse occurs, and the
radius of the remnant is about $377{R}_{\odot}$ at that moment. To investigate
the effect of ejecta mass on the light curves, we remove different central
sections from the initial model. We assume three mass cuts, i.e.,
$1.42{M}_{\odot}$ (just above the CO burning layer), $1.65{M}_{\odot}$ and
$1.89{M}_{\odot}$, corresponding to the ejecta masses equal to
$0.83{M}_{\odot}$, $0.60{M}_{\odot}$ and $0.36{M}_{\odot}$, respectively. The
mass profile and density structure of the initial model as well as the
positions of different mass cuts are shown in panel (a) and (b) of Fig. 1.

Secondly, we put the structures above the mass cuts into the explosion
calculations. According to the numerical explosion simulations, the explosion
energy of ECSNe is in an order of magnitude of ${10}^{50}\,{\rm{erg}}$ (e.g.,
Kitaura et al. 2006; Moriya & Eldridge 2016; Zha et al. 2022). To explore the
influence of explosion energy on the light curves, we inject different
energies (i.e., $0.05\,{\rm B}$, $0.10\,{\rm B}$ and $0.20\,{\rm B}$; $1\,{\rm
B}\equiv{10}^{51}\,{\rm{erg}}$) into a thin layer of about $0.01{M}_{\odot}$
at the inner boundary of the remaining envelope within $5$ milliseconds.

Thirdly, we calculate the shock propagating process from the innermost shell
of the envelope towards the surface. Since it is important to have enough
unshocked material for the $\tt{STELLA}$ code to calculate the light curves,
we stop our $\tt{MESA}$ simulations once there is $\sim 0.05{M}_{\odot}$
unshocked material remaining between the shock and the surface. During the
shock propagating process, we ignore all the nuclear reactions to keep the
elemental abundance distributions in the envelope remain unchanged. Note that
${}^{\rm 56}{\rm Ni}$ can be generated by the intense nuclear burning behind
the shock near the core-envelope boundary just after core collapse, and
numerical simulations indicated that the mass of ${}^{\rm 56}{\rm Ni}$
produced in ECSNe and accretion-induced-collapse events would be in the order
of magnitude of ${10}^{-3}$ to ${10}^{-4}\,{M}_{\odot}$ (e.g., Dessart et al.
2006; Moriya & Eldridge 2016). To explore how the mass of ${}^{\rm 56}{\rm
Ni}$ effects on the light curves, we artificially add different masses of
${}^{\rm 56}{\rm Ni}$ (i.e., $0.0001{M}_{\odot}$, $0.001{M}_{\odot}$ and
$0.01{M}_{\odot}$) in the envelope. This goal is achieved by reducing the mass
of ${}^{\rm 16}{\rm O}$ and replacing with ${}^{\rm 56}{\rm Ni}$. Besides,
previous works suggested that the ${}^{\rm 56}{\rm Ni}$ mixing process may
influence the sharp of light curves (e.g., Moriya et al. 2020). In our
simulations, we considered two different degrees of mixing. In the “full
mixing” model, ${}^{\rm 56}{\rm Ni}$ is uniformly mixed in the entire ejecta,
while in the “half mixing” model, the outer edge of the layer where ${}^{\rm
56}{\rm Ni}$ is mixed is in the center of the envelope. The elemental
abundance distribution profile of the ejecta with $0.001{M}_{\odot}$ of
${}^{\rm 56}{\rm Ni}$ is shown in panel (c) of Fig. 1.

Finally, after the outer edge of the shock reached the desired position, we
save the $\tt{STELLA}$ data by adopting “stella$\\_$num$\\_$points=$400$”,
where $400$ is a recommended value in $\tt{MESA}$ default. Note that the
stellar winds of the merger remnants prior to their explosions may form CSM
surround the progenitors, and the ejecta interacting with CSM could transform
kinetic energy into thermal energy and thus influence the light curves. To
investigate this effect, we attach the CSM on the top of the progenitors, and
use “stella$\\_$nz$\\_$extra=$40$” to save the data of CSM properties. We
consider the wind-like CSM. Based on the evolutionary results from the
previous work (e.g. Wu et al. 2023), the merger remnants could have enormous
wind mass-loss rates just prior to the SN explosions (the time intervals
between the significant wind mass-loss processes and the SN explosions are
within tens of days to several years, which are influenced by the total masses
of the merger remnants and the wind mass-loss prescriptions), which means that
the outer radii of the CSM surrounding the SNe are not too large. For the
wind-like CSM, the wind density parameter ($w$) is very useful in directly
estimating the density of the CSM (${\rho}_{\rm{CSM}}$), i.e.,

${w}=\frac{\dot{M}}{{v}_{\rm{wind}}},$ (1)

and

${\rho}_{\rm{CSM}}({\rm r})=\frac{w}{4\pi{r}^{-2}},$ (2)

where, $\dot{M}$ is the mass-loss rate of the progenitors, ${v}_{\rm{wind}}$
is the velocity of the stellar wind, and ${r}$ is the radius of CSM (e.g.,
Blinnikov 2017). Since the wind velocities and mass-loss rates of such H/He-
deficient giants are not well constrained from the observations, in the
current work, we consider various of wind velocities (${v}={10}\,{\rm km}/{\rm
s}$, ${100}\,{\rm km}/{\rm s}$ and ${1000}\,{\rm km}/{\rm s}$, corresponding
to the typical wind velocities of RGB/AGB stars, yellow supergiants and WR
stars, respectively) (e.g., Smith 2017), wind density parameters
(${3.15}\times{10}^{17}-{3.15}\times{10}^{19}\,{\rm g}/{\rm{cm}}$, depending
on the evolution of the progenitor stars under different wind mass-loss
parameters) and total masses of CSM (${M}_{\rm{CSM}}={0.05}$ and
$0.25{M}_{\odot}$). The density and velocity profiles of the remnants and CSM
(an example that we assume the enormous wind mass-loss rate, i.e.,
$0.5{M}_{\odot}/{\rm{yr}}$, occurs $0.1$ year prior to the SN explosion) are
provided in panel (d) of Fig. 1.

Figure 1: Properties of the merger remnant before light curves calculations.
Panel (a) and (b) represent the density profile and the enclosed mass as the
function of radius of the merger remnant before explosion. Blue squares, green
circles and red triangles represent the different mass cuts before energy
injection. Panel (c) and (d) represent the elemental abundance profile and the
density and velocity profiles of the $0.6{M}_{\odot}$ ejecta at the moment
when shock wave reached the position at $0.05{M}_{\odot}$ beneath the surface.
Panel (c) also shows the elemental abundance distribution in considering
${}^{\rm 56}{\rm Ni}$ with half and full mixing. Panel (d) also shows the
properties of CSM with different wind velocities (wind mass-loss rate of
$0.5{M}_{\odot}/{\rm{yr}}$ lasts for $0.1$ year prior to the SN explosion).

##  3 Bolometric Light Curves

We use the multi-group radiation hydrodynamics code $\tt{STELLA}$ to follow
the evolution of the explosion models described above through shock breakout
to $70$ days after explosion. $\tt{STELLA}$ is able to simulate SN evolution
at early times, and can also solve the shock breakout and interaction between
ejecta and CSM. In the calculation, SN explosion is regarded as a thermal bomb
locating at the inner boundary of the ejecta. $\tt{STELLA}$ solves the
radiative transfer equations in the intensity momentum approximation in each
frequency bin. We use $40$ frequency bins, which is enough to generate
bolometric light curves. The opacity is computed based on over $153000$
spectral lines from Kurucz & Bell (1995) and Verner et al. (1996). Some
simulated outputs from our $\tt{STELLA}$ simulations are publicly available on
Zenodo at doi:10.5281/zenodo.11066927.

We present the light curves from different progenitors described above in Fig.
2, while the evolution of blackbody temperature, density, radius, velocity and
opacity of a specific model (i.e., ${M}_{\rm{ej}}=0.60{M}_{\odot}$,
${E}=0.1\,{\rm B}$, without ${}^{\rm 56}{\rm{Ni}}$, without CSM) are presented
in Fig. 3. The overall evolutionary characteristics of the luminosity are as
follows (four evolutionary phases). (1) Rapid increasing phase: at the
beginning of the evolution, the ejecta is heated and ionized by the outwardly
propagating explosion shock wave. (2) Rapid declining phase: after the shock
breakout, the evolution of the bolometric luminosity is powered by the
diffusion of thermal energy depositing in the ejecta, and during this stage,
the luminosity declines quickly. (3) Plateau phase: as the ejecta expands and
cools, the temperature of the material gradually drops below the recombination
temperature and become largely transparent. The composition of the ejecta in
our model is mixed mainly by ${}^{\rm 12}{\rm C}$, ${}^{\rm 16}{\rm O}$,
${}^{\rm 20}{\rm{Ne}}$ and ${}^{\rm 24}{\rm{Mg}}$. For this mixture, the
recombination temperature is about $6000\,{\rm K}$. For the higher temperature
(${T}>6000\,{\rm K}$), the opacity of the mixture is greater than about
${0.04}\,{\rm{cm}}^{2}\,{\rm{g}}^{-1}$. When the temperature drops below
$6000\,{\rm K}$, oxygen recombines to neutral and the opacity drops
dramatically, which is similar to the O-Ne-Mg mixture (e.g., Kleiser & Kasen
2014). As the outer layers of the ejecta begin to recombines, a sharp
ionization front develops in the ejecta. At this moment, the photosphere stays
at the recombination front and the luminosity decline rate becomes
significantly smaller than before, resulting in the appearance of plateau-like
feature. (4) Late phase: the plateau phase lasts for about $10$-$25$ days
depending on the initial structures of the progenitors, and after which the
recombination wave reaches the center of the ejecta, exhausting the energy
stored in the ejecta and the light curve drops off extremely rapid, making the
end of the plateau phase.

The duration of the plateau phase is dominated by the dissipation rate of the
energy stored in the ejecta. From Fig. 2, we can see that the ejecta mass and
the explosion energy can significantly effect the plateau luminosity and the
duration of plateau phase (panel a and b). Under the same explosion energy, a
more massive ejecta mass leads to a lower plateau luminosity. This is because
when the recombination front developed in the envelope, the ejecta usually
have similar photosphere temperature, but a more massive ejecta has lower
expansion velocity, which causes a smaller photosphere radius. Besides, a
lower expansion velocity makes the energy deposited in the envelope needs more
time to be exhausted, hence, the plateau phase is longer for a more massive
ejecta mass. The effect of explosion energy on the light curves is similar to
that of ejecta mass, i.e., a higher explosion energy leads to a more luminous
but shorter plateau phase, and vice versa.

After the SN explosion, the nuclear burning behind the shock wave is able to
produce iron-group elements such as ${}^{\rm 56}{\rm Ni}$ in the ejecta. The
radioactive decay of ${}^{\rm 56}{\rm Ni}$ could heat the envelope and thus
increase the bolometric luminosity. In panel (c) of Fig. 2, we compare the
light curves of the progenitor models with different masses of ${}^{\rm
56}{\rm Ni}$. Various masses of ${}^{\rm 56}{\rm Ni}$ cannot influence the
light curves during the plateau phase significantly. About $15$ days after
explosion, the radioactive decay of isotopes gradually dominates the evolution
of the luminosity, and a ten times of ${}^{\rm 56}{\rm Ni}$ mass will cause an
order of magnitude increasing in the tail bolometric luminosity.
Theoretically, if ${}^{\rm 56}{\rm Ni}$ is entirely mixed in the ejecta, the
heating caused by radioactive decay in the outer layers of the ejecta becomes
efficient, which leads to the increase in luminosity occurs earlier than that
if ${}^{\rm 56}{\rm Ni}$ is only located at the inner portion of the ejecta.
However, due to the low ${}^{\rm 56}{\rm Ni}$ mass generated from ECSNe, the
plateau luminosity is still dominated by the energy deposited in the envelope,
which means that the degree of ${}^{\rm 56}{\rm Ni}$ mixing may not effect the
light curves significantly in the current situation.

After the merger of double WDs, the remnant can evolve to giant phase, and
during which the H/He-deficient CSM is formed due to the stellar wind. The
interaction between ejecta and CSM can transform kinetic energy into thermal
energy and thus alter the features of the light curves. The interacting
process is mainly influenced by the properties of the CSM such as its radius
and its density distribution. In panel (d) of Fig. 2, we present the light
curves of the explosion remnants with different CSM properties. Since the wind
mass-loss rates and wind velocities of such merger remnants are not well
constrained, we artificially added $0.05{M}_{\odot}$ or $0.25{M}_{\odot}$ of
CSM surrounding the progenitor. The density profile of the CSM is dominated by
the wind density parameters, and in these simulations, the parameter $w$ is in
the range between $3.15\times{10}^{17}$ and $3.15\times{10}^{19}\,{\rm
g}/{\rm{cm}}$, which is related to the wind mass-loss histories of the
progenitors. We consider three different wind velocities, i.e., ${v}=10\,{\rm
km}/{\rm s}$, ${v}=100\,{\rm km}/{\rm s}$ and ${v}=1000\,{\rm km}/{\rm s}$,
corresponding to the typical wind velocity of RGB/AGB stars, yellow
supergiants and WR stars, respectively (e.g., Smith 2017). From the panel (d)
of Fig. 2, we can see that the early bolometric luminosity becomes higher when
the wind density parameter of CSM is higher because of the extra radiation
from CSM interaction. A longer CSM interaction resulted from the larger radius
of the CSM causes the delayed appearance of peak luminosity in the light
curves. Afterwards, all the light curves restore to the similar evolutionary
paths, indicating the end of the CSM interaction.

Figure 2: Light curves from the explosions of double WD merger remnants with
different initial parameters, demonstrating how the parameters alter the light
curves. Panel (a): the effect of ejecta mass (${M}_{\rm ej}$). Panel (b): the
effect of explosion energy (${E}$). Panel (c): the effect of ${}^{\rm 56}{\rm
Ni}$ mass (${M}_{{}^{\rm 56}{\rm Ni}}$) and mixing. Panel (d): the effect of
CSM properties. Number (1)-(4) in panel (a) represent the four evolutionary
phases mentioned in Sect. 3. Figure 3: Additional properties of the merger
remnant (${M}_{\rm ej}=0.60{M}_{\odot}$, ${E}=1.0\times{10}^{50}\,{\rm erg}$)
after explosion. Panel (a) and (b): blackbody temperature, density, radius and
velocity as the function of time. Panel (c): Rosseland mean opacity at the
photosphere as a function of temperature.

##  4 Comparing to observation

To date, five SNe Icn have been reported. According to the analyses from the
previous works, it seems that this new type of SNe may have multiple formation
channels (e.g., Gal-Yam et al. 2022; Perley et al. 2022; Pellegrino et al.
2022). In order to investigate whether any SNe Icn may have originated from
the post-merger remnants resulted from the coalescence of ONe and CO WD pairs,
we obtain some observational data of three SNe Icn, i.e., SN 2019hgp, SN
2019jc and SN 2021csp from Pellegrino et al. (2022), and compare our models to
the observational data. The properties of the progenitors and the CSM are
based on the evolutionary results of $2.25{M}_{\odot}$ double WD merger
remnant provided in Wu et al. (2023). During the final evolutionary stages of
the merger remnants, there exist strong wind mass-loss processes, which lead
to the decrease in the radii of the progenitors and the formation of dense
CSM. However, note that it is difficult to accurately estimate the progenitor
radius and the mass of the CSM because of the uncertainty in wind mass-loss
rate. Hence, in our models, we assume the enormous stellar wind decreases the
radius of the progenitor down to $100{R}_{\odot}$, making the ejecta mass
equals to $0.468{M}_{\odot}$, and the dense CSM mass equals to
$0.025{M}_{\odot}$. Besides, we assume that the wind velocity equals to
$1000\,{\rm km}/{\rm s}$, which is consistent with the observed CSM velocity
of SNe Icn. We consider different explosion energies (from
$E=1.0\times{10}^{50}$ to $4.0\times{10}^{50}\,{\rm{erg}}$) and
$0.01{M}_{\odot}$ of ${}^{\rm 56}{\rm Ni}$ in the ejecta to simulate the
corresponding light curves. For comparison, we also simulate the light curves
with the same progenitor/CSM properties (${R}=100{R}_{\odot}$,
${M}_{\rm{ej}}=0.468{M}_{\odot}$ and ${E}=1.0\times{10}^{50}\,{\rm{erg}}$) but
without ${}^{\rm 56}{\rm Ni}$ in the ejecta or even without the interaction
between ejecta and CSM.

Figure 4: The comparison between light curves from the explosions of double WD
merger remnants and the bolometric light curves of three SNe Icn. In the light
curve simulations, the radius of progenitor is ${R}=100{R}_{\odot}$, the
ejecta mass is ${M}_{\rm{ej}}\approx 0.468{M}_{\odot}$, and the mass of nickel
${M}_{{}^{\rm 56}{\rm{Ni}}}=0.01{M}_{\odot}$ fully mixed into the ejecta. For
the property of CSM, the wind velocity is ${v}=1000{\rm{km}}/{\rm s}$, wind
density parameter is ${w}\approx 3.15\times{10}^{17}\,{\rm g}/{\rm{cm}}$, and
the mass of CSM is ${M}_{\rm{CSM}}=0.025{M}_{\odot}$. Solid lines with
different colors represent the light curves from our simulations with
different explosive energies (the CSM interaction and radioactive decay from
${}^{\rm 56}{\rm{Ni}}$ are considered). Dotted lines represent the simulated
light curves with explosive energy ${E}=0.1\,{\rm B}$, but without considering
radioactive decay from ${}^{\rm 56}{\rm{Ni}}$ (orange dotted line) and without
considering both CSM interaction and radioactive decay from ${}^{\rm
56}{\rm{Ni}}$ (black dotted line). For some of the models, the simulations
terminate at $23$ days after explosions in order to save the calculation time.
Green circles, blue squares and red triangles with error bars represent the
bolometric light curves of SN 2019jc, SN 2019hgp and SN 2021csp, respectively.

The comparison of our simulated light curves with these three SNe is presented
in Fig. 4, while the comparison of the evolution of blackbody temperature and
radius are presented in panel (a) and (b) of Fig. 5, and the multi-band light
curves from two of our models (${E}=2.0\times{10}^{50}$ and
$4.0\times{10}^{50}\,{\rm{erg}}$) are provided in panel (c) and (d) of Fig. 5.
The bolometric luminosities and blackbody properties of each objects are
estimated through the extensive multi-band observations (for more detail, see
Pellegrino et al. 2022). From Fig. 4, we can see that the bolometric
luminosities from the explosion of double WD merger remnants are significantly
lower that that of SN 2019hgp and SN 2021csp. This implies that these two SNe
should have more massive ejecta masses and higher explosion energies,
indicating that they may related to the explosions of massive WR stars.
Interestingly, the light curves from our models with explosion energy
${E}=2.0\times{10}^{50}-4.0\times{10}^{50}\,{\rm{erg}}$ provide reasonable
fits to the bolometric luminosity of SN 2019jc111Note that the progenitor/CSM
properties in our models are different with the estimated values from
Pellegrino et al. (2022). However, it is difficult to constrain some
properties such as ejecta mass, CSM mass, progenitor radius from the
observation directly, but the CSM velocity can be inferred from the spectra.
Thus, in our models, we assume the CSM velocity is $1000\,{\rm km}/{\rm s}$,
which is consistent with SNe Icn, but other parameters are related to the
evolution of the progenitors., which means that this SN may have a different
progenitor from other SNe Icn. Furthermore, SN 2019jc exploded on the
outskirts of the nearby spiral galaxy, where the star formation rate density
is low. The low star formation rate density argues against the progenitor with
high zero-age main-sequence mass. The progenitors of WDs are low mass main-
sequence stars and the SNe from the double WD mergers could form in any
positions of their host galaxies, which implies that the host galaxy
environment of SN 2019jc is also consistent with double WD merger channel.
However, it is still too early to make a conclusion that SN 2019jc is related
to double WD merger. On one hand, although the luminosities of B, V, R, I
bands are consistent with SN 2019jc, the evolution of multi-band light curves
obtained from our models are slower than that of SN 2019jc, which is probably
because of the discrepancy in the photosphere opacity resulted from the
difference of ejecta compositions or the radius of the progenitor. On the
other hand, SN 2019jc was not observed until approximately 3 days after
explosion and disappeared approximately 12 days after explosion, which cannot
constrain its early-time and subsequent evolution. For the early-time
evolution, our results imply that ultraviolet (UV) or X-ray radiation may have
major contribution to the bolometric light curves. For the late-time
evolution, the existence of ${}^{\rm 56}{\rm Ni}$ in the ejecta can
significantly influence the luminosity after the plateau phase. Hence, the
early-time UV or X-ray observations and the late-time observational data are
needed in constraining the progenitor channel of SN 2019jc. Currently, even
though we cannot provide any better limitations on the origin of SN 2019jc
depending on the limited data, we suggest that the explosion of double WD
merger remnant may be a new model in explaining the observational properties
of SN 2019jc, but other progenitor channels such as ultra-stripped SN or
stripped envelope of ECSN (e.g., Tauris et al. 2013; Tauris et al. 2015;
Moriya & Eldridge 2016) are still possible.

Figure 5: Panel (a) and (b): similar to Fig. 4, but for the evolution of
blackbody temperature (panel a) and blackbody radius (panel b). Panel (c) and
(d): multi-band light curves for double WD merger remnants with different
explosion energies.

##  5 Summary

In this work, we simulated the light curves from the explosion of ONe+CO WD
merger remnants. We found that the light curves decline quickly after the
shock breakout and enter into the plateau phases resulted from the
recombination of ionized carbon and oxygen. After the energy stored in the
envelope are exhausted, the light curves drop off rapidly. In the simulations,
we analysed the influences of different parameters such as ejecta mass,
explosion energy, ${}^{\rm 56}{\rm Ni}$ mass and different properties of CSM.
The ejecta mass and explosion energy can effect the plateau bolometric
luminosities and their durations. A higher explosion energy or a lower ejecta
mass leads to a more luminous but shorter plateau phase, and vice versa. The
${}^{\rm 56}{\rm Ni}$ mass can effect the tail bolometric luminosity. A ten
times of ${}^{\rm 56}{\rm Ni}$ mass will cause an order of magnitude
increasing in the tail bolometric luminosity. The density and the size of CSM
can effect the peak luminosity due to the CSM interaction. Furthermore, we
compare the simulated light curves with some SNe Icn. We found that the light
curves from the explosion of double WD merger remnant can fit the bolometric
luminosity of SN 2019jc within $12$ days after explosion, but significantly
fainter than that of SN 2019hgp and SN 2021csp. We proposed that double WD
merger may be an alternative progenitor channel of at least one SN Icn. In
order to better understand the progenitor systems of SN Icn and their
explosion mechanisms, more SN Icn samples with high-quality data (such as
early-time UV or X-ray observations) are needed in the future.

We thank the anonymous referee for his/her very helpful suggestions on the
manuscript. We thank Dongdong Liu, Zhenwei Li for helpful discussion. This
study is supported by the National Natural Science Foundation of China (Nos
12225304, 12288102, 12090040/12090043, 12303054), the National Key R&D Program
of China (No. 2021YFA1600404), the Yunnan Revitalization Talent Support
Program – Young Talent project, the Western Light Project of CAS (No. XBZG-
ZDSYS-202117), the science research grant from the China Manned Space Project
(No. CMS-CSST-2021-A12), the Frontier Scientific Research Program of Deep
Space Exploration Laboratory (No. 2022-QYKYJH-ZYTS-016), and the Yunnan
Fundamental Research Project (Nos 202301AU070039, 202201BC070003,
202301AU070109, 202401AU070063).

## References

  * Blinnikov (2017) Blinnikov, S. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 843, doi: 10.1007/978-3-319-21846-5_31
  * Blinnikov et al. (2000) Blinnikov, S., Lundqvist, P., Bartunov, O., Nomoto, K., & Iwamoto, K. 2000, ApJ, 532, 1132, doi: 10.1086/308588
  * Blinnikov & Bartunov (1993) Blinnikov, S. I., & Bartunov, O. S. 1993, A&A, 273, 106
  * Blinnikov et al. (1998) Blinnikov, S. I., Eastman, R., Bartunov, O. S., Popolitov, V. A., & Woosley, S. E. 1998, ApJ, 496, 454, doi: 10.1086/305375
  * Blinnikov & Panov (1996) Blinnikov, S. I., & Panov, I. V. 1996, Astronomy Letters, 22, 39
  * Blinnikov et al. (2006) Blinnikov, S. I., Röpke, F. K., Sorokina, E. I., et al. 2006, A&A, 453, 229, doi: 10.1051/0004-6361:20054594
  * Blinnikov & Tolstov (2011) Blinnikov, S. I., & Tolstov, A. G. 2011, Astronomy Letters, 37, 194, doi: 10.1134/S1063773711010051
  * Dan et al. (2014) Dan, M., Rosswog, S., Brüggen, M., & Podsiadlowski, P. 2014, MNRAS, 438, 14, doi: 10.1093/mnras/stt1766
  * Davis et al. (2023) Davis, K. W., Taggart, K., Tinyanont, S., et al. 2023, MNRAS, 523, 2530, doi: 10.1093/mnras/stad1433
  * Dessart et al. (2006) Dessart, L., Burrows, A., Ott, C. D., et al. 2006, ApJ, 644, 1063, doi: 10.1086/503626
  * Filippenko (1997) Filippenko, A. V. 1997, ARA&A, 35, 309, doi: 10.1146/annurev.astro.35.1.309
  * Fraser et al. (2021) Fraser, M., Stritzinger, M. D., Brennan, S. J., et al. 2021, arXiv e-prints, arXiv:2108.07278, doi: 10.48550/arXiv.2108.07278
  * Gal-Yam (2017) Gal-Yam, A. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 195, doi: 10.1007/978-3-319-21846-5_35
  * Gal-Yam et al. (2022) Gal-Yam, A., Bruch, R., Schulze, S., et al. 2022, Nature, 601, 201, doi: 10.1038/s41586-021-04155-1
  * Kitaura et al. (2006) Kitaura, F. S., Janka, H. T., & Hillebrandt, W. 2006, A&A, 450, 345, doi: 10.1051/0004-6361:20054703
  * Kleiser & Kasen (2014) Kleiser, I. K. W., & Kasen, D. 2014, MNRAS, 438, 318, doi: 10.1093/mnras/stt2191
  * Kurucz & Bell (1995) Kurucz, R. L., & Bell, B. 1995, Atomic line list
  * Moriya & Eldridge (2016) Moriya, T. J., & Eldridge, J. J. 2016, MNRAS, 461, 2155, doi: 10.1093/mnras/stw1471
  * Moriya et al. (2020) Moriya, T. J., Suzuki, A., Takiwaki, T., Pan, Y.-C., & Blinnikov, S. I. 2020, MNRAS, 497, 1619, doi: 10.1093/mnras/staa2060
  * Nagao et al. (2023) Nagao, T., Kuncarayakti, H., Maeda, K., et al. 2023, A&A, 673, A27, doi: 10.1051/0004-6361/202346084
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pellegrino et al. (2022) Pellegrino, C., Howell, D. A., Terreran, G., et al. 2022, ApJ, 938, 73, doi: 10.3847/1538-4357/ac8ff6
  * Perley et al. (2022) Perley, D. A., Sollerman, J., Schulze, S., et al. 2022, ApJ, 927, 180, doi: 10.3847/1538-4357/ac478e
  * Schwab (2021) Schwab, J. 2021, ApJ, 906, 53, doi: 10.3847/1538-4357/abc87e
  * Smith (2017) Smith, N. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 403, doi: 10.1007/978-3-319-21846-5_38
  * Tauris et al. (2013) Tauris, T. M., Langer, N., Moriya, T. J., et al. 2013, ApJ, 778, L23, doi: 10.1088/2041-8205/778/2/L23
  * Tauris et al. (2015) Tauris, T. M., Langer, N., & Podsiadlowski, P. 2015, MNRAS, 451, 2123, doi: 10.1093/mnras/stv990
  * Verner et al. (1996) Verner, D. A., Verner, E. M., & Ferland, G. J. 1996, Atomic Data and Nuclear Data Tables, 64, 1, doi: 10.1006/adnd.1996.0018
  * Wu et al. (2023) Wu, C., Xiong, H., Lin, J., et al. 2023, ApJ, 944, L54, doi: 10.3847/2041-8213/acb6f3
  * Zha et al. (2022) Zha, S., O’Connor, E. P., Couch, S. M., Leung, S.-C., & Nomoto, K. 2022, MNRAS, 513, 1317, doi: 10.1093/mnras/stac1035