# A Bernoulli Mixture Model to Understand and Predict Children’s Longitudinal
Wheezing Patterns

Pierre G. B. Moutounet-Cartan<EMAIL_ADDRESS>
Department of Mathematics
Imperial College London
London, SW7 2AZ, United Kingdom

###### Abstract

In this research, we perform discrete unsupervised machine learning through a
Bernoulli Mixture Model on data representing the expression of the wheeze
phenotype of patients at different stages of their childhood up to age 16.
Wheeze is a distinct noise produced while breathing due to narrowed airways,
such as asthma or viral chest infections. Due to a study from Henderson et
al., (2008), it has been estimated that around 23.5% of U.K. children had at
least wheezed once by six years of age, and 6.9% had persistent wheezing
problems. The usage of a Bernoulli Mixture Model is new in the field, where
previous classification methods used classical unsupervised learning such as
$K$-means, $K$-medoids, or Latent Class Analyses (Loza et al.,, 2016;
Kurukulaaratchy et al.,, 2014; Deliu et al.,, 2016; Brew et al.,, 2019). In
particular, Oksel et al., (2019) found that the Latent Class Analysis has
resulted majorly dependent to the sample size, and $K$-means is largely
dependent on the distance measure and so the data-set.

In this research, we estimate that around $27.99(\pm 2.15)\%$111The values in
brackets are the differences between the estimated values for the measures
listed and the left or right boundary of the 95% confidence interval. The
inverse cumulative distribution function value is taken from a
$t$-distribution with $n-1$ degrees of freedom. of the population has
experienced wheezing before turning 1 in the United Kingdom. Furthermore, the
Bernoulli Mixture Model classification is found to work best with $K=4$
clusters in order to better balance the separability of the clusters with
their explanatory nature, based on a cohort of $N=1184$. The probability of
the group of parents in the $j$th cluster to say that their children have
wheezed during the $i$th age is assumed $P_{ij}\sim\text{Beta}(1/2,1/2)$, the
probabilities of assignment to each cluster is
$R\sim\text{Dirichlet}_{K}(\alpha)$, the assignment of the $n$th patient to
each cluster is $Z_{n}\ |\ R\sim\text{Categorical}(R)$, and the $n$th patient
wheezed during the $i$th age is $X_{in}\ |\
P_{ij},Z_{n}\sim\text{Bernoulli}(P_{i,Z_{n}})$; where $i\in\\{1,\dots,6\\}$,
$j\in\\{1,\dots,K\\}$, and $n\in\\{1,\dots,N\\}$. The classification is then
performed through the E-M optimization algorithm (Bishop,, 2006; Saeed et
al.,, 2013). We found that this clustering method groups efficiently the
patients with late-childhood wheezing, persistent wheezing, early-childhood
wheezing, and none or sporadic wheezing. Furthermore, we found that this
method is not dependent on the data-set, and can include data-sets with
missing entries.

It is hoped this study will give medical staff an understanding of the
wheezing patterns in children up to age 16, and so provide an effective
treatment.

Keywords: Bayesian Statistics, E-M Optimization, Mixture Models, Bernoulli
Mixture Model, Wheeze, Respiratory Problems, Classification

###### Contents

  1. 1 The Cohort & Introduction
  2. 2 Theory of Bernoulli Mixture Model
  3. 3 Application to the data-set
  4. 4 Results
  5. 5 Conflicts of interest

##  1 The Cohort & Introduction

Parents of 1184 children residing in the United Kingdom were asked to tell if
they had wheeze in the previous year at ages 1, 3, 5, 8, 11, and 16. A
positive answer was recorded as a 1 in the registry, while a negative answer
was recorded as a 0. Each patient is given a fixed identification number going
from 0 to 1183. The data of 537 children was incomplete (45.4% of the data),
i.e., the parent did not answer to the question at least once. There were 647
full entries, i.e., parents of 647 children provided the information at every
stage described above. The information was later put into a spreadsheet, where
each row represented the data collected for each child, and each column the
answer for each of the periods as explained above, under ”Age 1,” ”Age 3,”
”Age 5,” ”Age 8,” ”Age 11,” and ”Age 16.” Whenever no information was
provided, NaN was recorded.

From this data, we can therefore estimate the number of children within the
U.K. population who had wheeze at most a year before turning 1, 3, 5, 8, 11,
and 16. We estimate the means by the sample means, and produce a 95%
confidence interval with the sample mean and the corresponding critical value
from a $t$-distribution with 647 degrees of freedom. This algorithm can be
found in Appendix.

Age period | Est. mean | LHS CI at 95% | RHS CI at 95%
---|---|---|---
Age 1 | 27.99% | 25.84% | 30.14%
Age 3 | 23.74% | 21.70% | 25.78%
Age 5 | 23.04% | 21.02% | 25.06%
Age 8 | 18.05% | 16.20% | 19.89%
Age 11 | 18.89% | 17.01% | 20.76%
Age 16 | 17.02% | 15.22% | 18.83%
Table 1: Estimate of the U.K. population having wheezed at different ages with
the LHS and RHS of the 95% confidence interval using a $t$-distribution with
647 degrees of freedom.

From Table 1, we can deduce that approximately $27.99(\pm 2.15)\%$ (with a 95%
confidence) of the population has experienced wheezing before turning 1. This
seems to corroborate Henderson et al., (2008) findings, where they estimated
that 23.5% of U.K. children had at least wheezed once by 6.

The challenge with this data is that 45.5% of it is incomplete. From Fig. 1,
the number of parents stopping to provide a status on their child wheeziness
increases whenever the latter gets older.

Figure 1: Number of NaN entries per age tranche within the cohort data.

In particular, other papers such as Loza et al., (2016); Kurukulaaratchy et
al., (2014); Deliu et al., (2016); Brew et al., (2019); Oksel et al., (2019)
have used well known clustering methods such as K-means, K-medoids, as well as
Latent Class Analysis. These methods require that no entry be missing or
classified as NaN in the data-set, which we have for a large proportion in our
case. Furthermore, Oksel et al., (2019) found that the Latent Class Analysis
has resulted majorly dependent to the sample size, the frequency, and the
timing of data collection.

Therefore, in this paper, we suggest a novel method for classifying the
longitudinal wheeze phenotype diagnosed in children aged 1 to 16, using a
Bernoulli Mixture Model based on the E-M optimization algorithm. This method
can allow missing entries in the data-set, and is generally less dependent on
the sample size than other methods.

##  2 Theory of Bernoulli Mixture Model

In this subsection, we consider the general case of the Bernoulli Mixture
Model so that we can apply it in our case later on.

Let $X_{1},\dots,X_{M}$ be a set of $M$ binary variables where
$X_{i}\sim\mathrm{Bernoulli}(\lambda_{i})$ for $\lambda_{i}\in(0,1)$, $\forall
i\in\\{1,2,\dots,M\\}$. Note that in this case, the marginals of $X_{i}$ are
given by $f_{X_{i}}(x_{i})=\lambda_{i}^{x_{i}}(1-\lambda_{i})^{1-x_{i}}$ where
$x_{i}\in\\{0,1\\}$, $\forall i\in\\{1,2,\dots,M\\}$.

Hence, we have that

$\mathbb{P}(\mathbf{X}|\Lambda)=\prod_{i=1}^{M}f_{X_{i}}(x_{i})=\prod_{i=1}^{M}\lambda_{i}^{x_{i}}(1-\lambda_{i})^{1-x_{i}}$

where $\mathbf{X}=(X_{1},\dots,X_{M})^{t}$ and
$\Lambda=(\lambda_{1},\dots,\lambda_{M})^{t}$.

As the above expression are written as a separate product of the different
probability mass functions of $X_{1},\dots,X_{M}$, we have that $\mathbf{X}$
is independent of $\Lambda$. We have
$\mathbb{E}(X_{i}|\lambda_{i})=\lambda_{i},\forall i\in\\{1,\dots,M\\}$ by
characteristic of the Bernoulli distribution, so that

$\mathbb{E}(\mathbf{X}|\Lambda)=\Lambda\in(0,1)^{M}\ldotp\ (\dagger)$

Similarly, by independence we have
$\mathrm{cov}(X_{i},X_{j}|\lambda_{i},\lambda_{j})=0$, $\forall
i,j\in\\{1,\dots,M\\}$ such that $i\neq j$, and by property of the Bernoulli
distribution, $\mathrm{var}(X_{i}|\lambda_{i})=\lambda_{i}(1-\lambda_{i})$,
$\forall i\in\\{1,\dots,M\\}$. Therefore,

$\mathbf{\Sigma}_{ik}:=\mathrm{cov}(\mathbf{X}|\Lambda)_{ik}=\delta_{ik}\lambda_{i}(1-\lambda_{i})\in(0,1)^{M\times
M}\ (\ddagger)$

where $\delta_{ik}$ is the Kronecker delta, i.e.,

$\delta_{ik}=\cases{1}\text{if}i=k,\\\ 0\text{if}i\neq k\ldotp$

We now build a finite mixture, with new parameters
$\Pi=(\pi_{1},\dots,\pi_{K})^{t}$, where $K\leqslant M$ will be the number of
clusters, by

$\mathbb{P}(\mathbf{X}|\tilde{\Lambda},\Pi):=\sum_{j=1}^{K}\pi_{j}\mathbb{P}(\mathbf{X}|\Lambda_{j})\
(\star)$

where $\tilde{\Lambda}=(\Lambda_{1},\dots,\Lambda_{K})^{t}$, the set of
parameters of each component (Saeed et al.,, 2013), where we have

$\mathbb{P}(\mathbf{X}|\Lambda_{j})=\prod_{i=1}^{M}\lambda_{ji}^{x_{i}}(1-\lambda_{ji})^{1-x_{i}}\ldotp$

For this mixture model, we have

$\mathbb{E}(\mathbf{X}|\tilde{\Lambda},\Pi)=\sum_{j=1}^{K}\pi_{i}\mathbb{E}(\mathbf{X}|\Lambda_{j})=\sum_{j=1}^{K}\pi_{j}\Lambda_{j}\in(0,1)^{M}\text{by}(\dagger)\ldotp$

Similarly,

$\mathrm{cov}(\mathbf{X}|\tilde{\Lambda},\Pi)=-\mathbb{E}(\mathbf{X}|\tilde{\Lambda},\Pi)\mathbb{E}(\mathbf{X}|\tilde{\Lambda},\Pi)^{t}+\sum_{j=1}^{K}\pi_{k}(\mathbf{\Sigma}_{j}+\Lambda_{j}\Lambda_{j}^{t})$

where $\mathbf{\Sigma}_{j}$ is as defined in $(\ddagger)$, i.e.,

$(\mathbf{\Sigma}_{j})_{ik}=\delta_{ik}\lambda_{ji}(1-\lambda_{ji})\ldotp$

For a random sample $\mathbf{X}_{1},\dots,\mathbf{X}_{N}$ which is distributed
as in $(\star)$ with respect to $\tilde{\Lambda}$ and $\Pi$, the likelihood
function for $\mathbf{X}_{1},\dots,\mathbf{X}_{N}|\tilde{\Lambda},\Pi$ is
given by

$\mathcal{L}(\mathbf{X}_{1},\dots,\mathbf{X}_{N},\tilde{\Lambda},\Pi)=\prod_{n=1}^{N}\sum_{j=1}^{K}\pi_{j}\mathbb{P}(\mathbf{X}_{n}|\Lambda_{j})$

so that the log-likelihood $\ell$ is

$\ell(\mathbf{X}_{1},\dots,\mathbf{X}_{N},\tilde{\Lambda},\Pi)=\sum_{n=1}^{N}\log\left(\sum_{j=1}^{K}\pi_{j}\mathbb{P}(\mathbf{X}_{n}|\Lambda_{j})\right)=\sum_{n=1}^{N}\log\left(\sum_{j=1}^{K}\left\\{\pi_{j}\prod_{i=1}^{M}\lambda_{ji}^{x_{ni}}(1-\lambda_{ji})^{1-x_{ni}}\right\\}\right)\ldotp$

The aim is to generate an algorithm that finds $\Pi$, $\tilde{\Lambda}$ which
maximize $\ell$, which we will call $\Pi_{\text{max}}$ and
$\tilde{\Lambda}_{\text{max}}$ respectively. Unfortunately, due to the shape
of $\ell$, $\Pi_{\text{max}}$ and $\tilde{\Lambda}_{\text{max}}$ cannot be
found in closed form in that case.

Hence, we introduce a latent binary variable
$\mathbf{Z}=(z_{1},\dots,z_{K})^{t}$ associated with each instance of
$\mathbf{X}$ (Bishop,, 2006). Therefore, the conditional distribution of
$\mathbf{X}$, given the latent variable $\mathbf{Z}$, is given by

$\mathbb{P}(\mathbf{X}|\mathbf{Z},\tilde{\Lambda})=\prod_{j=1}^{K}\mathbb{P}(\mathbf{X}|\Lambda_{j})^{z_{j}},$

where the prior for the latent variable is

$\mathbb{P}(\mathbf{Z}|\Pi)=\prod_{j=1}^{K}\pi_{j}^{z_{j}}\ldotp$

By considering the random sample $\mathbf{X}_{1},\dots,\mathbf{X}_{N}$, we
have the likelihood function $\mathcal{G}$ for
$\mathbf{X}_{1},\dots,\mathbf{X}_{N}|\mathbf{Z}$ given by

$\mathcal{G}(\mathbf{X}_{1},\dots,\mathbf{X}_{N},\mathbf{Z},\tilde{\Lambda},\Pi)=\prod_{n=1}^{N}\mathbb{P}(\mathbf{X}_{n},\mathbf{Z}|\tilde{\Lambda},\Pi)=\prod_{n=1}^{N}\prod_{j=1}^{K}(\pi_{j}\mathbb{P}(\mathbf{X}_{n}|\Lambda_{j}))^{z_{nj}},$

so that, once expanded, we have

$\mathcal{G}(\mathbf{X}_{1},\dots,\mathbf{X}_{N},\mathbf{Z},\tilde{\Lambda},\Pi)=\prod_{n=1}^{N}\prod_{j=1}^{K}\left(\pi_{j}\prod_{i=1}^{M}\lambda_{ji}^{x_{ni}}(1-\lambda_{ji})^{1-x_{ni}}\right)^{z_{nj}}\ldotp$

Therefore, the log-likelihood $g$ is given by

$g(\mathbf{X}_{1},\dots,\mathbf{X}_{N},\mathbf{Z},\tilde{\Lambda},\Pi)=\sum_{n=1}^{N}\sum_{j=1}^{K}z_{nj}\left(\log(\pi_{j})+\sum_{i=1}^{M}\left\\{x_{ni}\log(\lambda_{ji})+(1-x_{ni})\log(\lambda_{ji})\right\\}\right)\ldotp$

Hence, the expectation of the log-likelihood with respect to the marginal
distribution of $\mathbf{Z}$ is given by

$\mathbb{E}_{\mathbb{P}(\mathbf{Z}|\Pi)}(g(\mathbf{X}_{1},\dots,\mathbf{X}_{N},\mathbf{Z},\tilde{\Lambda},\Pi))=\sum_{n=1}^{N}\sum_{j=1}^{K}\mathbb{E}(z_{nj})\left(\log(\pi_{j})+\sum_{i=1}^{M}\left\\{x_{ni}\log(\lambda_{ji})+(1-x_{ni})\log(\lambda_{ji})\right\\}\right)\ldotp(\star\star)$

By Bayes’ Theorem and Saeed et al., (2013); Bishop, (2006), we have

$\mathbb{E}(z_{nj})=\dfrac{\pi_{j}\mathbb{P}(\mathbf{X}_{n}|\Lambda_{j})}{\sum_{k=1}^{K}\pi_{k}\mathbb{P}(\mathbf{X}_{n}|\Lambda_{k})}\ldotp\
(\diamond)$

Let

$N_{j}=\sum_{n=1}^{N}\mathbb{E}(z_{nj})\text{and}\mathbf{\overline{X}}_{j}=\dfrac{1}{N_{j}}\sum_{n=1}^{N}\mathbb{E}(z_{nj})\mathbf{X}_{n}\ldotp$

Then one can show that

$\tilde{\Lambda}_{\text{max}}=(\hat{\Lambda}_{1},\dots,\hat{\Lambda}_{K})^{t}\text{where}\hat{\Lambda}_{j}=\mathbf{\overline{X}}_{j}$
$\text{and}\hat{\Pi}_{\text{max}}=(\hat{\pi}_{1},\dots,\hat{\pi}_{K})^{t}\text{where}\hat{\pi}_{j}=N_{j}/N\ldotp$

Indeed, such a $\tilde{\Lambda}_{\text{max}}$ makes the derivative with
respect to the $\lambda_{j}$s vanish, and such a $\hat{\Pi}_{\text{max}}$
maximizes $(\star\star)$ through a Lagrange multiplier as seen in Saeed et
al., (2013); Bishop, (2006).

Therefore, the expectation-maximization algorithm for a Bernoulli mixture
model first gives initialization values to $\tilde{\Lambda}$ and $\Pi$. The
algorithm then computes the value of the log-likelihood at the initial values
$\tilde{\Lambda}_{0}$, $\Pi_{0}$.

On the next step, the algorithm does a loop by evaluating $\mathbb{E}(z_{nj})$
as in $(\diamond)$, and reevaluates $\hat{\Lambda}_{j},\hat{\pi}_{j}$ as found
above. Then, we evaluate the log-likelihood at these values. We stop the loop
whenever the log-likelihood meets a convergence criterion.

##  3 Application to the data-set

Each cluster has a certain probability of a yes answer at the different ages
of the children. We assume that they follow a Beta distribution. Hence, let

$P_{ij}\sim\mathrm{Beta}\left(\dfrac{1}{2},\dfrac{1}{2}\right)\text{for}i\in\\{1,\dots,6\\}\text{and}j\in\\{1,\dots,K\\}$

where $K$ is the number of clusters, as the Beta distribution is well known to
represent a distribution of probabilities, and $\mathrm{Beta}(1/2,1/2)$ is
more dense around its support boundaries $0$ and $1$. Hence, this is our
prior. That is, given a cluster $j\in\\{1,\dots,K\\}$ and the age tranche of
the children $i\in\\{1,\dots,6\\}$, the group of parents in the $j$th cluster
has a probability of $P_{ij}$ to say that their children have wheezed at (or
within a year of) the $i$th age.

Now, to provide the cluster assignments, we use the Categorical distribution,
which is a generalized Bernoulli distribution (in higher dimensions), as each
vector in the data either belongs to one of the $K$ clusters or does not.
Hence, let

$Z_{n}\ |\ R\sim\mathrm{Categorical}_{K}(R)\text{for}n\in\\{1,\dots,N\\},$

where $R$ represents the probabilities of the cluster assignments, which we
assume follows a Dirichlet distribution, i.e., a generalized Beta distribution
(in higher dimensions), so that

$R\sim\mathrm{Dirichlet}_{K}(\alpha)$

for some constant and positive vector $\alpha$, which we are free to choose as
it is an uninformative distribution.

Therefore, for $X_{in}$ the random variable representing the fact that the
$n$th patient wheezed at the $i$th age tranche, we have

$X_{in}\ |\
P_{ij},Z_{n}\sim\mathrm{Bernoulli}(P_{i,Z_{n}})\text{for}i\in\\{1,\dots,6\\}\text{and}n\in\\{1,\dots,N\\}\ldotp$

To sum up, the model is given by P_ij &∼Beta(1/2,1/2),
R ∼Dirichlet_K(α),
Z_n — R ∼Categorical_K(R),
X_in — P_ij, Z_n ∼Bernoulli(P_i,Z_n), for $i\in\\{1,\dots,6\\}$,
$j\in\\{1,\dots,K\\}$, and $n\in\\{1,\dots,N\\}$, where each variable
describes the following: P_ij&: probability of the group of parents in the jth
cluster to say that their children have
wheezed during the ith age,
R: probabilities of assignment to each cluster 1,…,K,
Z_n: assignment of the nth patient to each cluster 1,…,K,
X_in: the nth patient wheezed during the ith age.

The plate diagram for such a configuration can be found in Fig. 2. Note again
that $\alpha$, a $K$-dimensional vector with $\alpha_{j}>0$ for all
$j\in\\{1,\dots,K\\}$, can be chosen arbitrarily since the distribution of $R$
is uninformative ($Z_{n}$ is latent).

Figure 2: Plate diagram of the considered Bernoulli Mixture Model for the
considered cohort.

##  4 Results

We then perform the expectation-maximization algorithm on the log-likelihood
as explained in the previous section. We can then plot the Hinton plots for
$R$, the $P_{ij}$s, and the $Z_{n}$s. These are visualized in Fig. 3.

(a) The Hinton plot for $4$ clusters of $R$.

(b) The Hinton plot for $5$ clusters of $R$.

(c) The Hinton plot for $6$ clusters of $R$.

(d) The Hinton plot for $7$ clusters of $R$.

Figure 3: The Hinton plots of $R$ for different number of clusters. Each graph
shows the probabilities of being assigned to each cluster. Here, the data used
excludes the rows with NaN entries.

Hinton diagrams for $R$ show the number of elements per cluster. The bigger
the white square is, the bigger the cluster is. From Fig. 3, we can see that
for either setting $4,5,6$ or $7$ groups, there are $4$ dominant clusters.
Now, if we set the number of clusters to $64=2^{6}$, which is the number of
all possibilities to arrange 0s and 1s in a 6-dimensional vector, i.e., all
ways parents can answer to the questionnaire, then we get the Hinton diagram
for $R$ as shown in Fig. 4.

Figure 4: The Hinton diagram by setting $64$ clusters of $R$. Figure 5: Hinton
diagram of the $P_{ij}$s for $i\in\\{1,\dots,6\\}$ and $j\in\\{1,\dots,4\\}$,
i.e., the probabilities of the parents saying that their child wheezed at each
age tranche separated by cluster. The squares in the lower rectangle show the
”sizes” of each cluster (more formally, the probability of belonging to each
cluster).

We can see from Fig. 4 that there are $4$ dominant clusters. Hence, later on
for the Bernoulli mixture model we will take 4 different clusters. Note that
under the implemented algorithm, although the clusters – of course – do not
change by launching again the code, the order of these clusters do change. A
choice of 4 clusters give the following Hinton diagram for the $P_{ij}$s, as
seen in Fig. 5, superimposed with the Hinton diagram for $R$. According to
Fig. 5, Cluster 2 is made of the children with high probabilities of late
wheezing, Cluster 4 of the children with high probabilities of early wheezing,
Cluster 3 of the children with high probabilities of persisting wheezing, and
Cluster 1 of the children with probabilities of sporadic and benign wheezing.
To give ourselves an idea of the assignments of each cluster for each patient,
we can plot the cluster heat map as shown in Fig. 6(a). As seen in Fig. 6(a),
this clustering method seems to be grouping quite efficiently the patients
with late-childhood wheezing (blue group), persistent wheezing (green group),
early-childhood wheezing (red group), and none or sporadic wheezing (purple
group).

(a) Rows with missing entries disregarded.

(b) Rows with missing entries kept.

Figure 6: Heat map of the answers of the parents and their respective cluster
(the cluster they have been assigned with highest probability) where (6(a))
the rows with missing entries were disregarded, and where (6(b)) the rows with
missing entries were all kept. Each dark blue entry represents a positive
answer by the parents that their child had wheezed, while a light blue entry
is a negative answer. From Fig. 5, Cluster 1 is here the purple group, Cluster
2 is the blue group, Cluster 3 is the green group and Cluster 4 is the red
group. The gray-scaled bar represents the probability of assignment to this
cluster, going from $0.4$ (white) to $1$ (black). Each missing answer for Fig.
6(b) from the parents is shown in gray.

The answers in each group that look more like outliers do have a lower
probability to be assigned this cluster, as seen classified in Fig. 6(a). For
example, the set of answers (No, No, Yes, No, No, Yes) belongs to the purple
cluster with smallest probability of assignment in this cluster, at
approximately 46% (but still higher than the probabilities of assignment in
other clusters). The highest in the purple cluster being 97% for the set of
answers (No, No, No, No, No, No). For the blue cluster, the smallest
assignment probability is approximately 51% with set of answers (Yes, No, Yes,
No, Yes, Yes) while the highest is 97% with answers (No, No, No, Yes, Yes,
Yes) – which makes sense as it is the cluster of late-wheezing children. For
the green cluster, the smallest is 51% with answers (Yes, No, Yes, Yes, No,
No), the highest being at 99.9% with answers (Yes, Yes, Yes, Yes, Yes, Yes)
since it is the cluster of persistent wheeze issues children. Finally, for the
red cluster, the smallest probability is 61% for the answers (No, Yes, No, No,
No, Yes), the highest being 98% for the answers (Yes, Yes, No, No, No, No) as
it is the cluster of early-wheezing children.

Figure 7: Sankey diagram showing the patients flow between Bernoulli mixture
models with the complete data only (left) and all the data including rows with
missing entries (right). Note that the clusters here may not have the same
labels and color as before since the algorithm was called another time here.

We can also perform the Bernoulli process on the full data, i.e., also on the
parents that did not say if their child wheezed at all age tranches, since
$P_{ij}$ is a probability assigned to each answer individually. The aim is to
compare whether the change of grouping of the data with answers for each age
tranche changes whether or not we add up this new incomplete data. Ideally, it
does not so that the prediction is accurate. The heat map of the clustering
via Bernoulli mixture of all the data, including the missing one, is shown in
Fig. 6(b).

Comparing the sizes of the clusters in Fig. 6(a) and Fig. 6(b), adding the
data with missing entries has increased the sizes of the green and blue
clusters. We can here plot a Sankey diagram, as seen in Fig. 7 showing to
which clusters the entries with complete data are mapped from a Bernoulli
mixture with only complete data to a Bernoulli mixture with all the data
available. From Fig. 7, we can see that the model for full entries remains
stable with only seven swing patients, of which one of type $\alpha$ (Yes, No,
Yes, No, Yes, Yes), two of type $\beta$ (Yes, Yes, No, No, Yes, No), one of
type $\gamma$ (Yes, No, Yes, No, No, Yes), two of type $\xi$ (No, No, No, Yes,
No, Yes), and one of type $\lambda$ (No, Yes, No, No, No, Yes).

Those swing patients are shown in Table 2. The ”biggest” swing is made by the
two patients of type $\beta$ (see Table 2), where they have above 60% of
probability to be assigned different clusters when the missing entries from
the full set of data is added up. This is due to the fact that the expression
of the wheeze phenotype of these patients is very unpredictable and do not
show any pattern – the parents answered (Yes, Yes, No, No, Yes, No).

Patient type | Cl L | Pr L | Cl R | Pr R
---|---|---|---|---
$\alpha$ | 3 | 51% | 2 | 55%
$\beta$ | 1 | 67% | 2 | 65%
$\gamma$ | 4 | 64% | 2 | 40%
$\xi$ | 4 | 56% | 3 | 57%
$\lambda$ | 1 | 41% | 4 | 54%
Table 2: Allocations of the swing patients from a Bernoulli mixture with only
the clean data (assigned to the cluster Cl L with probability Pr L) to all the
data (assigned to the cluster Cl R with probability Pr R), where each cluster
number is as shown in Fig. 7.

The Bernoulli Mixture Model is advantageous in practice because it allows
practitioners to add up data throughout the childhood of a patient. E.g., if
the parents have a child aged 3 and gave answers to the questionnaire (whether
or not the child wheezed before 1, and around 3), the data can still be
clustered through the Bernoulli Mixture Model, so that the practitioner can
look at complete past entries in this group (”neighboring data”) and predict
future wheezing patterns of this child. As seen in Fig. 6(b), the Bernoulli
Mixture Model generally assumes higher probabilities of wheeze if no entry is
given, unless if it is sporadic. This method might be more accurate if the
parents are asked if their children wheezed during smaller age intervals, such
as yearly.

##  5 Conflicts of interest

The author(s) note no known conflict of interest by undertaking this research.

## References

  * Bishop, (2006) Bishop, C. M. (2006).  pages 444–455.  Information Science and Statistics. Springer (New York, NY).  ISBN: 978-03-873-1073-2.
  * Brew et al., (2019) Brew, B. K., Chiesa, F., Lundholm, C., Örtqvist, A., and Almqvist, C. (2019).  A modern approach to identifying and characterizing child asthma and wheeze phenotypes based on clinical data.  Plos One, 14(12).
  * Deliu et al., (2016) Deliu, M., Sperrin, M., Belgrave, D., and Custovic, A. (2016).  Identification of Asthma Subtypes Using Clustering Methodologies.  Pulmonary Therapy, 2(1):19–41.
  * Henderson et al., (2008) Henderson, A. J., Granell, R., Heron, J., et al. (2008).  Associations of wheezing phenotypes in the first 6 years of life with atopy, lung function and airway responsiveness in mid-childhood.  Thorax, 63(11):974–980.
  * Kurukulaaratchy et al., (2014) Kurukulaaratchy, R. J., Zhang, H., Raza, A., Patil, V., Karmaus, W., Ewart, S., and Arshad, S. H. (2014).  The Diversity of Young Adult Wheeze; a Cluster Analysis in a Longitudinal Birth Cohort.  Journal of the British Society for Allergy and Clinical Immunology, 44(5):724–735.
  * Loza et al., (2016) Loza, M. J., Djukanovic, R., Chung, K. F., Horowitz, D., Ma, K., Branigan, P., Barnathan, E. S., Susulic, V. S., and et al (2016).  Validated and longitudinally stable asthma phenotypes based on cluster analysis of the ADEPT study.  Respiratory Research, 17(165).
  * Oksel et al., (2019) Oksel, C., Granell, R., Mahmoud, O., Custovic, A., Henderson, A. J., Investigators, S., and Investigators, B. T. (2019).  Causes of variability in latent phenotypes of childhood wheeze.  The Journal of Allergy and Clinical Immunology, 143(5):1783–1790.
  * Saeed et al., (2013) Saeed, M., Javed, K., Babri, H. A., et al. (2013).  Machine learning using Bernoulli mixture models: Clustering, rule extraction and dimensionality reduction.  Neurocomputing, 119(7):366–374.


# Multispacecraft Remote Sensing and In Situ Observations of the 2020 November
29 Coronal Mass Ejection and Associated Shock: From Solar Source to
Heliospheric Impacts

Chong Chen  State Key Laboratory of Space Weather, National Space Science
Center, Chinese Academy of Sciences, Beijing 100190, People’s Republic of
China<EMAIL_ADDRESS>University of Chinese Academy of Sciences, Beijing
100049, People’s Republic of China  Ying D. Liu  State Key Laboratory of Space
Weather, National Space Science Center, Chinese Academy of Sciences, Beijing
100190, People’s Republic of China<EMAIL_ADDRESS>University of Chinese
Academy of Sciences, Beijing 100049, People’s Republic of China  Bei Zhu
Space Engineering University, Beijing 101416, People’s Republic of China

###### Abstract

We investigate the source eruption, propagation and expansion characteristics,
and heliospheric impacts of the 2020 November 29 coronal mass ejection (CME)
and associated shock, using remote sensing and in situ observations from
multiple spacecraft. A potential–field source–surface model is employed to
examine the coronal magnetic fields surrounding the source region. The CME and
associated shock are tracked from the early stage to the outer corona using
extreme ultraviolet and white light observations. Forward models are applied
to determine the structures and kinematics of the CME and the shock near the
Sun. The shock shows an ellipsoidal structure, expands in all directions, and
encloses the whole Sun as viewed from both SOHO and STEREO A, which results
from the large expansion of the CME flux rope and its fast acceleration. The
structure and potential impacts of the shock are mainly determined by its
radial and lateral expansions. The CME and shock arrive at Parker Solar Probe
and STEREO A. Only based on the remote sensing observations, it is difficult
to predict whether and when the CME/shock would arrive at the Earth. Combining
Wind in situ measurements and WSA-ENLIL simulation results, we confirm that
the far flank of the CME (or the CME leg) arrives at the Earth with no shock
signature. These results highlight the importance of multipoint remote sensing
and in situ observations for determining the heliospheric impacts of CMEs.

Interplanetary shocks (829) — Solar wind (1534) — Solar coronal mass ejections
(310)

##  1 Introduction

Coronal mass ejections (CMEs) are large–scale magnetized plasma ejected from
the Sun into interplanetary space, which are responsible for major geomagnetic
storms in the terrestrial environment. They are called interplanetary CMEs
(ICMEs), when traveling into interplanetary space. A fast CME can drive a
shock ahead of it, when the CME speed relative to the ambient solar wind is
greater than the magnetosonic speed or Alfvén speed of the ambient solar wind.
CME–driven shocks can accelerate solar energetic particles (SEPs) and enhance
the geo–effectiveness of CMEs. Understanding the process of CME/shock
propagation and evolution is of critical importance for space weather
forecasting.

Multispacecraft remote sensing and in situ observations are required to study
the complete chain of the evolution of CMEs/shocks from their solar sources
through heliospheric propagation to their impacts. As far as we know, the
study of such a complete chain is still lacking. First, the use of
multi–perspective remote sensing observations can better constrain the
three–dimensional (3D) structure of the CME and its associated shock. To
determine the 3D morphology of a CME, Thernisien et al. (2006, 2009) develop a
graduated cylindrical shell (GCS) model based on the coronagraph images from
different vantage points. As for CME–driven shocks, Kwon et al. (2014, 2015)
propose an ellipsoid model to reconstruct the 3D structure of a shock. Second,
separated multipoint in situ measurements can better assess the structure and
the heliospheric impacts of CMEs and shocks (e.g., Liu et al. 2008, 2019; Hu
et al. 2017; de Lucas et al. 2011). The multispacecraft measurements at
different locations can determine the propagation and expansion
characteristics of CMEs/shocks along different directions, including the
evolution anisotropy of CMEs/shocks resulting from different background solar
wind. Multipoint in situ measurements are also needed to verify the
forecasting accuracy along different directions of a magnetohydrodynamics
model in a same event simulation (e.g., Reinard et al. 2012; Biondo et al.
2021).

Finally, the combination of remote sensing and in situ observations from
multiple spacecraft can provide a more complete picture of CME/shock
evolution, and improve our capabilities to interpret events and to forecast
space weather effects (e.g., Liu et al. 2008, 2011, 2012; Nieves-Chinchilla et
al. 2012; Hu et al. 2016; Palmerio et al. 2021). Below we show examples of how
the multi–perspective remote sensing observations in combination with in situ
measurements are used to investigate the 3D propagation and evolution of
CMEs/shocks. Based on the wide–angle imaging observations from Solar
Terrestrial Relations Observatory (STEREO; Kaiser et al. 2008) and in situ
measurements at 1 AU, a geometric triangulation technique has been developed
to track CMEs (Liu et al., 2010a, b), and the Sun–to–Earth kinematics of fast
and slow CMEs have been revealed (Liu et al., 2013, 2016). Analyzing
multipoint imaging observations and widely separated in situ measurements, Liu
et al. (2017, 2019) find the fading and persistence of CME–driven shocks along
different directions in the heliosphere. In addition, using multipoint remote
sensing and in situ observations, the interactions of successive CMEs, which
may cause intense geomagnetic storms, are examined during the whole
propagation process (e.g., Liu et al. 2012, 2014b, 2014a; Webb et al. 2013;
Mishra et al. 2015; Lugaz et al. 2017). The geometrical and magnetic
relationships between ICMEs at 1 AU, CMEs near the Sun and their solar sources
have also been investigated, using coordinated imaging and in situ
observations from multiple vantage points (e.g., Liu et al. 2010b, 2011; Xie
et al. 2021; Marubashi et al. 2015; Syed Ibrahim et al. 2019; Savani et al.
2015). However, the connections between in situ signatures and solar source
characteristics are still not fully understood, because CMEs can have a very
complex evolution from the Sun to 1 AU, including deceleration, deflection,
rotation, erosion, and interaction with ambient solar wind structures and
other CMEs (e.g., Manchester et al. 2017 and references therein). Therefore,
the comparisons between multi–perspective remote sensing observations and
multipoint in situ measurements are necessary and critical to make the correct
connections between the solar source eruptions, CMEs near the Sun and their in
situ counterparts and to understand their heliospheric impacts.

On 2020 November 29, a large CME erupted, which is associated with a shock and
an M4.4 class flare that peaked at about 13:11 UT. The CME caused the first
type II radio burst and the first widespread SEP event of solar cycle 25
(Lario et al., 2021; Kollhoff et al., 2021). The eruption is a limb event as
viewed from the Earth, and is observed by a fleet of spacecraft at different
vantage points. The positions of the spacecraft in the ecliptic plane on 2020
November 29 are shown in Figure 1. STEREO A and Parker Solar Probe (PSP; Fox
et al. 2016) are 0.96 AU and 0.81 AU from the Sun, and 57.8° and 96.2° east of
the Earth, respectively. The red arrow indicates the location of the source
region (E97°S25°), which points to PSP in the ecliptic plane. This event
provides a unique opportunity to investigate the solar source, propagation,
and heliospheric impacts of the CME and its shock, as multi–perspective remote
sensing observations and multipoint in situ measurements are available for
this event. We examine source eruption signatures in Section 2, the
propagation and expansion characteristics of the CME/shock in Section 3, and
in situ measurements at different locations in Section 4. The results are
summarized in Section 5. The results of this work provide a more complete
propagation and expansion picture of the CME and its shock, and highlight the
importance of multispacecraft remote sensing and in situ observations for
determining the heliospheric impacts of CMEs.

##  2 source eruption signatures

The CME was launched around 12:45 UT on 2020 November 29 associated with
significant features, such as an M4.4 class flare, coronal dimmings, and EUV
waves as shown in Figure 2. Panel (a) shows an overview of the source region
in a 195 Å image from the Extreme–Ultraviolet Imager (EUVI) of the Sun Earth
Connection Coronal and Heliospheric Investigation (SECCHI; Howard et al. 2008)
aboard STEREO A with potential–field source–surface (PFSS) extrapolation
results mapped on it. The viewpoint of the image deviates a few degrees from
STEREO A (see the black part on the left of the EUVI image) to show the
coronal magnetic field configuration around the source region. The closed
magnetic field arches are toroidal and surrounded by open magnetic field lines
of the same polarity. The source region is beneath the right part of the
toroidal closed magnetic field arches as marked by the white arrow in panel
(a).

Panels (b)-(c) display STEREO A/EUVI 195 Å running–difference images at two
moments during the time of the eruption. The solar flare is visible at the
beginning of the eruption followed by a significant dimming region, which is
indicative of the removal of the coronal plasma. The dimming region
continuously extends. Panel (b) also shows the significant rise and expansion
of the plasma loop. The most dramatic signature in the EUV observations is the
EUV wave around the dimming, which propagates away from the source region and
sweeps a large area of the solar disk, as shown in panel (c). The continuously
extended areas of the dimming and the EUV wave indicate the expansion of the
CME as it rises up.

Panel (d) shows a hot channel at its nascent stage, which has been interpreted
as a flux rope (e.g., Zhang et al. 2012; Cheng et al. 2013, 2014) although the
magnetic field is not observed, in a 131 Å image observed by the Atmosphere
Imaging Assembly (AIA; Lemen et al. 2012) on board Solar Dynamics Observatory
(SDO; Pesnell et al. 2012). This is a limb event as viewed from the Earth, and
the source region is just behind the solar limb, so we can track the hot
channel at the early stage using EUV observations (see below). The axis of the
hot channel possibly has a largely inclined angle, because the two legs of the
hot channel separate with a large distance. SDO/AIA 193 Å running–difference
images are displayed in panels (e)-(f). In panel (e), the plasma loop is
bubble–shaped ahead of the hot channel shown in panel (d). There is an
ejection propagating toward the south as marked by the red arrow in panel (f).
This is not the eruption we study, because of its propagation direction. This
eruption may be caused by the EUV wave, which sweeps across and destabilizes
it.

These observations suggest that the hot channel has a large tilt angle, and,
as the hot channel rises and expands, it causes a dimming region and a
large–scale EUV wave. The EUV wave propagates across a significant portion of
the solar disk and is recognized as the footprint of the shock (see below).

##  3 structure and kinematics of the CME/shock

EUV and coronagraph observations of the eruption from GOES 16, STEREO A, and
Solar and Heliospheric Observatory (SOHO; Domingo et al. 1995) are displayed
in Figure 3. Panels (a)-(c) show the eruption at the early stage in EUV and
white light images. The primary EUV wave is considered to be the footprint of
the CME–driven shock (e.g., Patsourakos & Vourlidas 2012; Cheng et al. 2012;
Kwon et al. 2014), and at the early time the fronts of the CME and shock are
not well separated (see panel (c)). Therefore, they can be used together to
determine the shock structure. The well–developed CME and shock in coronagraph
images from two vantage points are shown in panels (d) and (g). The shock can
be seen in the images as a faint edge around the CME leading edge (e.g., Liu
et al. 2008; Hess & Zhang 2014; Zhu et al. 2021). We can see the backward
propagation of the shock on the opposite side of the eruption, because it
expands in all directions. Previous studies suggest that the shock can be
modeled well using a simple spheroidal structure (e.g., Hess & Zhang 2014;
Kwon et al. 2014; Liu et al. 2017, 2019). Here, we use a geometrical ellipsoid
model developed by Kwon et al. (2014) to determine the 3D morphology of the
shock, based on the EUV and white light observations from GOES 16, STEREO A,
and SOHO. The ellipsoid model has seven free parameters: the height,
longitude, and latitude of the center of the ellipsoid; the lengths of the
three semiprincipal axes; and the rotation angle of the ellipsoid. We assume
the cross section of the shock ellipsoid perpendicular to the propagation
direction to be a circle, which reduces two free parameters (one semiprincipal
axis and the rotation angle) in our fitting. We start to fit the shock at
12:56 UT when the EUV wave can be distinguished from the brightenings in the
dimming region, and get a good visual consistency between the model and
observations. At the initial stage, the footprint of the ellipsoid model on
the solar disk is consistent with the EUV wave front (see panels (a)-(b)). In
coronagraph images, as shown in panels (c), (f), and (i), the shock front is
well represented by the ellipsoid model. As for the CME, we employ the GCS
model proposed by Thernisien et al. (2006) to fit the CME based on
running–difference coronagraph images from STEREO A and SOHO. The GCS model
can determine the direction of propagation, tilt angle of CME flux rope and
height. By adjusting the free parameters, the GCS model fits the CME
observations from the two viewpoints very well, as shown in panels (e) and
(h).

Figure 4 shows the 3D modeled CME and shock structures at 15:18 UT on 2020
November 29. At this moment, the propagation longitude of the CME and shock is
about 85° and 80° east of the Sun–Earth line as marked by the blue and red
arrows, respectively. Application of the GCS model gives an average
propagation direction of about 90° east of the Sun–Earth line and about 17°
south of the ecliptic plane. The propagation direction of the CME is roughly
consistent with the location of the source region (E97°S25°), and does not
change much within 30 $R_{\sun}$. The tilt angle of the CME flux rope obtained
from the GCS model is about 75° with respect to the ecliptic plane, which is
consistent with the AIA 131 Å observations (see Figure 2 (d)). When we employ
the GCS model to fit the CME (e.g., Thernisien et al. 2006, 2009; Liu et al.
2010b), the half angle increases from 15° at 2.4 $R_{\sun}$ to 70° at 15.8
$R_{\sun}$, demonstrating the large expansion of the flux rope at its initial
stage. The speed of the CME leading edge is accelerated from $\sim$800 km s-1
to $\sim$2100 km s-1 within 10 $R_{\sun}$ (see below). Therefore, the
ellipsoidal structure of the shock is produced by the large expansion of the
CME flux rope and its fast acceleration. The front of the shock nose is close
to the CME nose, while the flank has farther distances from the CME. The shock
surrounds the whole Sun, which helps to explain the detection of SEPs at Solar
Orbiter behind the location of the eruption (see Figure 1). The CME and shock
could arrive at PSP and STEREO A, according to their propagation directions
and large expansions.

We can track the CME flux rope and the shock from their initial stage to the
outer corona with projection effects minimized, as the propagation direction
is almost perpendicular to the Sun–Earth line. The translation and expansion
distances of the shock can be obtained from the ellipsoid model fitting. As
shown in Figure 4, the lateral and radial expansion distances are represented
by “b” and “c”, respectively. The “d” and “d + c” denote the translation
distance of the shock center and shock nose height, respectively. The motion
of the CME flux rope can be tracked by stacking the EUV running–difference
images within a slit along the radial direction, as marked by the white line
in the panel (d) of Figure 2. The distance–time map of 131 Å
running–difference images from SDO/AIA is shown in the left panel of Figure 5.
The heights of the hot channel are extracted along the track, as indicated by
the red dashed curve. The right panel of Figure 5 shows the time–elongation
map produced by stacking running–difference images of C2 and C3 from
SOHO/LASCO within a slit along the ecliptic plane. The track associated with
the CME is marked by the red dashed curve. The track in C2 is not clear,
because the CME is so fast that C2 only captures three images. Elongation
angles of the CME leading edge in the ecliptic plane are extracted along the
track. Given that SOHO is about 90° away from the propagation direction of the
CME, we use a harmonic mean (HM) approximation (Lugaz et al., 2009), which
assumes that CMEs are attached to the Sun as a spherical front and move along
a fixed radial direction, to derive the CME kinematics. Readers are directed
to Liu et al. (2013) for discussions of selection of CME geometry depending on
the observation angle of the spacecraft.

The height and speed of the CME obtained from the GCS modeling, HM
approximation, and EUV observations, and the kinematics of the shock are
displayed in Figure 6. Note that we have used the GCS propagation longitude
(E90°) as input for the HM approximation. As shown in the top panel, the CME
height from the HM approximation can connect with the GCS model. The height of
the shock nose is slightly larger than the CME height, because the front of
the shock nose is ahead of and close to the CME leading edge (see Figure 4).
The radial distance of the CME flux rope obtained from the EUV observations
could be connected with the height of the CME from the GCS model. The speeds
of the CME and shock nose are shown in the middle panel. The speed profiles of
the CME and shock nose are similar. The CME and shock nose are accelerated to
about 2100 km s-1 within only 1 hour below $\sim$6 $R_{\sun}$. Then the speeds
have a deceleration. According to the propagation characteristics of fast CMEs
described by Liu et al. (2013), the speeds probably have a gradual longer
deceleration in interplanetary space to meet the shock speed at STEREO A
(about 700 km s-1), as marked by a horizontal dashed line in the middle panel.
The bottom panel shows the speeds of the lateral expansion, radial expansion
and translation of the shock. The expansion speeds of the shock are much
larger than the translation speed of the shock center. Therefore, the shock
radial expansion provides a major contribution to the propagation of the shock
nose, and the structure and potential impacts of the shock are mainly
determined by its radial and lateral expansions. We also overlap the GOES 1–8
Å X–ray flux on the middle panel scaled in arbitrary units. The liftoff of the
CME flux rope and the rise of the X–ray flux are almost at the same time. All
the speeds increase during the flare rising phase, and reach maxima after the
peak value of the X–ray flux.

Figure 7 shows the radio dynamic spectra associated with the 2020 November 29
CME from PSP, STEREO A, and Wind. The intense type III radio busts, which
start at about 12:55 UT on November 29 within several minutes after the
liftoff of the CME, are observed by all three spacecraft. An intense type III
radio bust portends the occurrence of a major CME on the Sun (e.g., Reiner et
al. 2001). All the spacecraft observe a short–duration type II radio burst in
the initial phase, which starts at about 13:10 UT almost near the peak of the
X–ray flux, as shown in the panels zoomed in. We use the Leblanc density model
(Leblanc et al., 1998) with an electron density of 15 cm-3 at 1 AU to convert
the shock nose distances obtained from the ellipsoid model to frequencies. The
Leblanc density model describes the average radial variation of the density of
the medium where the shock was propagating, so the electron density of 15 cm-3
is the nominal density and not necessary to be observed at 1 AU. The
frequencies are doubled to their harmonic frequencies and overlapped on the
spectra. The frequencies roughly match the observed type II radio burst when
the corresponding shock nose distances are below $\sim$ 5 $R_{\sun}$. The
vertical dashed lines in the PSP spectrum represent the shock arrival times at
PSP from in situ measurements. The previous one is associated with the small
CME on 2020 November 26. The second shock is the one we study in this paper,
and is also observed at STEREO A. The arrival times of the second shock at PSP
and STEREO A correspond to the suddenly enhanced, diffuse intensity in the
spectra. Note that there is no shock signature observed at Wind.

To determine the propagation and heliospheric impacts of the CME/shock, we
request a WSA–ENLIL simulation (Arge et al., 2004; Odstrcil et al., 2004) run
from the Community Coordinated Modeling Center (CCMC). CME parameters derived
from the GCS model, such as the time across the inner boundary at 21.5
$R_{\sun}$, half angular width, radial speed, and propagation direction, are
input in the simulation. Note that, as mentioned above, a small CME occurred
on 2020 November 26 and was also observed by PSP. Here, we only insert the CME
of interest in the simulation, and the results are good compared with the in
situ measurements. The CME/shock arrival times at PSP and Wind predicted by
the WSA–ENLIL model are consistent with the observed arrival times from in
situ measurements. As for STEREO A, the predicted arrival time is only about 5
hours earlier than the observed shock arrival time. The density and velocity
profiles from the simulation generally agree with the in situ measurements,
although the simulation overestimates the density. The simulation results are
displayed in Figure 8. The top panels show the density and radial speed
distributions at 01:01:06 UT on 2020 December 1. The CME arrives at PSP and
STEREO A with a large angle width. The density and radial speed of the CME are
much larger than the ambient solar wind. From the bottom panels, the far flank
of the CME (or the CME leg) could arrive at the Earth with relatively low
density and speed about one day after arriving at STEREO A.

##  4 in situ measurements

Figure 9 displays the solar wind magnetic field measurements from the FIELDS
instruments (Bale et al., 2016) aboard PSP, and the WSA–ENLIL simulation
results at PSP. From PSP in situ measurements, there are two shocks and two
ICME structures observed at PSP. The first shock passes PSP around 23:05 UT on
November 29, which may be associated with the small CME event on November 26
propagating toward PSP. The second shock is the one of interest, which is
observed at PSP around 18:35 UT on November 30, followed by a long–duration
sheath region and an ICME structure. The long–duration sheath region is from
18:35 UT on November 30 to 02:23 UT on December 1 lasting about 8 hours. This
indicates that PSP likely observes the flank of the CME, because the standoff
distance between the shock front and the CME flux rope increases from the nose
to the flank (see Figure 4). The propagation direction of the CME in latitude
and the location of the source region also imply that the northern flank of
the CME would arrive at PSP and STEREO A. The ICME interval is from 02:23 UT
to 11:40 UT on December 1, which is determined by the bidirection of the
electron pitch angle distribution (PAD) and magnetic field profiles, such as
the smooth and strong magnetic field and the coherent rotation of the field.
Although the resolution of the electron PAD is low, there are some indications
of the bidirectional streaming electron (BDE) signatures. We can see the
mainly northward $B_{N}$ component in the ICME implying a large tilt angle of
the CME flux rope. This is consistent with the results obtained from the GCS
model. The WSA–ENLIL simulation results are plotted as red dashed lines.
Because there are no plasma data from PSP during the period of time, we can
not compare the simulation results with the in situ plasma measurements. The
shock arrival time predicted by the WSA–ENLIL model is consistent with the
observed shock arrival time. Therefore, we use the WSA–ENLIL simulation
results to estimate the shock speed at PSP, which is about 850 km s-1. The
simulated magnetic field is much lower than the observed, because in the
WSA–ENLIL model, the inserted CME has no internal magnetic field and carries
the field from the ambient solar wind. From the in situ measurements at PSP,
the second shock has not propagated into the first ICME yet.

The in situ measurements from STEREO A associated with the 2020 November 29
CME are presented in Figure 10, with the simulation results at STEREO A shown
as red dashed lines. A shock is observed at STEREO A around 07:20 UT on
December 1, which is $\sim$5 hours later than the predicted arrival time from
the simulation. The shock speed predicted by the WSA–ENLIL model is about 850
km s-1, which is larger than the observed shock speed at STEREO A (about 700
km s-1). The shock is followed by a sheath region lasting about 4 hours, and
the ICME interval is determined by the electron PAD and magnetic field
profiles. The magnetic field strength of the ICME from STEREO A shows a
declining profile similar to the observations from PSP. However, the magnetic
field components are more dynamic than that observed by PSP.

If we only consider the location of the source region and the propagation
direction of the CME/shock, the CME and shock may not arrive at the Earth,
because the Earth is about 90° away from the propagation direction of the
CME/shock. Figure 11 shows the Wind in situ measurements probably associated
with the CME, which seem to be not like a typical ICME (e.g., Zurbuchen &
Richardson 2006). Given the large expansion of the CME and shock near the Sun,
and the simulation suggests that the CME may arrive at the Earth (see Figure
8), we infer that the structure observed at Wind may be a part of the CME. The
lack of typical ICME signatures (e.g., no clear BDE signature) probably
suggests that the far flank of the CME (or the CME leg) arrives at the Earth.
According to the proton temperature and magnetic field strength profiles, we
give an approximate interval of the ICME, as indicated by two vertical dashed
lines. The density profile from the simulation generally agrees with the in
situ measurements, although the WSA–ENLIL simulation overestimates the density
about 5 times. The CME arrival time (07:33 UT on December 2) predicted by the
WSA–ENLIL model is also in good agreement with the measurements. These
indicate that the CME indeed arrives at the Earth. According to the CME and
shock structures in the coronagraph images and the simulation results, we
suggest that the far flank of the CME (or the CME leg) arrives at the Earth
due to the vast expansion of the CME. However, there is no shock signature in
the measurements. This supports the perspective of Liu et al. (2017) that, at
some point, the shock would be just a wave and lose the nonlinear steepening
character especially near the wake.

##  5 summary and conclusions

We have investigated the source eruption, propagation and expansion
characteristics, and heliospheric impacts of the 2020 November 29 CME/shock,
combining multi–perspective remote sensing observations and multipoint in situ
measurements. EUV observations together with PFSS modeling are used to examine
the source eruption signatures and the coronal magnetic field configuration.
The structures and kinematics of the CME and associated shock are analyzed
using the forward models and remote sensing observations from multiple
spacecraft. The WSA–ENLIL simulation and the multipoint in situ measurements
are used to study the heliospheric impacts of the CME/shock. Key results are
obtained concerning the source eruption, interplanetary propagation, and
heliospheric impacts of the CME and associated shock.

The source region is beneath toroidal closed magnetic field arches, which are
surrounded by open magnetic field lines of the same polarity. The CME eruption
is associated with a significant dimming and an EUV wave sweeping a large area
of the solar disk. This indicates a large CME eruption associated with a
shock, and implies a rapid expansion of the CME/shock. This is a limb event as
viewed from the Earth, and the CME flux rope can be seen in AIA 131 Å images,
whose axis has a large tilt angle. Therefore, we can track the flux rope at
its early stage using EUV observations.

Due to the large expansion of the CME flux rope and its fast acceleration, the
shock shows an ellipsoidal structure, expands in all directions, and encloses
the whole Sun as viewed from both SOHO and STEREO A. We use the GCS model and
the ellipsoid model to fit the CME and the shock, respectively. The tilt angle
and propagation direction of the CME flux rope obtained from the GCS model are
consistent with the EUV observations. The CME flux rope quickly expands and
accelerates at its initial stage. Therefore, the ellipsoidal structure of the
shock is produced by the large expansion of the CME flux rope and its fast
acceleration. The ellipsoid model fitting results suggest that the shock
expands in all directions and encloses the whole Sun as viewed from both SOHO
and STEREO A, which helps to explain the detection of SEPs at Solar Orbiter
behind the location of the eruption. The structure and kinematics of a shock
can be used to study the relationship between SEPs arriving at different
locations and their sources (e.g., Rouillard et al. 2011; Lario et al. 2017;
Zhu et al. 2021). The study of the 2020 November 29 widespread SEP event has
been performed by Kollhoff et al. (2021), but they did not give a detailed
analysis of the CME and associated shock kinematics. Therefore, the shock
ellipsoid fitting results (see Table 1) in our paper can be a complement and
provide some key parameters for the study of the SEP acceleration.

The structure and potential impacts of the shock are mainly determined by its
radial and lateral expansions. We track the CME by producing jmaps with
projection effects minimized, as the propagation direction is almost
perpendicular to the Sun–Earth line. The kinematics of CME/shock from
different methods agree with each other. The speeds of the CME and shock nose
are rapidly accelerated to maxima within only 1 hour below $\sim$6 $R_{\sun}$,
and then show a deceleration, which is a typical speed profile for fast CMEs
(Liu et al., 2013). The shock expansion speeds are much larger than the
translation speed, so we suggest that the structure and potential impacts of
the shock are mainly determined by its radial and lateral expansions. The
liftoff of the CME flux rope and the rise of the X–ray flux are almost at the
same time, and all the speeds of the CME and shock increase during the flare
rising phase and peak after the maximum value of the X–ray flux.

The type III and type II radio bursts associated with the eruption are
observed at PSP, STEREO A, and Wind. All three spacecraft observe an intense
type III radio burst, which indicates a major CME on the Sun, and a
short–duration type II radio burst in the initial phase. We convert the shock
nose distances obtained from the ellipsoid model to frequencies using the
Leblanc density model with an electron density of 15 cm-3 at 1 AU. The shock
kinematics seem consistent with the frequency drift of the associated type II
burst.

According to the separated in situ measurements and simulation results, the
CME and shock arrive at PSP and STEREO A, and the Earth observe the far flank
of the CME (or the CME leg) due to the large expansion of the CME. The shock
is observed at PSP around 18:35 UT on November 30, followed by a long–duration
sheath region and an ICME structure. The ICME has a mainly northward $B_{N}$
component implying a large tilt angle of the CME flux rope, which is
consistent with the results obtained from the GCS model. The shock passes
STEREO A around 07:20 UT on December 1. The magnetic field strength of the
ICME from STEREO A shows a declining profile similar to the observations from
PSP, but the magnetic field components are more dynamic than that observed by
PSP. The Earth is about 90° away from the CME/shock nose, so according to the
location of the source region and the propagation direction of the CME and
shock, they may not arrive at the Earth. Moreover, it is difficult to predict
whether and when the CME/shock would arrive at the Earth only based on the
remote sensing observations. Combining the in situ measurements from Wind and
WSA–ENLIL simulation results, we finally find that, because of the large
expansion of the CME, the far flank of the CME (or the CME leg) arrives at the
Earth with no shock signature. This is surprising because the shock is
generally thought to be wider than the CME. The shock may have already decayed
before reaching Wind along this far flank direction. These results highlight
that multispacecraft remote sensing and in situ observations are important for
determining the heliospheric impacts of CMEs.

The research was supported by NSFC under grant 41774179, Beijing Municipal
Science and Technology Commission (Z191100004319003), the Specialized Research
Fund for State Key Laboratories of China, and the Strategic Priority Research
Program of Chinese Academy of Sciences (XDA15018500). We acknowledge the use
of data from Parker Solar Probe, STEREO, SDO, SOHO, GOES, and Wind. The
WSA–ENLIL simulations are provided by CCMC through their public Runs on
Request system (http//ccmc.gsfc.nasa.gov).

Figure 1: Positions of the spacecraft and planets in the ecliptic plane on
2020 November 29. The “A”, “B”, and “SO” represent the satellites of STEREO A,
STEREO B, and Solar Orbiter respectively. The black circle and red ellipse
mark the orbits of the Earth and PSP respectively. The dotted curves show
Parker spiral magnetic fields created with a solar wind speed of 450 km s-1.
The red arrow indicates the direction of the source region in the ecliptic
plane.

Figure 2: Magnetic field configuration and EUV observations of the source
region. (a) PFSS modeled coronal magnetic fields surrounding the source region
mapped onto the EUVI 195 Å image at 13:14 UT on 2020 November 29. The blue
lines represent the closed field lines, and the green lines are the open
magnetic field lines. The white arrow marks the source region. (b)-(c)
Running–difference images of STEREO A/EUVI at 195 Å showing the low coronal
signatures associated with the CME. (d) Hot channel at 131 Å from SDO/AIA. The
white line indicates a slice along the radial direction to create a
distance–time diagram for the hot channel. (e)-(f) Running–difference images
of SDO/AIA at 193 Å showing the plasma loop and the EUV wave from the
viewpoint of the Earth. The red arrow marks the eruption caused by the EUV
wave.

Figure 3: Running–difference images from GOES 16, STEREO A, and SOHO with
corresponding modeling of the CME and shock. Panels (a)-(c) show the same
shock modeling result for 13:15 UT. Panels (d)-(i) display the
running–difference coronagraph images and corresponding GCS modeling (green)
and ellipsoid shock modeling (red).

Figure 4: Geometry of the CME (blue lines) and the shock (red lines) at 15:18
UT on 2020 November 29. The blue and red arrows indicate the propagation
directions of the CME and shock, respectively. The directions of the Earth,
STEREO A, PSP, and Solar Orbiter are marked by the black arrows. The center of
the shock ellipsoid (O) and its distances from the center of the Sun (d), from
the nose of the shock (c) and from the flank perpendicular to the propagation
direction (b) are also indicated.

Figure 5: Left: distance–time map of the hot channel by stacking SDO/AIA 131 Å
running–difference images along the white line in the panel (d) of Figure 2.
The red curve indicates the track of the hot channel. Right: time–elongation
map constructed from running–difference images of C2 and C3 from SOHO/LASCO
along the ecliptic plane. The red curve marks the track of the CME leading
edge, along which the elongation angles are extracted.

Figure 6: Kinematics of the CME and shock derived from different methods. Top:
the distances of the CME derived from EUV observations, the GCS model, and the
HM approximation, and the height of the shock nose derived from the ellipsoid
model. Middle: speeds of the CME and shock nose derived from the numerical
differentiation of the distances in the top panel. GOES 1–8 Å X–ray flux is
overlapped in the panel scaled in arbitrary units. The horizontal dashed line
represents the shock speed of $\sim$700 km s-1 at STEREO A. Bottom: expansion
speeds of the shock along the radial and lateral directions and translation
speed of the shock center. The distances and speeds of the CME front are
binned to reduce the uncertainties, so the CME parameters are average values
and standard deviations within the bins. Following Kwon et al. (2014), the
uncertainty of the ellipsoid model parameters is estimated to be 8%.

Figure 7: Radio dynamic spectra associated with the 2020 November 29 eruption
from PSP, STEREO A, and Wind. The white arrows mark the type II radio burst
observed by the spacecraft. The shock nose distances derived from the
ellipsoid model are converted to frequencies by using the Leblanc density
model with an electron density of 15 cm-3 at 1 AU. The frequencies are
overlapped in the spectra as red crosses. The white curve in the middle panel
is the GOES X–ray flux scaled in arbitrary units. The vertical dashed lines
represent the shock arrival times from in situ measurements at PSP and STEREO
A. The areas with radio bursts are expanded and plotted over each image for
clarity.

Figure 8: Density and radial speed distributions from the WSA–ENLIL simulation
in the ecliptic plane at 01:01:06 UT on 2020 December 1 (top) and 07:00:21 UT
on 2020 December 2 (bottom). The values of the density and speed are scaled by
the color bar at the right corner of each panel. The positions of the Earth,
STEREO A, and PSP are marked by the circle, triangle, and square,
respectively. The first time indicates the moment when the CME have just
passed PSP and would arrive at STEREO A. The second time is close to the CME
arrival time at the Earth.

Figure 9: In situ solar wind measurements and simulation results at PSP. From
top to bottom, the panels show the normalized 314 eV electron PAD, proton
density, bulk speed, proton temperature, magnetic field strength, and
components, respectively. Note that there are no plasma data from PSP during
the time of interest. The red dashed lines represent the simulation results.
The blue vertical dashed lines indicate the observed arrival times of the
shocks marked with “S0” and “S”, respectively. The shaded regions indicate the
intervals of two ICMEs.

Figure 10: In situ solar wind measurements and simulation results at STEREO A.
From top to bottom, the panels show the normalized 93–247 eV electron PAD,
proton density, bulk speed, proton temperature, magnetic field strength, and
components, respectively. The red dashed lines represent the simulation
results, and the density from the simulation is divided by a factor of 10. The
dotted line in the fourth panel denotes the expected proton temperature
calculated from the observed speed (Lopez, 1987). The blue vertical dashed
line indicates the observed arrival time of the shock marked with “S”, and the
shaded region indicates the interval of the ICME.

Figure 11: In situ solar wind measurements and simulation results at Wind.
From top to bottom, the panels show the normalized 265 eV electron PAD, proton
density, bulk speed, proton temperature, magnetic field strength, and
components, respectively. The red dashed lines represent the simulation
results, and the density from the simulation is divided by a factor of 5. The
dotted line in the fourth panel denotes the expected proton temperature
calculated from the observed speed (Lopez, 1987). The black vertical dashed
lines indicate the interval of the ICME–like structure.

Table 1: The shock ellipsoid fitting parameters and corresponding times.

Time (UT) | $\theta$ | $\phi$ | $h$ | $a$ | $c$
---|---|---|---|---|---
STEREO A | GOES/SOHO | (°) | (°) | $R_{\sun}$ | $R_{\sun}$ | $R_{\sun}$
EUVI 12:55:30 | SUVI 12:56:37 | -95 | -18 | 1.31 | 0.51 | 0.51
COR1 13:10:18 | SUVI 13:10:37 | -90 | -15 | 1.79 | 1.27 | 1.06
COR1 13:15:18 | SUVI 13:15:27 | -90 | -15 | 2.00 | 1.52 | 1.35
COR1 13:20:18 | SUVI 13:20:37 | -85 | -15 | 2.34 | 2.12 | 1.83
COR2 13:24:00 | C2 13:23:38 | -80 | -15 | 2.73 | 3.10 | 2.38
COR2 13:39:00 | C3 13:39:54 | -80 | -23 | 3.16 | 5.49 | 4.77
COR2 13:54:00 | C3 13:51:55 | -80 | -25 | 3.52 | 7.65 | 6.51
– | C3 14:03:54 | -80 | -25 | 3.83 | 9.43 | 8.21
– | C3 14:15:38 | -80 | -27 | 4.10 | 11.34 | 9.84
COR2 14:24:00 | C3 14:27:39 | -80 | -30 | 4.41 | 13.10 | 11.53
COR2 14:39:00 | C3 14:39:38 | -80 | -30 | 4.68 | 14.85 | 13.16
– | C3 15:07:59 | -80 | -33 | 5.49 | 18.53 | 16.90
– | C3 15:18:34 | -80 | -35 | 5.71 | 19.88 | 18.24

Note. — $\theta$, $\phi$ and $h$ are the longitude, latitude, and height of
the center of the ellipsoid, respectively. $a$ and $c$ are the lengths of the
two semiprincipal axes. Because we assume the cross section of the shock
ellipsoid perpendicular to the propagation direction to be a circle, there are
only 5 free parameters in our fitting.

## References

  * Arge et al. (2004) Arge, C. N., Luhmann, J. G., Odstrcil, D., Schrijver, C. J., & Li, Y. 2004, Journal of Atmospheric and Solar-Terrestrial Physics, 66, 1295, doi: 10.1016/j.jastp.2004.03.018
  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., 204, 49, doi: 10.1007/s11214-016-0244-5
  * Biondo et al. (2021) Biondo, R., Pagano, P., Reale, F., & Bemporad, A. 2021, A&A, 654, L3, doi: 10.1051/0004-6361/202141892
  * Cheng et al. (2013) Cheng, X., Zhang, J., Ding, M. D., Liu, Y., & Poomvises, W. 2013, ApJ, 763, 43, doi: 10.1088/0004-637X/763/1/43
  * Cheng et al. (2012) Cheng, X., Zhang, J., Olmedo, O., et al. 2012, ApJ, 745, L5, doi: 10.1088/2041-8205/745/1/L5
  * Cheng et al. (2014) Cheng, X., Ding, M. D., Guo, Y., et al. 2014, ApJ, 780, 28, doi: 10.1088/0004-637X/780/1/28
  * de Lucas et al. (2011) de Lucas, A., Schwenn, R., dal Lago, A., Marsch, E., & Clúa de Gonzalez, A. L. 2011, Journal of Atmospheric and Solar-Terrestrial Physics, 73, 1281, doi: 10.1016/j.jastp.2010.12.011
  * Domingo et al. (1995) Domingo, V., Fleck, B., & Poland, A. I. 1995, Sol. Phys., 162, 1, doi: 10.1007/BF00733425
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7, doi: 10.1007/s11214-015-0211-6
  * Hess & Zhang (2014) Hess, P., & Zhang, J. 2014, ApJ, 792, 49, doi: 10.1088/0004-637X/792/1/49
  * Howard et al. (2008) Howard, R. A., Moses, J. D., Vourlidas, A., et al. 2008, Space Sci. Rev., 136, 67, doi: 10.1007/s11214-008-9341-4
  * Hu et al. (2016) Hu, H., Liu, Y. D., Wang, R., Möstl, C., & Yang, Z. 2016, ApJ, 829, 97, doi: 10.3847/0004-637X/829/2/97
  * Hu et al. (2017) Hu, H., Liu, Y. D., Wang, R., et al. 2017, ApJ, 840, 76, doi: 10.3847/1538-4357/aa6d54
  * Kaiser et al. (2008) Kaiser, M. L., Kucera, T. A., Davila, J. M., et al. 2008, Space Sci. Rev., 136, 5, doi: 10.1007/s11214-007-9277-0
  * Kollhoff et al. (2021) Kollhoff, A., Kouloumvakos, A., Lario, D., et al. 2021, A&A, 656, A20, doi: 10.1051/0004-6361/202140937
  * Kwon et al. (2014) Kwon, R.-Y., Zhang, J., & Olmedo, O. 2014, ApJ, 794, 148, doi: 10.1088/0004-637X/794/2/148
  * Kwon et al. (2015) Kwon, R.-Y., Zhang, J., & Vourlidas, A. 2015, ApJ, 799, L29, doi: 10.1088/2041-8205/799/2/L29
  * Lario et al. (2017) Lario, D., Kwon, R. Y., Richardson, I. G., et al. 2017, ApJ, 838
  * Lario et al. (2021) Lario, D., Richardson, I. G., Palmerio, E., et al. 2021, ApJ, 920, 123, doi: 10.3847/1538-4357/ac157f
  * Leblanc et al. (1998) Leblanc, Y., Dulk, G. A., & Bougeret, J.-L. 1998, Sol. Phys., 183, 165, doi: 10.1023/A:1005049730506
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17, doi: 10.1007/s11207-011-9776-8
  * Liu et al. (2010a) Liu, Y., Davies, J. A., Luhmann, J. G., et al. 2010a, ApJ, 710, L82, doi: 10.1088/2041-8205/710/1/L82
  * Liu et al. (2011) Liu, Y., Luhmann, J. G., Bale, S. D., & Lin, R. P. 2011, ApJ, 734, 84, doi: 10.1088/0004-637X/734/2/84
  * Liu et al. (2010b) Liu, Y., Thernisien, A., Luhmann, J. G., et al. 2010b, ApJ, 722, 1762, doi: 10.1088/0004-637X/722/2/1762
  * Liu et al. (2008) Liu, Y., Luhmann, J. G., Müller-Mellin, R., et al. 2008, ApJ, 689, 563, doi: 10.1086/592031
  * Liu et al. (2016) Liu, Y. D., Hu, H., Wang, C., et al. 2016, ApJS, 222, 23, doi: 10.3847/0067-0049/222/2/23
  * Liu et al. (2017) Liu, Y. D., Hu, H., Zhu, B., Luhmann, J. G., & Vourlidas, A. 2017, ApJ, 834, 158, doi: 10.3847/1538-4357/834/2/158
  * Liu et al. (2013) Liu, Y. D., Luhmann, J. G., Lugaz, N., et al. 2013, ApJ, 769, 45, doi: 10.1088/0004-637X/769/1/45
  * Liu et al. (2014a) Liu, Y. D., Yang, Z., Wang, R., et al. 2014a, ApJ, 793, L41, doi: 10.1088/2041-8205/793/2/L41
  * Liu et al. (2019) Liu, Y. D., Zhu, B., & Zhao, X. 2019, ApJ, 871, 8, doi: 10.3847/1538-4357/aaf425
  * Liu et al. (2012) Liu, Y. D., Luhmann, J. G., Möstl, C., et al. 2012, ApJ, 746, L15, doi: 10.1088/2041-8205/746/2/L15
  * Liu et al. (2014b) Liu, Y. D., Luhmann, J. G., Kajdič, P., et al. 2014b, Nature Communications, 5, 3481, doi: 10.1038/ncomms4481
  * Lopez (1987) Lopez, R. E. 1987, J. Geophys. Res., 92, 11189, doi: 10.1029/JA092iA10p11189
  * Lugaz et al. (2017) Lugaz, N., Temmer, M., Wang, Y., & Farrugia, C. J. 2017, Sol. Phys., 292, 64, doi: 10.1007/s11207-017-1091-6
  * Lugaz et al. (2009) Lugaz, N., Vourlidas, A., & Roussev, I. I. 2009, Annales Geophysicae, 27, 3479, doi: 10.5194/angeo-27-3479-2009
  * Manchester et al. (2017) Manchester, W., Kilpua, E. K. J., Liu, Y. D., et al. 2017, Space Sci. Rev., 212, 1159, doi: 10.1007/s11214-017-0394-0
  * Marubashi et al. (2015) Marubashi, K., Akiyama, S., Yashiro, S., et al. 2015, Sol. Phys., 290, 1371, doi: 10.1007/s11207-015-0681-4
  * Mishra et al. (2015) Mishra, W., Srivastava, N., & Chakrabarty, D. 2015, Sol. Phys., 290, 527, doi: 10.1007/s11207-014-0625-4
  * Nieves-Chinchilla et al. (2012) Nieves-Chinchilla, T., Colaninno, R., Vourlidas, A., et al. 2012, Journal of Geophysical Research (Space Physics), 117, A06106, doi: 10.1029/2011JA017243
  * Odstrcil et al. (2004) Odstrcil, D., Riley, P., & Zhao, X. P. 2004, Journal of Geophysical Research (Space Physics), 109, A02116, doi: 10.1029/2003JA010135
  * Palmerio et al. (2021) Palmerio, E., Kilpua, E. K. J., Witasse, O., et al. 2021, Space Weather, 19, e2020SW002654, doi: 10.1029/2020SW002654
  * Patsourakos & Vourlidas (2012) Patsourakos, S., & Vourlidas, A. 2012, Sol. Phys., 281, 187, doi: 10.1007/s11207-012-9988-6
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Reinard et al. (2012) Reinard, A. A., Lynch, B. J., & Mulligan, T. 2012, ApJ, 761, 175, doi: 10.1088/0004-637X/761/2/175
  * Reiner et al. (2001) Reiner, M. J., Kaiser, M. L., & Bougeret, J. L. 2001, J. Geophys. Res., 106, 29989, doi: 10.1029/2000JA002228
  * Rouillard et al. (2011) Rouillard, A. P., Odstřcil, D., Sheeley, N. R., et al. 2011, ApJ, 735, 7, doi: 10.1088/0004-637X/735/1/7
  * Savani et al. (2015) Savani, N. P., Vourlidas, A., Szabo, A., et al. 2015, Space Weather, 13, 374, doi: 10.1002/2015SW001171
  * Syed Ibrahim et al. (2019) Syed Ibrahim, M., Joshi, B., Cho, K. S., Kim, R. S., & Moon, Y. J. 2019, Sol. Phys., 294, 54, doi: 10.1007/s11207-019-1443-5
  * Thernisien et al. (2009) Thernisien, A., Vourlidas, A., & Howard, R. A. 2009, Sol. Phys., 256, 111, doi: 10.1007/s11207-009-9346-5
  * Thernisien et al. (2006) Thernisien, A. F. R., Howard, R. A., & Vourlidas, A. 2006, ApJ, 652, 763, doi: 10.1086/508254
  * Webb et al. (2013) Webb, D. F., Möstl, C., Jackson, B. V., et al. 2013, Sol. Phys., 285, 317, doi: 10.1007/s11207-013-0260-5
  * Xie et al. (2021) Xie, H., Gopalswamy, N., & Akiyama, S. 2021, ApJ, 922, 64, doi: 10.3847/1538-4357/ac23cc
  * Zhang et al. (2012) Zhang, J., Cheng, X., & Ding, M.-D. 2012, Nature Communications, 3, 747, doi: 10.1038/ncomms1753
  * Zhu et al. (2021) Zhu, B., Liu, Y. D., Kwon, R.-Y., et al. 2021, ApJ, 921, 26, doi: 10.3847/1538-4357/ac106b
  * Zurbuchen & Richardson (2006) Zurbuchen, T. H., & Richardson, I. G. 2006, Space Sci. Rev., 123, 31, doi: 10.1007/s11214-006-9010-4