# Similarity of multiple dose-response curves in interlaboratory studies in
regulatory toxicology

Ludwig A. Hothorn,
Im Grund 12, D-31867 Lauenau, Germany
(retired from Leibniz University Hannover)
Dimitrios Spiliotopoulos,
Xenometrix AG, CH-4123 Allschwil, Switzerland

###### Abstract

To claim similarity of multiple dose-response curves in interlaboratory
studies in regulatory toxicology is a relevant issue during the assay
validation process. Here we demonstrated the use of dose-by-laboratory
interaction contrasts, particularly Williams-type by total mean contrasts.
With the CRAN packages statint and multcomp in the open-source software R the
estimation of adjusted p-values or compatible simultaneous confidence
intervals is relatively easy. The interpretation in terms of global or partial
equivalence, i.e. similarity, is challenging, because thresholds are not
available a-priori. This approach is demonstrated by selected in-vitro Ames
MPF assay data.

##  1 Introduction

A relevant objective of interlaboratory studies in regulatory toxicology is to
demonstrate the similarity of multiple dose-response curves in $l$
participating laboratories $Lab_{l}$. Almost all of these bioassays are based
on the design $[C-,D_{1},...,D_{k}]$, assuming a primary normally distributed
endpoint. This provides a comparison of dose-response curves, whereby the dose
can be modeled as a qualitative factor using contrast tests or as a
quantitative covariate using nonlinear regression models. Both approaches have
advantages and disadvantages in this specific design (comparisons to C- are
meaningful, $k$ is low, e.g. 3), see for details e.g. [9]. Similarity is
statistically translated into equivalence, which is demonstrated by inclusion
in $(1-2\alpha)$ confidence intervals [1]. The crux of the matter is that the
required tolerability limits are a-priori endpoint/assay specific not defined.
Either a post-hoc comparative interpretation of the estimated intervals with
regard to their tolerability is carried out or, alternatively, the proof-of-
hazard standard used in toxicology is used as an equivalent, where $p>0.10$.
One can demonstrate equivalence for a nonlinear function and thus the whole
curve, or only a part of the curve (e.g. the linearized one) or for a single
estimator, such as benchmark dose (BMD), LD50 or no-observed-adverse-effect-
concentration (NOAEC).
Here we are guided by the evidence of a significant trend as an effectiveness
criterion in the guidelines, e.g. for in-vivo micronucleus assay [17]. And as
a trend test we choose the Williams test ([26]) because it is geared to
comparisons with C-, as a multiple contrast test it is robust for several
forms of dose-response dependence, modeling the dose as a qualitative factor
level and recommended in guidance, e.g. US-NTP for continuous endpoints. If
the dose factor is modeled together with the laboratory factor, a lack of two-
way interaction is a criterion for similarity. However, this global criterion
cannot indicate that at least one, any interaction exists, but it cannot
indicate which one exactly. It could be that 1 out of 10 laboratories behave
differently (but the other 9 are similar). Or only the difference between
$D_{3}-D_{2}$ is different in at least a single lab, which is toxicologically
not relevant. Furthermore, confidence intervals are not readily available for
the F-test statistics. A better alternative are interaction contrasts. For the
factor dose the Williams contrast is described here, for the factor laboratory
the total mean contrast (to leave the number of comparisons at $k$); other
definitions are possible.

To characterize the Ames fluctuation test [19] used the lowest ineffective
concentration as criterion, the fold change increase to compare standard
miniaturized and Ames II and MPF Assay [4, 21], where the confidence intervals
for benchmark dose estimates were used for the blood Pig-a gene mutation assay
[3] (requiring the same underlying nonlinear model for any condition).

##  2 A motivating example: Ames MPF in-vitro assay

In an extension interlaboratory study, the new MPF-assay was considered for 6
selected compounds, in 7 particular laboratories, using 5 strains and 2 kinds
of metabolization [22]. Substance 4 was used as data example, where $n_{i}=6$
(pooled over $n_{i}=3$ per metabolization S9+, S9-) per concentration and
laboratory were available. The number of revertants were transformed to
Nishiyama [16] to achieve approximate normal distribution in this small sample
design. The boxplots in Figure 1 reveal monotonic increasing curves with
similar shapes in the different labs, still additive shifts between the labs
(e.g. all comparable values in lab 3 are less than in lab 7), tied values,
variance heterogeneity (commonly higher variance in higher concentrations, but
not always), and small sample sizes per concentration and laboratory (pooled
over $n_{i}=3$ for S9+, S9-).

Figure 1: Boxplot for concentration-by-laboratory interaction in compound 4

##  3 Williams-by-total mean interaction contrasts

In the following, the interaction contrasts are derived for the two-way layout
dose-by-laboratory. First, the Williams-contrast test is derived for the
primary factor ’dose’, then the contrast test for comparison with the total
mean for the secondary factor ’laboratories’, and finally the resulting
interaction contrast.

###  3.1 Factor dose: The Williams-type multiple contrast test

Williams original procedure [26] is a rather complex approach based on maximum
likelihood estimators under total order restriction, hard to generalize. Here
we use its re-formulation as multiple contrast test [2], assuming order
restriction with respect to C-. A multiple contrast test is a maximum test on
$q$ elementary tests (where $q$ depends on the kind of multiple contrast)
(i.e. a special version of an union-intersection test):
$t_{MCT}=max(t_{1},...,t_{q})$ which follows jointly
$(t_{1},\ldots,t_{q})^{\prime}$ a $q$-variate $t$\- distribution with a common
degree of freedom $df$ and the correlation matrix $R$, which is depending on
$c_{i},n_{i}$ under simple assumptions (or more complex for other MCT’s [9]).
The included single contrast tests
$t_{SC}=\sum_{i=0}^{k}c_{i}\bar{x}_{i}/S\sqrt{\sum_{i}^{k}c_{i}^{2}/n_{i}}$
and $\sum_{i=0}^{k}c_{i}=0$ differ by their weights $c_{i}$, i.e. particular
contrast matrix definitions determine the particular test version. For the
simple balanced design with k=3 the contrast matrix for the one-sided
Williams-type test is:

$c_{i}$ | C | $D_{1}$ | $D_{2}$ | $D_{3}$
---|---|---|---|---
$c_{a}$ | -1 | 0 | 0 | 1
$c_{b}$ | -1 | 0 | 1/2 | 1/2
$c_{c}$ | -1 | 1/3 | 1/3 | 1/3

The simple Williams contrast illustrates the idea and interoperability of an
MCT: either the comparison of $D_{3}$ vs. C is significant (i.e. strict
monotone), or to the pooled $(D_{3}+D_{2})/2$, or $(D_{3}+D_{2}+D_{1})/3$,
i.e. a plateau (or all three, or neither ($H_{0}$). Either multiplicity-
adjusted p-values or simultaneous confidence limits (here one-sided lower
limits) are available:
$[\sum_{i=0}^{k}c_{i}\bar{x}_{i}-S*t_{q,df,R,2-sided,1-\alpha}\sqrt{\sum_{i}^{k}c_{i}^{2}/n_{i}}]$.

This makes the Williams test the recommended trend test in
pharmacology/toxicology: sensitive to some monotonic and partially non-
monotonic forms, a comparison to the control, easily interpretable confidence
intervals for the required effect sizes. In regulatory toxicology not only the
difference of means is relevant as effect size [23], such as proportions or
counts. Therefore, the Williams test is available for difference [10], risk
ratio or odds ratio of proportions [9], ratio-to-control estimates [11], the
nonparametric relative effect sizes [14], hazard rates [7], multiple endpoints
[6], heteroscedastic error terms [8] and poly-k-adjusted tumor rates [20]).

###  3.2 Factor laboratory: The total mean multiple contrast test

Qualitative levels of a factor (without order restriction) are commonly
compared by Tukey’s all-pairs comparison procedure [25]. Their numerous, exact
$(l*(l-1)/2)*2$ comparisons make this procedure conservative and especially
difficult to interpret. One or more deviating laboratories can also be
identified by total mean comparisons [18], but with only $l$ to be
interpreted. For $l=4$ laboratories the contrast matrix is (in a balanced
design) simply

$c_{i}$ | $Lab_{A}$ | $Lab_{B}$ | $Lab_{C}$ | $Lab_{D}$
---|---|---|---|---
$c_{a}$ | -1 | 1/3 | 1/3 | 1/3
$c_{b}$ | 1/3 | -1 | 1/3 | 1/3
$c_{c}$ | 1/3 | 1/3 | -1 | 1/3
$c_{d}$ | 1/3 | 1/3 | 1/3 | -1

###  3.3 Interaction dose-by-laboratory: Williams-by-total mean interaction
contrast test

When reformulating the common ANOVA interaction model:
$y_{klj}=\mu+\alpha_{k}+\beta_{l}+(\alpha\beta)_{kl}+\epsilon_{klj}$ (with $j$
replicates) into a cell means model $y_{klj}=\mu_{kl}+\epsilon_{klj}$ the null
hypothesis of interaction can be written:
$H_{0}:(\mu_{kl}-\mu_{kl^{\prime}})-(\mu_{k^{\prime}l}-\mu_{k^{\prime}l^{\prime}})=0$
for all $(k,k^{\prime})$ and $(l,l^{\prime})$ [13] as a special form of
product-type interactions [5]. Using the both contrast matrices for Williams-
type between doses $C_{Dose}$ and between laboratories total mean $C_{l}ab$,
the interaction contrast is the Kronecker product
$C_{Dose,lab}=C_{Dose}\otimes C_{lab}$. This can easily be realized for
different two-way contrast types with the CRAN-package statint [24]. Explicit
for the above Ames assay example Table 1 contains the 36 individual contrasts
for Williams-by-total mean interaction with k=7
$([0,0.015625,0.03125,0.0625,0.125,0.25,0.5]$ and l=6
$[Lab_{1},Lab_{3},Lab_{4},Lab_{5},Lab_{6},Lab_{7}])$:

No. | IA-contrast lab:dose abbreviation | Contrast type
---|---|---
1 | ((1 - 2,3,4,5,6,7):0.5) - ((1 - 2,3,4,5,6,7):0) |  (Lab1 vs. all Labs) by Williams contrast 1 $(0.5-C-)$
2 | ((1 - 2,3,4,5,6,7):0.25,0.5) - ((1 - 2,3,4,5,6,7):0) |  (Lab1 vs. all Labs) by Williams contrast 2 $((0.5+0.25)/2-C-)$
3 | ((1 - 2,3,4,5,6,7):0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) |  (Lab1 vs. all Labs) by Williams contrast 3 $((0.5+0.25+0.125)/3-C-)$
4 | ((1 - 2,3,4,5,6,7):0.0625,0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) |  (Lab1 vs. all Labs) by Williams contrast 4 $((0.5+0.25+0.125+0.0625)/4-C-)$
5 | ((1 - 2,3,4,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) |  (Lab1 vs. all Labs) by Williams contrasts 5 $((0.5+0.25+0.125+0.0625+0.031)/5-C-)$
6 | ((1 - 2,3,4,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) |  (Lab1 vs. all Labs) by Williams contrast 6 $((0.5+0.25+0.125+0.0625+0.031+0.015)/6-C-)$
7 | ((2 - 1,3,4,5,6,7):0.5) - ((2 - 1,3,4,5,6,7):0) | (Lab2 vs. all Labs) by Williams contrast 1
8 | ((2 - 1,3,4,5,6,7):0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | (Lab2 vs. all Labs) by Williams contrast 2
9 | ((2 - 1,3,4,5,6,7):0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | (Lab2 vs. all Labs) by Williams contrast 3
10 | ((2 - 1,3,4,5,6,7):0.0625,0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | (Lab2 vs. all Labs) by Williams contrast 4
11 | ((2 - 1,3,4,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | (Lab2 vs. all Labs) by Williams contrast 5
12 | ((2 - 1,3,4,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | (Lab2 vs. all Labs) by Williams contrast 6
13 | ((3 - 1,2,4,5,6,7):0.5) - ((3 - 1,2,4,5,6,7):0) | (Lab3 vs. all Labs) by Williams contrast 1
14 | ((3 - 1,2,4,5,6,7):0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | (Lab3 vs. all Labs) by Williams contrast 2
15 | ((3 - 1,2,4,5,6,7):0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | (Lab3 vs. all Labs) by Williams contrast 3
16 | ((3 - 1,2,4,5,6,7):0.0625,0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | (Lab3 vs. all Labs) by Williams contrast 4
17 | ((3 - 1,2,4,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | (Lab3 vs. all Labs) by Williams contrast 5
18 | ((3 - 1,2,4,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | (Lab3 vs. all Labs) by Williams contrast 6
19 | ((4 - 1,2,3,5,6,7):0.5) - ((4 - 1,2,3,5,6,7):0) | (Lab4 vs. all Labs) by Williams contrast 1
20 | ((4 - 1,2,3,5,6,7):0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | (Lab4 vs. all Labs) by Williams contrast 2
21 | ((4 - 1,2,3,5,6,7):0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | (Lab4 vs. all Labs) by Williams contrast 3
22 | ((4 - 1,2,3,5,6,7):0.0625,0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | (Lab4 vs. all Labs) by Williams contrast 4
23 | ((4 - 1,2,3,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | (Lab4 vs. all Labs) by Williams contrast 5
24 | ((4 - 1,2,3,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | (Lab4 vs. all Labs) by Williams contrast 6
25 | ((5 - 1,2,3,4,6,7):0.5) - ((5 - 1,2,3,4,6,7):0) | (Lab5 vs. all Labs) by Williams contrast 1
26 | ((5 - 1,2,3,4,6,7):0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | (Lab5 vs. all Labs) by Williams contrast 2
27 | ((5 - 1,2,3,4,6,7):0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | (Lab5 vs. all Labs) by Williams contrast 3
28 | ((5 - 1,2,3,4,6,7):0.0625,0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | (Lab5 vs. all Labs) by Williams contrast 4
29 | ((5 - 1,2,3,4,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | (Lab5 vs. all Labs) by Williams contrast 5
30 | ((5 - 1,2,3,4,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | (Lab5 vs. all Labs) by Williams contrast 6
31 | ((6 - 1,2,3,4,5,7):0.5) - ((6 - 1,2,3,4,5,7):0) | (Lab6 vs. all Labs) by Williams contrast 1
32 | ((6 - 1,2,3,4,5,7):0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | (Lab6 vs. all Labs) by Williams contrast 2
33 | ((6 - 1,2,3,4,5,7):0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | (Lab6 vs. all Labs) by Williams contrast 3
34 | ((6 - 1,2,3,4,5,7):0.0625,0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | (Lab6 vs. all Labs) by Williams contrast 4
35 | ((6 - 1,2,3,4,5,7):0.03125,0.0625,0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | (Lab6 vs. all Labs) by Williams contrast 5
36 | ((6 - 1,2,3,4,5,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | (Lab6 vs. all Labs) by Williams contrast 6
37 | ((7 - 1,2,3,4,5,6):0.5) - ((7 - 1,2,3,4,5,6):0) | (Lab7 vs. all Labs) by Williams contrast 1
38 | ((7 - 1,2,3,4,5,6):0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | (Lab7 vs. all Labs) by Williams contrast 2
39 | ((7 - 1,2,3,4,5,6):0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | (Lab7 vs. all Labs) by Williams contrast 3
40 | ((7 - 1,2,3,4,5,6):0.0625,0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | (Lab7 vs. all Labs) by Williams contrast 4
41 | ((7 - 1,2,3,4,5,6):0.03125,0.0625,0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | (Lab7 vs. all Labs) by Williams contrast 5
42 | ((7 - 1,2,3,4,5,6):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | (Lab7 vs. all Labs) by Williams contrast 6
Table 1: Interaction contrast structure for Williams trend-by-total means

This table of this already ample example makes clear the advantages and
disadvantages of this approach: on the one hand, 36 individual decisions
instead of one single F-test, on the other hand, the individual decisions
allow the statement between which laboratories and which trend contrasts there
is similarity and where not. In terms of an union-intersection test (UIT), the
approach is conservative, but not so extreme because of the high correlation
between the contrasts. If you want to know which interactions are equivalent
and which are not, use a combined intersection-union-by-union-intersection
test (IUT-UIT): IUT for the equivalence tests in both directions, and UIT
between the IA contrasts. It is also possible to achieve a global statement by
using an IUT-IUT (then the IA-test at the marginal $\alpha=0.05$ level), where
all individual tests should fulfill $p>0.10$.

###  3.4 A possible simplification: $(D_{k}-C-)$-by-total mean interaction
contrasts

Assuming a monotonic dose-response dependence, the entire William contrast can
be replaced by $(D_{k}-C-)$ comparison in a simplified way. This is
surprising, but $(D_{k}-C-)$ is the most important single contrast in the
William test, is part of the corresponding closure test [9] and part of the
strict trend test according [15]. This reduction especially simplifies the
interpretation to only $k$ interaction contrasts. This simplification should
be avoided if downturns at higher dose(s) are possible or if you are more
interested in comparing the NOAEC instead of a trend itself. The
simplification becomes clear looking on Table 2 for the above Ames assay
example: only 7 contrasts, claiming for local equivalence for laboratories
$[1,2,3,4,5,6]$ and but not for global equivalence because laboratory 7
behaves borderline non-equivalent.

Number | Interaction contrast | p-value
---|---|---
1 | ((1 - 2,3,4,5,6,7):0.5) - ((1 - 2,3,4,5,6,7):0) | 0.985
2 | ((2 - 1,3,4,5,6,7):0.5) - ((2 - 1,3,4,5,6,7):0) | 0.410
3 | ((3 - 1,2,4,5,6,7):0.5) - ((3 - 1,2,4,5,6,7):0) | 0.999
4 | ((4 - 1,2,3,5,6,7):0.5) - ((4 - 1,2,3,5,6,7):0) | 0.999
5 | ((5 - 1,2,3,4,6,7):0.5) - ((5 - 1,2,3,4,6,7):0) | 0.584
6 | ((6 - 1,2,3,4,5,7):0.5) - ((6 - 1,2,3,4,5,7):0) | 0.120
7 | ((7 - 1,2,3,4,5,6):0.5) - ((7 - 1,2,3,4,5,6):0) | 0.100
Table 2: Compound 4: interaction contrast for C vs. $D_{k}$ by grand means
between labs- S9+

###  3.5 A modification

Variance heterogeneity often occurs in these assays and their ignorance can
lead to significant bias. Therefore you should use the sandwich variance
estimator (see below using vcov=sandwich argument in the function glht()) [8].

##  4 Evaluation of the example using the CRAN packages statint and multcomp

The following R-code can be used for the Ames assay example for the data
object A984p:

    library(multcomp); library(statint); library(sandwich)
    InteractionContrastsP984t <- iacontrast(fa=A984p$lab,fb=A984p$Conc,
                                           typea="totalMean", typeb="Williams")
    A984p$labConct <- InteractionContrastsP984t$fab
    CellMeansModelP984t <- lm(Trans ~ labConct-1, data=A984p) # cell means model
    MultTestP984t <- glht(model=CellMeansModelP984t,
                     linfct = mcp(labConct=InteractionContrastsP984t$cmab), vcov=sandwich)
    piap984t<-summary(MultTestP984t)ΨΨΨΨΨΨ#calculating adjusted p-values

Number | Interaction contrast | p-value
---|---|---
1 | ((1 - 2,3,4,5,6,7):0.5) - ((1 - 2,3,4,5,6,7):0) | 0.9994
2 | ((1 - 2,3,4,5,6,7):0.25,0.5) - ((1 - 2,3,4,5,6,7):0) | 0.999
3 | ((1 - 2,3,4,5,6,7):0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) | 0.999
4 | ((1 - 2,3,4,5,6,7):0.0625,0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) | 0.9698
5 | ((1 - 2,3,4,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) | 0.8706
6 | ((1 - 2,3,4,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((1 - 2,3,4,5,6,7):0) | 0.9230
7 | ((2 - 1,3,4,5,6,7):0.5) - ((2 - 1,3,4,5,6,7):0) | 0.5901
8 | ((2 - 1,3,4,5,6,7):0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | 0.4831
9 | ((2 - 1,3,4,5,6,7):0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | 0.8529
10 | ((2 - 1,3,4,5,6,7):0.0625,0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | 0.8546
11 | ((2 - 1,3,4,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | 0.9306
12 | ((2 - 1,3,4,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((2 - 1,3,4,5,6,7):0) | 0.9719
13 | ((3 - 1,2,4,5,6,7):0.5) - ((3 - 1,2,4,5,6,7):0) | 0.999
14 | ((3 - 1,2,4,5,6,7):0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | 0.999
15 | ((3 - 1,2,4,5,6,7):0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | 0.9990
16 | ((3 - 1,2,4,5,6,7):0.0625,0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | 0.9887
17 | ((3 - 1,2,4,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | 0.999
18 | ((3 - 1,2,4,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((3 - 1,2,4,5,6,7):0) | 0.999
19 | ((4 - 1,2,3,5,6,7):0.5) - ((4 - 1,2,3,5,6,7):0) | 0.999
20 | ((4 - 1,2,3,5,6,7):0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | 0.999
21 | ((4 - 1,2,3,5,6,7):0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | 0.999
22 | ((4 - 1,2,3,5,6,7):0.0625,0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | 0.999
23 | ((4 - 1,2,3,5,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | 0.999
24 | ((4 - 1,2,3,5,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((4 - 1,2,3,5,6,7):0) | 0.999
25 | ((5 - 1,2,3,4,6,7):0.5) - ((5 - 1,2,3,4,6,7):0) | 0.7772
26 | ((5 - 1,2,3,4,6,7):0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | 0.8020
27 | ((5 - 1,2,3,4,6,7):0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | 0.9052
28 | ((5 - 1,2,3,4,6,7):0.0625,0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | 0.9928
29 | ((5 - 1,2,3,4,6,7):0.03125,0.0625,0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | 0.9828
30 | ((5 - 1,2,3,4,6,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((5 - 1,2,3,4,6,7):0) | 0.9999
31 | ((6 - 1,2,3,4,5,7):0.5) - ((6 - 1,2,3,4,5,7):0) | 0.1843
32 | ((6 - 1,2,3,4,5,7):0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | 0.2064
33 | ((6 - 1,2,3,4,5,7):0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | 0.7596
34 | ((6 - 1,2,3,4,5,7):0.0625,0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | 0.8450
35 | ((6 - 1,2,3,4,5,7):0.03125,0.0625,0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | 0.9857
36 | ((6 - 1,2,3,4,5,7):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((6 - 1,2,3,4,5,7):0) | 0.9952
37 | ((7 - 1,2,3,4,5,6):0.5) - ((7 - 1,2,3,4,5,6):0) | 0.1511
38 | ((7 - 1,2,3,4,5,6):0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | 0.4369
39 | ((7 - 1,2,3,4,5,6):0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | 0.3633
40 | ((7 - 1,2,3,4,5,6):0.0625,0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | 0.4856
41 | ((7 - 1,2,3,4,5,6):0.03125,0.0625,0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | 0.7422
42 | ((7 - 1,2,3,4,5,6):0.015625,0.03125,0.0625,0.125,0.25,0.5) - ((7 - 1,2,3,4,5,6):0) | 0.8666
Table 3: Compound 4: interaction contrast for Williams trend-by grand means
between labs- S9+

The related adjusted p-values for this IUT-UIT test in Table 3 are
predominantly large, which indicate equivalence, even global equivalence. Only
one contrast, No. 37 shows with $p=0.1511$ a tendency of non-equivalence; but
the single comparison $(D_{k}-C-)$ should not be relevant for the Ames assay
because the decision is based on the lower doses. This p-value of $0.1511$ is
part of the comparisons of Lab 7 against all others, which are not as large as
the others (but still equivalent). Notice, Lab 7 used a different source of
bacteria revealing a higher spontaneous reversion rate.

Even more informative is the plot of the simultaneous two-sided $(1-2\alpha)$
confidence intervals in Figure 2, demonstrating the problems with Lab 7 in the
scale of transformed revertant count differences. Notice, using the IUT-IUT
approach with marginal levels, global equivalence can be claimed for this
selected example.

Figure 2: 42 interaction contrasts: compatible simultaneous confidence
intervals

Notice, the (two-sided) p-value for the common F-test on dose-by-laboratory
interaction is $0.10$ indicating that an interaction cannot be excluded.
However, it is not clear for which laboratories and which dose comparisons - a
nice proof of the usefulness of the above IA-contrast method in this example.

##  5 Summary

Similarity of multiple dose-response curves in interlaboratory studies in
regulatory toxicology can be demonstrated using dose-by-laboratory interaction
contrasts [12, 9]. To illustrate a trend of the dose-response curve, Williams-
by-Laboratory interaction contrasts are proposed here. With help of the CRAN
packages statint and multcomp the estimation of adjusted p-values or
compatible simultaneous confidence intervals is relatively easy. The
interpretation in terms of global or partial equivalence (similarity) is
challenging, but impressive.

## References

  * [1] P. Bauer and M. Kieser.  A unifying approach for confidence intervals and testing of equivalence and difference.  Biometrika 83: 934 – 937, 1996.
  * [2] F. Bretz.  An extension of the Williams trend test to general unbalanced linear models.  Comput. Statist. and Data Analysis, 50(7):1735–1748, April 2006.
  * [3] S.D. Dertinger, et al.  Intra- and inter-laboratory reproducibility of the rat blood Pig-a gene mutation assay.  Environ. Mol. Mutagen. 61:500–507, 2020.
  * [4] I. Flueckiger-Isler and M. Kamber.  Direct comparison of the Ames microplate format (MPF) test in liquid medium with the standard Ames pre-incubation assay on agar plates by use of equivocal to weakly positive test compounds.  Mutat. Res./Genet. Toxic. and Environ. Mutagen. 747: 36–45, 2012.
  * [5] K.R. Gabriel and J. Putter and Y. Wax.  Simultaneous confidence intervals for product-type interaction contrasts.  J.R.S.S.-B 35: 234 – 244, 1973.
  * [6] M. Hasler and L.A. Hothorn.  A multivariate Williams-type trend procedure.  Statist. in Biopharm. Res., 4(1):57–65, 2012.
  * [7] E. Herberich and L.A. Hothorn.  Statistical evaluation of mortality in long-term carcinogenicity bioassays using a Williams-type procedure.  Regul. Toxicol. Pharm., 64:26–34, 2012.
  * [8] E. Herberich, J. Sikorski, and T. Hothorn.  A robust procedure for comparing multiple means under heteroscedasticity in unbalanced designs.  PLOS One, 5(3):e9788, March 2010.
  * [9] L.A. Hothorn.  Statistics in Toxicology- using R.  Chapman Hall, 2016.
  * [10] L.A. Hothorn and M. Sill and F. Schaarschmidt.  Evaluation of incidence rates in pre-clinical studies using a Williams-type procedure.  The Int. J. of Biostat. 6: 1,15, 2010.
  * [11] L.A. Hothorn and G.D. Djira.  A ratio-to-control Williams-type test for trend.  Pharm. Statist., 10(4):289–292, July 2011.
  * [12] L.A. Hothorn  Statistics of interlaboratory in vitro toxicological studies.  ATLA.31: 43–63.2003.
  * [13] A. Kitsche and F. Schaarschmidt.  Analysis of statistical interactions in factorial experiments.  J. Agr. Crop Sci. 201: 69 –79, 2014.
  * [14] F. Konietschke and L.A. Hothorn.  Evaluation of toxicological studies using a nonparametric Shirley-type trend test for comparing several dose levels with a control group.  Statist.in Biopharm. Res., 4(1):14–27, 2012.
  * [15] K.K. Lin and M.A. Rahman.  Comparisons of false negative rates from a trend test alone and from a trend test jointly with a control-high groups pairwise test in the determination of the carcinogenicity of new drugs.  J. Biopharm. Stat., 29(1):128–142, 2019.
  * [16] H. Nishiyama and T. Omori and I. Yoshimura.  A composite statistical procedure for evaluating genotoxicity using cell transformation assay data.  Environmetrics 14: 183 –192, 2003.
  * [17] OECD Guideline for testing of chemicals No. 474. In vivo micronucleus test, 2006.
  * [18] P. Pallmann and L.A. Hothorn.  Analysis of means (ANOM): a generalized approach using R.  J. Appl. Stat. doi:10.1080/02664763.2015.1117584, 2016.
  * [19] G. Reifferscheid and H.M. Maes et al.  International Round-Robin Study on the Ames Fluctuation Test.  Environm. and Mol. Mutagen. 53:185–197,2012.
  * [20] F. Schaarschmidt and M. Sill and L.A. Hothorn.  Poly-k-tend tests for survival adjusted analysis of tumor rates formulated as approximate multiple contrast test.  J. Biopharm. Stat. 18: 934 –948, 2008.
  * [21] D. Spiliotopoulos and C. Koelbert  Assessment of the miniaturized liquid Ames microplate format (MPF) for a selection of the test items from the recommended list of genotoxic and non-genotoxic chemicals.  Mutat. Res. 856–857:503218,2020.
  * [22] D. Spiliotopoulos et al.  An interlaboratory study to validate Ames MPF assay.  Mut. Res., to appear.
  * [23] E. Szoecs and R.B. Schaefer.  A comparison of statistical approaches for analysis of count and proportion data in ecotoxicology.  Environ. Sci. Pollut. Res. 22:13990–13999,2015.
  * [24] A. Kitsche and F. Schaarschmidt.  Generating labeled product type interaction contrasts.  CRAN package
  * [25] J. Tukey.  The problem of multiple comparisons.  Princeton University, 1953.
  * [26] D.A. Williams.  A test for differences between treatment means when several dose levels are compared with a zero dose control.  Biometrics, 27(1):103–117, 1971.


# Excitation of quasi-periodic fast-propagating waves in the early stage of
the solar eruption

Jialiang Hu Yunnan Observatories, Chinese Academy of Sciences, P.O. Box 110,
Kunming, Yunnan 650216, People’s Republic of China  University of Chinese
Academy of Sciences, Beijing 100049, People’s Republic of China  Yunnan Key
Laboratory of Solar Physics and Space Science, Kunming, Yunnan 650216,
People’s Republic of China  Jing Ye Yunnan Observatories, Chinese Academy of
Sciences, P.O. Box 110, Kunming, Yunnan 650216, People’s Republic of China
Yunnan Key Laboratory of Solar Physics and Space Science, Kunming, Yunnan
650216, People’s Republic of China  Yunnan Province China-Malaysia HF-VHF
Advanced Radio Astronomy Technology International Joint Laboratory, Kunming,
Yunnan 650216, People’s Republic of China Jing Ye, Zhixing Mei<EMAIL_ADDRESS><EMAIL_ADDRESS>Yuhao Chen Yunnan Observatories, Chinese Academy of
Sciences, P.O. Box 110, Kunming, Yunnan 650216, People’s Republic of China
University of Chinese Academy of Sciences, Beijing 100049, People’s Republic
of China  Yunnan Key Laboratory of Solar Physics and Space Science, Kunming,
Yunnan 650216, People’s Republic of China  Zhixing Mei Yunnan Observatories,
Chinese Academy of Sciences, P.O. Box 110, Kunming, Yunnan 650216, People’s
Republic of China  Yunnan Key Laboratory of Solar Physics and Space Science,
Kunming, Yunnan 650216, People’s Republic of China  Zehao Tang Yunnan
Observatories, Chinese Academy of Sciences, P.O. Box 110, Kunming, Yunnan
650216, People’s Republic of China  University of Chinese Academy of Sciences,
Beijing 100049, People’s Republic of China  Yunnan Key Laboratory of Solar
Physics and Space Science, Kunming, Yunnan 650216, People’s Republic of China
Jun Lin Yunnan Observatories, Chinese Academy of Sciences, P.O. Box 110,
Kunming, Yunnan 650216, People’s Republic of China  University of Chinese
Academy of Sciences, Beijing 100049, People’s Republic of China  Yunnan Key
Laboratory of Solar Physics and Space Science, Kunming, Yunnan 650216,
People’s Republic of China  Center for Astronomical Mega-Science, Chinese
Academy of Sciences, Beijing 100012, People’s Republic of China

###### Abstract

We propose a mechanism for the excitation of large-scale quasi-periodic fast-
propagating magnetoacoustic (QFP) waves observed on both sides of the coronal
mass ejection (CME). Through a series of numerical experiments, we
successfully simulated the quasi-static evolution of the equilibrium locations
of the magnetic flux rope in response to the change of the background magnetic
field, as well as the consequent loss of the equilibrium that eventually gives
rise to the eruption. During the eruption, we identified QFP waves propagating
radially outwards the flux rope, and tracing their origin reveals that they
result from the disturbance within the flux rope. Acting as an imperfect
waveguide, the flux rope allows the internal disturbance to escape to the
outside successively via its surface, invoking the observed QFP waves.
Furthermore, we synthesized the images of QFP waves on the basis of the data
given by our simulations, and found the consistence with observations. This
indicates that the leakage of the disturbance outside the flux rope could be a
reasonable mechanism of QFP waves.

waveguide, flux rope, QFP waves

##  1 Introduction

Solar atmosphere is a highly structured medium capable of supporting various
types of waves, which might belong to the families of the Magnetohydrodynamic
(MHD) waves in the solar corona including fast-mode (Williams et al., 2001,
2002; Nakariakov et al., 2003) and slow-mode waves (De Moortel et al., 2000;
Wang et al., 2009; Ofman & Wang, 2002) , and they have been extensively
investigated by many authors (e.g., see also Li et al. 2020; Nakariakov &
Kolotkov 2020; Banerjee et al. 2021; Nakariakov et al. 2021; Wang et al.
2021b). They impose additional constraints on the explanation or understanding
of fundamental physical processes occurring in solar activity, such as
particle acceleration and magnetic energy release (McLaughlin et al., 2018).
In addition, these perturbation may interact with the plasma and the magnetic
filed nearby, delivering the energy to the background atmosphere, which may
eventually compensate the radiative losses of thermal energy in the atmosphere
(Porter et al., 1994; Van Doorsselaere et al., 2020).

Global large scale EUV waves, usually observed in the EUV wavelengths, are
considered a ubiquitous phenomenon (Thompson et al., 1998; Wang et al., 2009;
Lin et al., 2010) associated with the disturbance to the atmosphere by the
eruption. They typically manifest as bright, ring-like structures that
gradually propagate outward, sweeping across a significant portion of the
solar surface and leaving behind a dimming area (Thompson & Myers, 2009; Mei
et al., 2012b). In their propagation, refraction and reflection may occur to
them as they go through regions of different densities or through the
boundaries separating magnetic fields of different topologies, like that
between active regions and coronal holes (Gopalswamy et al., 2009; Shen et
al., 2019; Wang et al., 2009, 2015; Xie et al., 2019). Furthermore, if the
wave is energetic enough, the nearby magnetic fields of opposite polarity
could even be pushed to approach each other, leading to the occurrence of
magnetic reconnection (Zhou et al., 2020).

The quasi-periodic wave train (QFP, Liu et al. 2011) is an intriguing
phenomenon of the disturbance to the corona that was observed for the first
time by the Atmospheric Imaging Assembly (AIA) onboard the Solar Dynamics
Observatory (SDO, Lemen et al. 2012). It has been interpreted as the fast-mode
magnetoacoustic wave (Cooper et al., 2003), characterized by continuous,
narrow arc-like structures that propagate rapidly at speed up to 1416 km
s-1(Shen et al., 2013; Kumar et al., 2017). A QFP is most apparent in AIA 171
Å (the associated emission produced by Fe ix ions and peaking at temperature
of around 0.8 MK), and is often closely related to solar flares. Typically,
QFPs first appear more than 100 Mm away from the flare center, possibly due to
the time it takes for the disturbance to grow and become detectable after
propagating a certain distance (Roberts et al., 1983). Some QFP waves exhibit
periods consistent with that of the oscillating component of the flare light
curve, which might imply that the flare invokes the QFP waves (Miao et al.,
2019). The mechanisms behind the formation of quasi-periodic pulsations in
flares include the dispersion of perturbations, magnetic configuration
resonance, nonlinear processes during magnetic reconnection, and oscillations
leaking from the other atmospheric layers into the solar corona (Nakariakov &
Melnikov, 2009).

The magnetic configuration resonance is typically associated with long-period
global oscillations of coronal loops and does not excite QFPs, while the other
three mechanisms can excite QFPs. However, some observations suggest that QFPs
with multiple narrow arc-like wavefronts are also related to the coronal mass
ejection (CME) and appear ahead of CME, indicating that CME is also a possible
trigger of the QFP wave (Zhou et al., 2021). Liu et al. (2012) discovered QFPs
on the inner side of a global-scale EUV wave, and noticed that QFPs propagated
on outside of CMEs, propagating a distance of over 0.5 solar radius on the
solar surface, at velocities decreasing from 1400 to 650 km s-1. They found
that QFPs originated from the same location, but more significant in the lower
corona at heights less than 110 Mm, and their brightness rapidly attenuated as
the altitude increases. Due to the limitation to present observational
techniques, the origin of QFPs remains unclear. Several mechanisms have been
proposed to explain the relationship between QFPs and CMEs by Wang et al.
(2021a) and Shen et al. (2019) or flares by Roberts et al. (1984), Nakariakov
et al. (2004), Nisticò et al. (2014), Ofman et al. (2011), Ofman & Liu (2018),
but still need to be further investigated.

Many authors have investigated QFP waves in the corona via numerical
simulations. Yang et al. (2015) found that the QFP wave could result from
collisions of the plasmoids produced in the current sheet by the tearing mode
instability with the magnetic field nearby, and propagated at speed up to 1000
km s-1. Takasao & Shibata (2016) revealed that QFPs are excited by the local
oscillation on the top of flare loops, which is governed by magnetic
reconnection that takes place in the periodical fashion. Ye et al. (2021)
discovered QFPs at the trailing edge of CME and suggested a potential
connection to the turbulence.

On the other hand, for MHD wave, magnetic field structures filled with plasma
may act as effective waveguides (Roberts et al., 1983). Based on the
difference in the speeds of the waves inside and outside the waveguide, the
structure of the waveguide as well as the internal wave could be classified
into two categories. The first one corresponds to the case in which the speed
of the internal wave is lower than that of the fast magnetoacoustic wave
outside the waveguide, so the internal wave is well confined in the waveguide.
The second one is of the case that the internal wave propagates faster than
the fast magnetoacoustic wave outside the waveguide, so the internal wave
could leak to outside of the waveguide, and propagates outward continuously.
Ofman & Liu (2018) performed simulations of QFPs using a three-dimensional MHD
model and successfully obtained various parameters for QFPs, including
waveform, amplitude, wavelength, and propagating speed, as well as duplicated
the reflection of QFP on the closed magnetic field structures. They also found
that the speed and energy flux density of QFPs rapidly decrease with the
altitude, which could be ascribed to the divergent characteristics of the
funnel-shaped magnetic structure and rapid weakening of the magnetic field at
higher altitude. Furthermore, Shen et al. (2022) noted an interesting feature
of QFPs, namely, wave amplitude grows initially, and decreases later on during
propagation. This phenomenon may result from the combination of the growth in
the amplitude of the wave and the expanding feature of the coronal magnetic
field configuration that serves as waveguides. Here the growth of the wave
amplitude is believed to result from the gravity stratification of the
atmosphere .

Regarding the leakage of the wave from the waveguide-like magnetic structure
in the corona, the seminal analytical work on leaky waves was conducted by
Cally (1986), and a comparative study of the leaky and trapped regimes was
performed by Nakariakov et al. (2012). Pascoe et al. (2013, 2014) studied the
leakage of the wave from the funnel-shaped magnetic structure. They saw the
leaking wave propagating outward, exhibiting multiple “wing”-shaped
wavefronts, which might account for the QFP wave observed in the eruption. Shi
et al. (2023) performed three-dimensional MHD simulations to study the impact
of the kink-like velocity perturbation to the flare loop (modeled by the
straight cylinder). Their results revealed that the straight cylinder could
behave as waveguides, allowing long-period (57 s period) and short-period (5.8
s period) oscillations to occur within the waveguide. They used a fast
gyrosynchrotron code to produce synthetic signals in the microwave band,
explaining the generation of quasi-periodic pulsations (QPPs) in solar flares.
It is worth noting that although not explicitly mentioned in their study,
their results suggest that when short-period oscillation signals pass through
the boundary of the cylinder, a portion of the oscillation leaks outside the
cylinder and forms multiple wavefronts, which is similar to the multiple QFP
waves observed in the solar eruption. Analyzing the observational data from
the 1.6-meter Goode Solar Telescope, Yuan et al. (2023) discovered persistent
transverse waves in the chromospheric umbral fibrils of a sunspot with strong
magnetic field. They reproduced these waves using a two-fluid MHD simulation.
Interestingly, their results indicate that perturbations confined within the
chromospheric fibrils (acting as waveguide structures) oscillated back and
forth, and leaked out at the boundaries, giving rise to multiple wavefront
signals.

In order to investigate the QFP phenomena and mechanisms for their generation
during solar eruptions, we will first examine the evolution in the coronal
magnetic structure including a magnetic flux ropes through a set of quasi-
static states to eventually losing the equilibrium. Then, we follow the
consequent dynamic evolution in the system, investigate the disturbance to the
background environment and the occurrence of QFP waves in this process, look
into the physical property and mechanism triggering the QFP wave. We will
provide a detailed introduction of the model and the relevant numerical
techniques used in the numerical simulations in Section 2, and display our
results in Section 3. Synthetic images deduced according to the numerical
experiments and the associated comparisons with observations will be given in
Section 4. Finally, we summarize this work in Section 5.

##  2 Descriptions of the numerical method

We utilize the MPI-AMRVAC code (Keppens et al., 2012; Porth et al., 2014; Xia
et al., 2018) to study the evolution in the coronal magnetic configuration
through a set of quasi-static equilibrium states to the loss of equilibrium.
Our calculations are involved in solving the following MHD equations, which
include the gravity and the thermal conduction:

$\partial_{t}\rho+\bigtriangledown\cdot(\rho\textbf{v})=0,$ (1)

$\partial_{t}(e)+\nabla\cdot[(e+p+\frac{1}{2\mu_{0}}|\textbf{B}|^{2})-\frac{1}{\mu_{0}}(\textbf{v}\cdot\textbf{B})\textbf{B}]=\rho\textbf{g}\cdot\textbf{v}+\nabla\cdot[\kappa_{\parallel}(\nabla
T\cdot\hat{\textbf{B}})\hat{\textbf{B}}],$ (2)

$\partial_{t}(\rho\textbf{v})+\bigtriangledown\cdot[\rho\textbf{v}\textbf{v}+(p+\frac{1}{2}|\textbf{B}|^{2})\textbf{I}-\textbf{B}\textbf{B}]=\rho\textbf{g},$
(3)

$\partial_{t}\textbf{B}=\bigtriangledown\times(\textbf{v}\times\textbf{B}-\eta\bigtriangledown\times\textbf{B}),$
(4)

where $\rho$, v, B, $p$, $e$, and $T$ are density, velocity, magnetic field,
pressure, energy, and temperature, respectively. In the description of thermal
conduction in the plasma, we employ the classical Spitzer model (Spitzer,
1962) with $\kappa_{\parallel}$= 9000 W m${}^{-1}~{}$K-1 for the coefficient
of the thermal conductivity parallel to the magnetic field. The gravity is
$\textbf{g}(y)=g_{0}/(1+y/r_{0})^{2}\hat{\textbf{y}}$ with $g_{0}$
($g_{0}=274$ m s-2) being the value of the gravity on the solar surface,
$r_{0}$ ($r_{0}=6.9551\times 10^{8}$ m) the radius of the Sun, and
$\hat{\mathbb{y}}$ the unit vector in the $y$-direction. Due to numerical
diffusion in calculations, we set $\eta$ to 0 for simplicity. We note here
that the thermal conduction parallel to the magnetic field is included in our
calculations. The point of doing so is that the thermal conduction in the
coronal environment is a dominant process of the energy transportation, and
that the field-aligned conduction is powerful compared to that perpendicular
to the field in the corona. Therefore including a field-aligned thermal
conduction in calculations here is of consideration for the model to be as
realistic as possible.

In addition, the energy $e$ and the pressure $p$ are given as:

$e=\frac{1}{2}\rho v^{2}+\frac{p}{\gamma-1}+\frac{B^{2}}{2\mu_{0}},$ (5)
$p=2k_{b}\rho(y)T(y)/m_{H},$ (6)

where $\gamma=5/3$. To suppress spurious oscillations in our computations, we
employ a third-order asymmetric TVD limiter and utilize a “three-step” scheme
of the third-order predictor-corrector. Regarding the divergence-free of
magnetic field, we apply the “typedivbfix=linde” condition in practical
calculations to maintaining the numerical conservation of the magnetic field
divergence ($\nabla\cdot\textbf{B}=0$). According to Equation (5), in low-beta
plasma atmospheres, numerical errors in the magnetic energy may potentially
exceed the internal energy, leading to negative values for the internal energy
or gas pressure when solving the complete energy equation. To address this
issue, we simultaneously solve both the complete energy equation and the
internal energy equation.

###  2.1 the gravitationally stratified atmosphere

Calculations are conducted in the Cartesian coordinate system $(x,y)$, with
the origin located at the bottom of the photosphere. The $x$-axis lies on
bottom of the photosphere, and the $y$-axis is the radial direction away from
the Sun. The simulation domain spans over $-1.6L_{0}\leq x\leq 1.6L_{0}$ and
$0\leq y\leq 3.2L_{0}$, where $L_{0}=2\times 10^{10}$ cm. In our simulation,
we employ a two-layer atmospheric structure to model the solar atmosphere,
where the region with higher density ($y<h_{p}$) represents the photosphere,
and the region with lower density ($y>h_{p}$) represents the corona. The
temperature of the gravitationally stratified atmosphere is set to:

$T(y)=\left\\{\begin{aligned} &T_{p}=2\times 10^{3}~{}\rm K,\quad&0\leq y\leq
h_{p}&,\\\ &T_{c}=10^{6}~{}\rm K,\quad&h_{p}<y&.&\end{aligned}\right.$ (7)

For the photospheric layer, we follow the setting of Mei et al. (2012b), and
take $h_{p}=2.5\times 10^{7}$ cm. The isothermal property of the atmosphere in
the low corona yields the gas pressure:

$p(y)=\left\\{\begin{aligned}
&p_{p}\exp\left[-\frac{y}{\lambda_{p}(1+y/r_{0})}\right],\quad&0\leq y\leq
h_{p}&,\\\
&p_{c}\exp\left[-\frac{y-h_{p}}{\lambda_{c}(1+y/r_{0})(1+h_{p}/r_{0})}\right],\quad&h_{p}<y&,\end{aligned}\right.$
(8)

where $p_{c}$ and $p_{p}$ are the pressure at the bottom of the corona and the
bottom of the photosphere, respectively, with
$p_{p}=p_{c}\exp\left[h_{p}/\lambda_{p}/(1+h_{p}/r_{0})\right]$. Here,
$\lambda_{p}$ and $\lambda_{c}$ are the pressure scale heights for the
photosphere and the corona, respectively:

$\lambda_{p}=\frac{2+3H_{ea}}{1+4H_{ea}}\frac{k_{b}T_{p}}{m_{p}g_{0}},$ (9)
$\lambda_{c}=\frac{2+3H_{ea}}{1+4H_{ea}}\frac{k_{b}T_{c}}{m_{p}g_{0}},$ (10)

where $H_{ea}$, $k_{b}$, and $m_{p}$ represent the helium abundance, Boltzmann
constant, and proton mass, respectively. The density distribution could be
determined according to the isothermal property of the atmosphere.

###  2.2 The structure of magnetic field and filament

Following Isenberg et al. (1993), we model the filament with a current-
carrying magnetic flux rope levitating in the corona. To ensure a smooth
transition of the parameters from the flux rope to the background corona, we
impose a thin shell outside the flux rope with a radius of $R$. This shell is
known as the PCTR (prominence-corona transition region, e.g., see also Mei et
al. 2012a), with a thickness of $\Delta=R/5$. Based on the distance $r_{-}$
from any point in space to the center of the flux rope
($r_{-}=\sqrt{x^{2}+(y-h)^{2}}$, where $h$ is the initial height of the flux
rope), the entire simulation domain is divided into three regions: the
interior of the flux rope (Ar1), for $0\leq r_{-}<(R-\Delta/2)$; the PCTR
(Ar2), for $(R-\Delta/2)\leq r_{-}\leq(R+\Delta/2)$; and the exterior of the
flux rope (Ar3), for $r_{-}>(R-\Delta/2)$. According to Mei et al. (2012a),
the current density of the flux rope is:

$j_{z}(r_{-})=\left\\{\begin{aligned} &j_{0},\quad&r_{-}\in Ar1,\\\
&\frac{j_{0}}{2}\left\\{\cos\left[\frac{\pi(r_{-}-R+\Delta/2)}{\Delta}\right]+1\right\\}&r_{-}\in
Ar2,\\\ &0,\quad&r_{-}\in Ar3,\end{aligned}\right.$ (11)

where $j_{0}$ is a constant.

The initial magnetic field results from three sources: the internal current of
the flux rope centered at $y=h$, the mirror current centered at $y=-h$, and
the buried magnetic quadrupoles beneath the photosphere at $y=-d$. So the
initial magnetic field is described by:

$B_{x}=B_{\phi}(r_{-})(y-h)/r_{-}-B_{\phi}(r_{+})(y+h)/r_{+}-B_{\phi}(R+\Delta/2)Md(R+\Delta/2)(y+d)(3x^{2}-(y+d)^{2})/r_{d}^{6},$
(12)

$B_{y}=-B_{\phi}(r_{-})x/r_{-}+B_{\phi}(r_{+})x/r_{+}-B_{\phi}(R+\Delta/2)Md(R+\Delta/2)x(-x^{2}+3(y+d)^{2})/r_{d}^{6},$
(13)

where

$B_{\phi}(r)=\left\\{\begin{aligned} &-j_{0}Jr/2,&\quad r\in Ar1\\\
&-j_{0}J((R-\Delta/2)^{2}/2-(\Delta/\pi)^{2}+r^{2}/2+(\Delta
R/\pi)sin[\pi(r-R+\Delta/2)/\Delta]+\\\
&(\Delta/\pi)^{2}cos[\pi(r-R+\Delta/2)/\Delta])/(2r),&r\in Ar2\\\
&-j_{0}J[R^{2}+(\Delta/2)^{2}-2(\Delta/\pi)^{2}]/(2r),&\quad r\in
Ar3\end{aligned}\right.$ (14)

$r_{+}=\sqrt{x^{2}+(y+h)^{2}},$ (15) $r_{d}=\sqrt{x^{2}+(y+d)^{2}},$ (16)
$M=125\sigma/32,$ (17)

with $r_{+}$ and $r_{d}$ being the distances from any point in space $(x,y)$
to the mirror current and the quadrupole, respectively, and $\sigma$ denoting
the relative strength of the background magnetic field (interested readers
refer to Forbes 1990).

In the magnetohydrostatic equilibrium, the initial gas pressure within the
flux rope is given by:

$p_{f}(y)=p(y)-\int_{r_{-}}^{\infty}j_{z}(r^{\prime})B_{\phi}(r^{\prime})dr^{\prime}+p_{Bz}.$
(18)

The second term on the right-hand side of equation (18) is the pressure
generated by the azimuthal magnetic field, and the third term on the right-
hand side is the pressure generated by the axial magnetic field. The external
equilibrium condition of the flux rope determines the intensity of the
electric current within the flux rope, while the pressure due to the azimuthal
magnetic field inside the flux rope governs its internal equilibrium. Prior to
numerical experiments, we adjust the axial magnetic pressure according to
specific requirements to regulate the internal pressure of the flux rope so
that its temperature and density could be as close as possible to those in the
true corona.

To smoothly connect the temperature inside the flux rope to that in the
background, we set the temperature distribution as follows:

$T(r_{-})=\left\\{\begin{aligned} &T_{f},\quad&r_{-}\in Ar1\\\
&(T_{c}-T_{f})(r_{-}-R+\Delta/2)/\Delta+T_{f},\quad&r_{-}\in
Ar2\end{aligned}\right.$ (19)

where $T_{f}=5\times 10^{4}$ K and $T_{c}=10^{6}$ K, the density distribution
inside the flux rope can be obtained through the ideal gas law.

Since we need to study the gradual evolution in the background magnetic field
driving the coronal magnetic field to evolve through a set of stable
equilibrium to a critical state in a quasi-static fashion, to ensure
consistency across all sets of numerical experiments, the simulation grid for
each experiment is uniformly set to $3840\times 3840$. Figure 1(a) illustrates
the initial magnetic configuration and the density distribution in the
gravitationally stratified atmosphere, Figure 1(b) displays the distribution
of fast speed. We notice that the speed around the flux rope is almost
uniform, which allows the wave around the flux rope to propagate almost
isotropically. This issue will be further discussed later.

Figure 1: Initial distributions of (a) the plasma density and (b) fast speed.
The black solid curves are the magnetic field lines and point A in (a) is
selected for analyzing the period of the QFP as shown later in Fig 6.

##  3 Results

Figure 2: Variations of the equilibrium height of the flux rope versus the
strength of the background field. The solid line in panel (a) is the
equilibrium curve that was duplicated from the analytical solution of Isenberg
et al. (1993), and the black cross marks the critical point. All the solid
dots on the equilibrium curve are used as the initial location of the flux
rope and the corresponding strength of the background field for numerical
experiments. Blue triangles mark positions where the flux rope finds the
equilibrium in numerical experiments, while red triangles outline the boundary
above which no equilibrium location exists for the flux rope in numerical
experiments. The inset in (a) displays more details of equilibrium locations
obtained in both theory (continuous curve) and numerical experiments (blue
triangles). The flux rope at locations marked by crosses on the curve are
still in equilibrium in theory, but not in numerical experiments. Instead the
red asterisk and the red triangle of the same value of $\sigma$ mark the
critical point for the numerical experiment. Panel (b) shows the temporal
evolution of the flux rope location in the corresponding experiments with the
initial locations marked by solid dots, asterisk, and crosses on the
theoretical curve displayed in the inset of (a).

Isenberg et al. (1993) studied analytically the change in the background
magnetic field driving the coronal magnetic field structure to evolve through
a set of stable equilibrium states to a critical state in a quasi-static
approach. They described the equilibrium height $h$ of the flux rope as a
function of the background magnetic field $\sigma$ (see the black solid line
in Figure 2(a)). Before investigating the quasi-periodic fast-propagating
(QFP) wave in the early stage of the solar eruption, it is essential to verify
whether the process investigated by Isenberg et al. (1993) can indeed occur.
In other words, we need to confirm that the coronal magnetic configuration
would indeed lose the equilibrium after evolving quasi-statically through a
set of equilibrium states to the critical state in response to the change in
the background field. This is important for two reasons: first, an unstable
state should evolve from the stable state, so the magnetic configuration in
reality needs to be in equilibrium for a while prior to the eruption. Second,
numerical experiments starts with a configuration in non-equilibrium might
bring unphysical issues into calculations, and yielding the unrealistic
results and conclusions.

Figure 2(a) shows the continuous curve taken from Isenberg et al. (1993),
depicting the equilibrium height $h$ of the magnetic flux rope as a function
of the strength of the background field $\sigma$. Due to discrepancies between
analytical and numerical solutions, when we tentatively take the parameters
from the continuous curve (black dots) for the initial magnetic configuration
in our numerical experiments, it is not surprising that the equilibrium
locations (as indicated by the row of blue inverted triangles located at the
bottom) of the flux rope in the numerical experiment are not located on the
continuous curve, instead the equilibrium locations determined in the
numerical experiment always deviate from those deduced from the analytical
solution. Moreover, the closer the tentative positions are to the critical
point on the analytic curve (near the black cross), the larger the deviation
is. The four black dots on the curve correspond to $\sigma=0.88,0.86,0.84$,
and 0.82 respectively. The results indicate that when the magnetic flux rope
is initially located at heights corresponding to these four points on the
equilibrium curve, it undergoes oscillations for a while and eventually finds
a stable equilibrium position at the height marked by the triangles closest to
the black equilibrium curve after self-adjustment. This demonstrates that the
stable equilibrium states and the associated quasi-static evolution in the
configuration described by Isenberg et al. (1993) can indeed be realized.
Driven by the change in the photospheric magnetic field and plasma, the system
evolves quasi-statically from one stable equilibrium state to another.

However, a change occurs near the critical point on the equilibrium curve
(marked by the cross). The critical point is the joint where the stable
section on the equilibrium curve connects to the unstable one. In reality and
numerical experiment, no magnetic configuration could exist if its parameters
are located on the unstable section of the equilibrium curve. So the loss of
equilibrium in the magnetic configuration inevitably occurs as the system
evolves to the critical point along the stable section of the equilibrium. In
this study, we found that when the initial position is located at point
$\sigma=0.80$ (the location marked by symbol “+” in Figure 2(a)), the flux
rope could no longer find equilibrium location through self-adjustment during
the consequent evolution, and the motion of the flux rope becomes highly
dynamic.

Evolution behaviors of the system at $\sigma=0.82$ and $\sigma=0.80$ suggest
the existence of a critical value between these two $\sigma$ values, to which
the corresponding magnetic configuration reaches to the critical state where
stable and unstable equilibria join. To locate the critical point, we
conducted four sets of test experiments between $\sigma=0.82$ and
$\sigma=0.80$, using $\sigma$= 0.815, 0.8125, 0.810, and 0.8075. The change in
the height of the flux rope versus time for each group of experiments is shown
in Figure 2(b). From the figure, we can see that when $\sigma$ = 0.820, 0.815,
and 0.8125, the flux rope quickly finds a new equilibrium several rounds of
oscillation. When $\sigma$ = 0.8075 and 0.805, the flux rope fails to find an
equilibrium and rapidly moves upward. When $\sigma=0.810$, the flux rope
oscillates and remains at a certain height for a while but eventually loses
its equilibrium. This indicates that the catastrophe occurs at this point and
it is the critical point found in our numerical experiments, denoted by the
red pentagram on the analytical curve in Figure 2(a). The height at which the
flux rope could stay for a while in the numerical experiment is marked by the
red triangle in the inset of Figure 2(a). Thus, connecting those triangles
closest to the equilibrium curve in Figure 2(a) yields the equilibrium curve
for our numerical experiments. Figure 2(a) also shows that the equilibrium
curve obtained from numerical experiments in the present work differs slightly
from the analytical results given by Isenberg et al. (1993). This difference
is primarily due to the existence of the gas pressure, gravity, and numerical
dissipation in simulations (e.g., see also Chen et al. 2022). However, as
expected, these issues do not introduce significant deviations between the
numerical and the analytical results because the lower corona is very close to
a force-free environment.

Moreover, in order to find the regions near the equilibrium curve where the
flux rope can stably stay, we performed a set of test calculations, and
roughly outlined a range of values of the relevant parameters that govern the
equilibrium property of the system. We find that the system could eventually
reach the equilibrium if the initial value of the relevant parameters were
taken in the region below the red triangles and above the solid curve shown in
Figure 2(a). If the initial values were taken somewhere below the solid curve,
the flux rope will move upward and find the equilibrium position in that
region as well. The system will never find the equilibrium as the initial
value is taken at or above the red triangle, or in the region at the left to
the red asterisk on the curve.

So far, we have confirmed that the magnetic configuration shown in Figure 1
(a) can evolve, in response to the change in the background field, quasi-
statically through a set of stable equilibrium configurations to the critical
configuration, and then loses its equilibrium, triggering the eruption. The
main purpose of this work is to investigate the wave phenomenon in the
eruption, so to save computational resources, we directly start our numerical
experiments at $\sigma=0.81$ and focus on the issues related to the waves
occurring in the eruption.

###  3.1 The global evolution

Figure 3: Evolutionary features of quasi-periodic fast propagating (QFP) waves
at different times. The rows represent density, pressure, horizontal velocity,
and vertical velocity, respectively. The arrows indicate the features
described in section 3.1, where “w1” and “w2” refer to the first and second
wavefronts, respectively.

Figure 3 illustrates evolutionary features in the system at different times.
In the early stages of evolution at $t=16$ s, we see an arc-shaped wavefront
(labeled as “w1”) that includes the flux rope propagating outward while
expanding and extending in various directions simultaneously. As it extends to
the bottom boundary, an echo occurs as indicated by black arrows in Figure
3(g), a phenomenon commonly reported in both numerical experiments and
observation (e.g., see also Xie et al. 2019; Wang et al. 2009, 2015; Mei et
al. 2012b; Liu et al. 2012; Asai et al. 2012). In addition to the density
distribution, we also examine the distributions of gas pressure $p$,
horizontal velocity component $v_{x}$, and vertical velocity component $v_{y}$
(see rows 2 to 4 in Figure 3). We notice that, comparing with the distribution
of the density, those of the velocity and the pressure reveal more fine
structures. These features may not be easily recognized in images that
basically display the brightness of targets. If the information about the
Doppler shift could be obtained, say via spectroscopic approach, they may
become apparent.

At $t=33$ s, as shown in the second column of Figure 3, two semicircular wave
fronts (“w${1}$” and “w${2}$”) are propagating forward of the flux rope. In
contour plots for velocity (Figure 3(e) and 3(h)), the wave patterns are quite
rich and impressive. Previous numerical experiments (e.g., see also Xie et al.
2019; Wang et al. 2009, 2015; Mei et al. 2012b) usually displayed only one
wave front, while the eruption process investigated in the present work
produces three sets of wave fronts, which can even be clearly distinguished in
the density distribution plot (Figure 3(c)).

In general, the generated waves seem to propagate isotropically, i.e., the
local propagation speed does not seem to depend on the direction. This is
governed by the fact that the fast speed is almost uniform in the region near
the flux as shown in Figure 1(b). These wave fronts appear strikingly similar
to the quasi-periodic fast magnetoacoustic wave trains displayed by Liu et al.
(2012). We shall discuss these results later in more detail and compare them
with observation. Furthermore, density contours in Figure 3(a) through Figure
3(c) showed that the wavefront at the low altitude is more significant than
that at the high altitude . This difference may be attributed to the
gravitationally stratified atmosphere, in which the plasma is dense at low
altitudes and tenuous at high altitudes, therefore, the response of low
atmosphere to the wave is more apparent than that of high atmosphere.

###  3.2 Characteristics of Wavefronts

Figure 4: Physical quantities variations in slices crossing wavefronts. (a)
displays a snapshot of density and magnetic filed at t = 40 s. Magnetic field
lines are depicted by solid gray lines. (b-d) and (e-g) illustrate changes in
density, pressure and temperature along slices s1 and s2, respectively.

According to Figures 3(a)-3(c), these fronts appear sequentially in roughly
symmetric shapes ahead of the flux rope. To look into their physical
properties, we further analyze their apparent behaviors. Figure 4(a) displays
the density distribution at $t=40$ s, with black solid lines for magnetic
field lines. From the figure, we see that these two wavefronts propagate
successively on the solar surface to a distance of a few $10^{4}$ km, and an
obvious deflection of the magnetic field occurs on the wavefront with an
apparent enhancement in the tangential component of the magnetic field. This
suggests that these waves are the fast-mode shock.

To quantitatively reveal the property of these fronts, we analyze the changes
in the relevant parameters before and after sweeping of the front along two
white lines, s${1}$ and s${2}$. Here line s${1}$ crosses the front
horizontally, with the upstream (undisturbed medium) on the right side of the
front and the downstream (disturbed medium) on the left side, s${2}$ crosses
the shock front vertically, with the upstream on the top side and the
downstream on the bottom side of the front. Figures 4(b)-4d show the
distribution of density, pressure, and temperature along s${1}$, while Figures
4(e)-4(g) show the distributions along s${2}$.

We found that the density and the pressure along s${1}$ exhibit clear
discontinuities, while the temperature change is relatively smooth. After the
passage of the first wavefront, the plasma density, pressure, and temperature
increase by 1.25, 1.23, and 1.1 times, respectively; after the passage of the
second wave, the increases become 1.17, 1.18, and 1.08 times, respectively.
This indicates that the first wavefront is stronger and impacts the
environment more significantly than the second one. The distributions of the
above parameters along s${2}$ are similar to those along s${1}$, and also
exhibit clear discontinuity layers in density and pressure. After the passage
of the first wavefront, the plasma density, pressure, and temperature increase
by 1.48, 2.25, and 1.7 times, respectively; and on the second wavefront, the
corresponding changes become 1.48, 2.14, and 1.41 times, respectively. In
addition, we also noticed that the properties of the plasma and magnetic field
along s${1}$ and s${2}$ show slight differences because of the effect of
gravitational stratification. Furthermore, apparent similarities could be
recognized when comparing the magnetoacoustic shock structure shown in Figure
4 with the formation of the shock as a result of the leakage of fast-mode
waves from a plasma slab studied by Pascoe et al. (2017). In summary, the
observed discontinuity layers along s${1}$ and s${2}$ suggest again that these
disturbances are fast-mode shock.

Figure 5: Propagation speeds for the first wavefront “w1” in (a) and the
second wavefront “w2” in (b) versus time. The blue symbols “x” are for the
velocities of the wavefronts, while the red dots for the local fast speed that
could be deduced from Figure 1(b) correspondingly.

To further confirm the fast-mode properties of these perturbations, we compare
their velocities with the local fast magnetosonic wave speeds. Figures 5(a)
and 5(b) show the variations of the velocity versus time for the first and
second wavefronts, respectively. The blue symbols “x” are for the velocities
of the wavefronts, while the red dots for the local fast speed that could be
deduced from Figure 1(b) accordingly. In general, as the propagation distance
increases, the velocities of the wavefronts decrease. From Figure 5(a), we can
see that in the early stages, the velocity of the first wavefront is greater
than the local fast magnetosonic wave speed, with an average Mach number of
1.12 (relative to the fast magnetosonic wave speed); on the other hand, as the
wave further propagates, it gradually degrades into the fast magnetosonic
wave. As for the second wavefront, its initial velocity (approximately 1700 km
s-1) is lower than that of the first wavefront (approximately 1900 km s-1),
and its velocity is approximately equal to the fast magnetosonic wave speed.

Figure 6: Periodicity analysis of quasi-periodic wave trains. (a) denotes the
periodicity analysis of the fixed point outside the flux rope (marked as “A”
in Fig 1), and (b) denotes that of the central point of the flux rope.

Moreover, the wave patterns shown in Figure 3 are suggestive of the periodical
generation of these waves. To investigate this process, we conducted a wavelet
analysis of the energy variation versus time at a specific location outside
the flux rope (labeled as “A” in Figure 1). The results are presented in
Figure 6(a), revealing a disturbance signal with a period of approximately 20
seconds. Simultaneously, we also performed a similar analysis of the energy
disturbance at the center of the flux rope, and the results are displayed in
Figure 6(b). Interestingly, we found that the internal disturbance also
exhibits clear periodical feature, and the period is roughly the same as that
of the external disturbance. This suggests the connection of the origins of
the two disturbances to each other, which we will further discuss below.

###  3.3 Origin of Wave Trains

Figure 7: The velocity divergence at different times. The red arrow “p” is the
inner perturbation, and the black arrow “b” the boundary of the flux rope.
“w${1}$” and “w${2}$” are the first and the second wavefronts, respectively.
The animated velocity divergence images display the propagation of internal
disturbances and the formation of QFP waves in the time interval from 0 s to
50.67 s with a duration of 2 s. Static image only presents information at a
few key moments, while the animation shows the evolution. An animated version
of this figure is available.

We calculate the velocity divergence at different times (Figure 7 and the
supplementary movie) to analyze the origin and property of disturbances.
According to Equation (1), the velocity divergence describes the compression
($\nabla\cdot\bf{v}<0$) and expansion ($\nabla\cdot\bf{v}>0$) of the plasma,
making it an effective indicator of the shock.

As mentioned earlier, in the numerical experiments, the initial position of
the flux rope is close to the critical point, so it stays in this position for
a short period before losing the equilibrium and then rapidly moving upward
after losing the equilibrium. Simultaneously, the internal equilibrium of the
flux rope is also disrupted. Because of the compressibility of the plasma,
once a disturbance occurs, it propagates outward in the form of the wave. With
losing the internal equilibrium, two waves are invoked inside the flux rope,
and propagate inward and outward, respectively. Some of the outward-
propagating waves eventually leave the flux rope and enter the surrounding
corona, while other part is reflected at the edge of the flux rope, forming
inward-propagating waves. The inward waves eventually converge at the center
of the flux rope and then rebound outward. As a result, we can see an image of
the flux rope constantly oscillating and waves continuously propagating
outward from its edges. The oscillations of the flux rope eventually dies down
as a result of the dissipation and the wave leakage into the external medium.

Figure 7 qualitatively describes such a process. At $t=1$ s (see Figure 7(a)),
an outward-propagating wave (labeled as “w${1}$”) is observed outside the flux
rope, corresponding to the first wave shown in Figure 3. It is caused by the
loss of equilibrium of the flux rope and its rapid upward motion. Meanwhile, a
circular inward-propagating disturbance (labeled as “p”) appears inside the
flux rope. Between w${1}$ and p, a circular boundary interface exists and does
not change significantly with time, which is actually the boundary of the flux
rope (labeled as “b”). From Figure 7(b) and 7(c), we realize that w1 continues
to propagate outward, while p contracts further toward the center of the flux
rope until it shrinks to a point.

Later, at $t=10$ s (see Figure 7(d)), the downward portion of w${1}$ has
propagated into the lower solar atmosphere, while p starts expanding outward
from the center of the flux rope. At $t=15$ s (see Figure 7(e)), the downward
portion of w${1}$ produces aforementioned echoes. Simultaneously, disturbance
p encounters boundary b, and a portion of it transmits through b and
propagates outward, forming the second wavefront (labeled as “w2” in Figure
7(e)), which is consistent with the results of Pascoe et al. (2017). Another
portion reflects and propagates back toward the center of the flux rope again.
The propagation of p as described above forms a cycle that continuously
triggers oscillations inside the flux rope, leaking waves on its surface until
the extra energy inside the flux rope totally dissipates. This scenario well
displays the creation and propagation of various waves/shocks around the flux
rope during the eruption, which could find observational counterparts easily.

##  4 Observational consequences

Figure 8: Comparison between observations and the simulation. The top row
displays the observations while the middle row the simulation. Specifically,
(a) depicts an event that occurred on May 7, 2012, where multiple wavefronts
were detected on both sides of a flux rope through the white-light coronagraph
COR1 (Kumar et al. 2017). (b) shows an event that occurred on September 8-9,
2010, where QFP waves were observed in the AIA image at 193 Å (Liu et al.
2012). (c) and (d) are the time-sequences of the brightness distribution along
curve s${3}$ and line s${4}$ in (b). (i) and (j) are magnified images of the
velocity divergence within the boxed regions of (g) and (h).

The quasi-periodic fast-mode magnetoacoustic waves propagating on the solar
surface triggered by CME have been reported by a few authors. Kumar et al.
(2017) reported several wavefronts around the flux rope observed by the white-
light coronagraph on board STERO-A (Figure 8(a)). Liu et al. (2012) performed
detailed studies of the similar phenomenon observed in EUV by SDO/AIA (Figure
8(b)).

Figures 8(c) and 8(d) show the time-sequences of the brightness distributions
along curve s${3}$ and line s${4}$ specified in Figure 8(b), revealing
information about the propagation of the wavefronts in various directions.
Figure 8(c) demonstrated multiple wavefronts (indicated by white arrows) are
successively ejected from the boundary of CME core. Box1 in the figure
highlights how the upper part of the wavefront leaks out. Figure 8(d)
describes the situation in the vertical direction along s${4}$, where multiple
downward wavefronts (indicated by white arrows) reach the solar surface, and
the echo could be recognized. Box2 in the figure outlines the region where the
waves interact with the lower solar atmosphere. Interested readers refer to
Xie et al. (2019) for more information regarding the identification and
analysis of the echo on the solar surface.

To compare with observations, we create synthetic images in white-light and
EUV wavelength on the base of the simulation results. A synthetic images is
basically constructed according to instrument response functions and atomic
databases. The instrument response function provides the response of the
detector in a given wavelength to the plasma radiation at a given temperature.
The detector responds to the radiation in the way of:

$R_{f}=\int n_{e}^{2}f_{i}(T_{e},n_{e})dl.\qquad(DN/s)$ (20)

Here, $dl$ is the length element along the line of sight (LOS), $DN$ the
number of electrons produced by the photons received by the instrument,
$n_{e}$ the electron density, and $f_{i}(T_{e},n_{e})$ the response function
of the instrument for density $n_{e}$ and temperature $T_{e}$. A synthetic
image then could be obtained provided the form of $f_{i}(T_{e},n_{e})$, and
distributions of $T_{e}$ as well as $n_{e}$ in LOS are known.

Figure 8(e) shows the synthetic image in white-light, and Figures 8(f) through
8(h) provide the information in 193 Å. Comparing observations (Figures 8(a)
and 8(b)) with the synthetic images (Figures 8(e) and 8(f)), we find that the
simulation results exhibit multiple wavefronts in both wavelengths, but
observations display little bit different scenario such that these waves are
almost invisible in front of CME, but apparent on the side of CME. This could
be due to the fact that the density on the CME side is higher than that in
front of CME, and the disturbance on the CME side may produce more apparent
observational consequences; in the numerical experiments, on the other hand,
manifestations of the synthetic image could be manipulated by adjusting the
range of color bar.

To further justify that the multiple wavefront scenario is the oscillation
inside the flux rope in origin, we deduce the time-sequences of the brightness
distributions in directions parallel and perpendicular to the solar surface as
indicated by Figure 8(f). The results are displayed in Figures 8(g) and 8(h).
Comparing Figure 8(c) with Figure 8(g), we see several wave fronts in both
observational and synthetic images, and these wavefront structures originate
from the core region of the eruption. Furthermore, our simulation, on which
the synthetic images are created, definitely confirms that these wavefronts
result from leaking of the oscillation inside the flux rope.

In addition, a set of X-like patterns inside the flux rope are observed. In
Figure 8(i), we plot velocity divergence in the region inside box1 in Figure
8(g). The oscillating features in the radial direction of the flux rope could
be recognized, and the leakage of the oscillation through the boundary of the
flux rope is also seen clearly. In this process, the waves associated the
oscillation propagate both inward and outward inside the flux rope.
Concentrating and bouncing of the the inward wave at the center of the flux
rope constitute the X-shape feature as shown by shown in Figure 8(g) and
Pascoe et al. (2013).

Comparing Figure 8(d) with Figure 8(h), we find multiple wavefronts
propagating downward in both observations and synthetic images. To study the
interaction of the downward wavefronts with the solar surface, we calculate
the corresponding velocity divergence in the region indicated by box2 in
Figure 8(h), and the result is plotted in Figure 8(j), which shows clearly the
pattern of the echo. The similar result was also reported by Xie et al.
(2019), with the only difference being that Xie et al. (2019) showed a single
echo, whereas several echos could be identified in Figure 8(j).

##  5 Summary And Discussions

On the basis of the analytical model of Isenberg et al. (1993), we
investigated the evolution in a magnetic configuration that includes an
electric current-carrying flux rope and possesses a relatively simple
background magnetic field. Our results confirmed that the coronal magnetic
structure described by the model can quasi-statically evolve through a set of
configurations in the stable equilibrium to an unstable configuration,
eventually leading to an eruption in a catastrophic fashion. Comparing with
the works by Lin et al. (2001), Kolotkov et al. (2016, 2018), as well as Chen
et al. (2022), we notice that the evolutionary behaviors of the magnetic
structure studied here is simple and straightforward such that the
displacement of the flux rope only moves in the $y$-direction during both
quasi-static and dynamic stages. Lin et al. (2001), Kolotkov et al. (2016,
2018), and Chen et al. (2022) found that the motion of the flux rope could be
in both $x$\- and $y$-directions(here the $y$-direction is the vertical
direction) as the background field is complex, and even the loss of
equilibrium in the flux rope may not be able to develop to a plausible
eruption eventually. This might be more realistic.

However, the purpose of this work is to look into the origin of the QFP wave
during the eruption no matter how the eruption is triggered. So a simpler
magnetic configuration was chosen and the straightforward evolutionary
behavior in the configuration was studied. After confirming the occurrence of
the loss of equilibrium of the flux rope, we subsequently investigated the
occurrence of QFP waves in the solar eruption. Our results showed that arc-
like wavefronts successively appear on the surface of the flux rope and
propagate outward at speed of a few $10^{3}$ km s-1 as the flux rope moves
out.

We analyzed the mechanism and propagating features of these wavefronts , and
found out they were moving at speeds faster or marginally faster than the
local fast magneto-acoustic wave speed, which implies their nature of the
fast-mode shock. Looking into the velocity divergence of the wavefront, we
discovered that these QFP waves originate from disturbances within the flux
rope as a result of the loss of the internal equilibrium, and part of the
disturbance leaks through the flux rope boundary out to the surrounding
corona, forming the QFP waves observed. When QFP waves reach the lower solar
atmosphere, they generate multiple echos. In summary, the role of flux rope in
the origin of the wave trains observed in the solar eruption is twofold.
First, the loss of the internal equilibrium of the flux rope ignites
oscillations inside the flux rope; and second, the flux rope behaves like a
leaking wave guide from which the oscillation successively leaks out to the
surrounding corona, invoking the wave train with period of several tens of
second. The numerical simulation of Wang et al. (2021a) showed that the period
of the wave train is around 30 s.

In addition to the above point regarding the origin of wave trains, Wang et
al. (2021a) suggested that the wave train is a component of the large-scale
EUV waves. But this mechanism is hard to account for the fact that various
wave fronts come from roughly the same region, and that velocities of
different wave fronts are approximately the same.

Numerical experiments performed here duplicated the large-scale QFP waves
observed in the lower solar corona by Liu et al. (2012) as well. In addition,
Kumar et al. (2017) also identified QFP waves in the high corona on the basis
of analyzing the data from the white-light coronagraph COR1 on STEREO-A. To
our knowledge, no author reports simultaneous identification of QFP waves in
both the low and the high solar corona regions in a single eruptive event.
Furthermore, because of the limit to the computational resources, our
numerical experiments could not investigate QFP in both regions
simultaneously, either; so it is very hard to reveal the relations of the wave
train in the low to that in the high corona for the time being. But the
possibility exists that the QFP wave in the low corona propagates upward,
reaches to the large altitudes, and evolves into that in the high corona.

Overall, to definitively determine the mechanisms for their origins, it is
necessary to identify and compare the processes of QFP wave generation and
propagation in both the low and higher solar corona regions. This requires
observational instruments with large FOV that cover both the low and the high
corona regions with sufficiently high time cadence. The existing observational
data may not fully meet these requirements, we hope that advanced technology
and observational instruments in the future would be able to help solve this
puzzle.

This work was supported by the National Key R&D Program of China
No.2022YFF0503804, the NSFC grants 11933009, 12273107, U2031141, and 12073073,
grants associated with the Yunling Scholar Project of the Yunnan Province, the
Yunnan Province Scientist Workshop of Solar Physics, the Applied Basic
Research of Yunnan Province (2019FB005, 202101AT070018), Yunnan Key Laboratory
of Solar Physics and Space Exploration of Code 202205AG070009. JY and ZM. also
acknowledge the support by grants associated with the Yunnan Revitalization
Talent Support Program, the Foundation of the Chinese Academy of Sciences
(Light of West China Program). The numerical computation in this paper was
carried out on the computing facilities of the Computational Solar Physics
Laboratory of Yunnan Observatories (CosPLYO).

## References

  * Asai et al. (2012) Asai, A., Ishii, T. T., Isobe, H., et al. 2012, ApJ, 745, L18, doi: 10.1088/2041-8205/745/2/L18
  * Banerjee et al. (2021) Banerjee, D., Krishna Prasad, S., Pant, V., et al. 2021, Space Sci. Rev., 217, 76, doi: 10.1007/s11214-021-00849-0
  * Cally (1986) Cally, P. S. 1986, Sol. Phys., 103, 277, doi: 10.1007/BF00147830
  * Chen et al. (2022) Chen, Y., Ye, J., Mei, Z., et al. 2022, ApJ, 933, 148, doi: 10.3847/1538-4357/ac73ef
  * Cooper et al. (2003) Cooper, F. C., Nakariakov, V. M., & Williams, D. R. 2003, A&A, 409, 325, doi: 10.1051/0004-6361:20031071
  * De Moortel et al. (2000) De Moortel, I., Ireland, J., & Walsh, R. W. 2000, A&A, 355, L23
  * Forbes (1990) Forbes, T. G. 1990, J. Geophys. Res., 95, 11919, doi: 10.1029/JA095iA08p11919
  * Gopalswamy et al. (2009) Gopalswamy, N., Yashiro, S., Temmer, M., et al. 2009, ApJ, 691, L123, doi: 10.1088/0004-637X/691/2/L123
  * Isenberg et al. (1993) Isenberg, P. A., Forbes, T. G., & Demoulin, P. 1993, ApJ, 417, 368, doi: 10.1086/173319
  * Keppens et al. (2012) Keppens, R., Meliani, Z., van Marle, A. J., et al. 2012, Journal of Computational Physics, 231, 718, doi: 10.1016/j.jcp.2011.01.020
  * Kolotkov et al. (2016) Kolotkov, D. Y., Nisticò, G., & Nakariakov, V. M. 2016, A&A, 590, A120, doi: 10.1051/0004-6361/201628501
  * Kolotkov et al. (2018) Kolotkov, D. Y., Nisticò, G., Rowlands, G., & Nakariakov, V. M. 2018, Journal of Atmospheric and Solar-Terrestrial Physics, 172, 40, doi: 10.1016/j.jastp.2018.03.005
  * Kumar et al. (2017) Kumar, P., Nakariakov, V. M., & Cho, K.-S. 2017, ApJ, 844, 149, doi: 10.3847/1538-4357/aa7d53
  * 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
  * Li et al. (2020) Li, B., Antolin, P., Guo, M. Z., et al. 2020, Space Sci. Rev., 216, 136, doi: 10.1007/s11214-020-00761-z
  * Lin et al. (2001) Lin, J., Forbes, T. G., & Isenberg, P. A. 2001, J. Geophys. Res., 106, 25053, doi: 10.1029/2001JA000046
  * Lin et al. (2010) Lin, J., Shen, C., & Wang, H. 2010, in 38th COSPAR Scientific Assembly, Vol. 38, 4
  * Liu et al. (2012) Liu, W., Ofman, L., Nitta, N. V., et al. 2012, ApJ, 753, 52, doi: 10.1088/0004-637X/753/1/52
  * Liu et al. (2011) Liu, W., Title, A. M., Zhao, J., et al. 2011, ApJ, 736, L13, doi: 10.1088/2041-8205/736/1/L13
  * McLaughlin et al. (2018) McLaughlin, J. A., Nakariakov, V. M., Dominique, M., Jelínek, P., & Takasao, S. 2018, Space Sci. Rev., 214, 45, doi: 10.1007/s11214-018-0478-5
  * Mei et al. (2012a) Mei, Z., Shen, C., Wu, N., et al. 2012a, MNRAS, 425, 2824, doi: 10.1111/j.1365-2966.2012.21625.x
  * Mei et al. (2012b) Mei, Z., Udo, Z., & Lin, J. 2012b, Science China Physics, Mechanics, and Astronomy, 55, 1316, doi: 10.1007/s11433-012-4752-3
  * Miao et al. (2019) Miao, Y. H., Liu, Y., Shen, Y. D., et al. 2019, ApJ, 871, L2, doi: 10.3847/2041-8213/aafaf9
  * Nakariakov et al. (2004) Nakariakov, V. M., Arber, T. D., Ault, C. E., et al. 2004, MNRAS, 349, 705, doi: 10.1111/j.1365-2966.2004.07537.x
  * Nakariakov et al. (2012) Nakariakov, V. M., Hornsey, C., & Melnikov, V. F. 2012, ApJ, 761, 134, doi: 10.1088/0004-637X/761/2/134
  * Nakariakov & Kolotkov (2020) Nakariakov, V. M., & Kolotkov, D. Y. 2020, ARA&A, 58, 441, doi: 10.1146/annurev-astro-032320-042940
  * Nakariakov & Melnikov (2009) Nakariakov, V. M., & Melnikov, V. F. 2009, Space Sci. Rev., 149, 119, doi: 10.1007/s11214-009-9536-3
  * Nakariakov et al. (2003) Nakariakov, V. M., Melnikov, V. F., & Reznikova, V. E. 2003, A&A, 412, L7, doi: 10.1051/0004-6361:20031660
  * Nakariakov et al. (2021) Nakariakov, V. M., Anfinogentov, S. A., Antolin, P., et al. 2021, Space Sci. Rev., 217, 73, doi: 10.1007/s11214-021-00847-2
  * Nisticò et al. (2014) Nisticò, G., Pascoe, D. J., & Nakariakov, V. M. 2014, A&A, 569, A12, doi: 10.1051/0004-6361/201423763
  * Ofman & Liu (2018) Ofman, L., & Liu, W. 2018, ApJ, 860, 54, doi: 10.3847/1538-4357/aac2e8
  * Ofman et al. (2011) Ofman, L., Liu, W., Title, A., & Aschwanden, M. 2011, ApJ, 740, L33, doi: 10.1088/2041-8205/740/2/L33
  * Ofman & Wang (2002) Ofman, L., & Wang, T. 2002, ApJ, 580, L85, doi: 10.1086/345548
  * Pascoe et al. (2017) Pascoe, D. J., Goddard, C. R., & Nakariakov, V. M. 2017, ApJ, 847, L21, doi: 10.3847/2041-8213/aa8db8
  * Pascoe et al. (2013) Pascoe, D. J., Nakariakov, V. M., & Kupriyanova, E. G. 2013, A&A, 560, A97, doi: 10.1051/0004-6361/201322678
  * Pascoe et al. (2014) —. 2014, A&A, 568, A20, doi: 10.1051/0004-6361/201423931
  * Porter et al. (1994) Porter, L. J., Klimchuk, J. A., & Sturrock, P. A. 1994, ApJ, 435, 502, doi: 10.1086/174831
  * Porth et al. (2014) Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4, doi: 10.1088/0067-0049/214/1/4
  * Roberts et al. (1983) Roberts, B., Edwin, P. M., & Benz, A. O. 1983, Nature, 305, 688, doi: 10.1038/305688a0
  * Roberts et al. (1984) —. 1984, ApJ, 279, 857, doi: 10.1086/161956
  * Shen et al. (2019) Shen, Y., Chen, P. F., Liu, Y. D., et al. 2019, ApJ, 873, 22, doi: 10.3847/1538-4357/ab01dd
  * Shen et al. (2022) Shen, Y., Zhou, X., Duan, Y., et al. 2022, Sol. Phys., 297, 20, doi: 10.1007/s11207-022-01953-2
  * Shen et al. (2013) Shen, Y. D., Liu, Y., Su, J. T., et al. 2013, Sol. Phys., 288, 585, doi: 10.1007/s11207-013-0395-4
  * Shi et al. (2023) Shi, M., Li, B., Chen, S.-X., Guo, M., & Yuan, S. 2023, ApJ, 943, L19, doi: 10.3847/2041-8213/acb3c6
  * Spitzer (1962) Spitzer, L. 1962, Physics of Fully Ionized Gases
  * Takasao & Shibata (2016) Takasao, S., & Shibata, K. 2016, ApJ, 823, 150, doi: 10.3847/0004-637X/823/2/150
  * Thompson & Myers (2009) Thompson, B. J., & Myers, D. C. 2009, ApJS, 183, 225, doi: 10.1088/0067-0049/183/2/225
  * Thompson et al. (1998) Thompson, B. J., Plunkett, S. P., Gurman, J. B., et al. 1998, Geophys. Res. Lett., 25, 2465, doi: 10.1029/98GL50429
  * Van Doorsselaere et al. (2020) Van Doorsselaere, T., Srivastava, A. K., Antolin, P., et al. 2020, Space Sci. Rev., 216, 140, doi: 10.1007/s11214-020-00770-y
  * Wang et al. (2021a) Wang, C., Chen, F., & Ding, M. 2021a, ApJ, 911, L8, doi: 10.3847/2041-8213/abefe6
  * Wang et al. (2015) Wang, H., Liu, S., Gong, J., Wu, N., & Lin, J. 2015, ApJ, 805, 114, doi: 10.1088/0004-637X/805/2/114
  * Wang et al. (2009) Wang, H., Shen, C., & Lin, J. 2009, ApJ, 700, 1716, doi: 10.1088/0004-637X/700/2/1716
  * Wang et al. (2021b) Wang, T., Ofman, L., Yuan, D., et al. 2021b, Space Sci. Rev., 217, 34, doi: 10.1007/s11214-021-00811-0
  * Williams et al. (2002) Williams, D. R., Mathioudakis, M., Gallagher, P. T., et al. 2002, MNRAS, 336, 747, doi: 10.1046/j.1365-8711.2002.05764.x
  * Williams et al. (2001) Williams, D. R., Phillips, K. J. H., Rudawy, P., et al. 2001, MNRAS, 326, 428, doi: 10.1046/j.1365-8711.2001.04491.x
  * Xia et al. (2018) Xia, C., Teunissen, J., El Mellah, I., Chané, E., & Keppens, R. 2018, ApJS, 234, 30, doi: 10.3847/1538-4365/aaa6c8
  * Xie et al. (2019) Xie, X., Mei, Z., Huang, M., et al. 2019, MNRAS, 490, 2918, doi: 10.1093/mnras/stz2576
  * Yang et al. (2015) Yang, L., Zhang, L., He, J., et al. 2015, ApJ, 800, 111, doi: 10.1088/0004-637X/800/2/111
  * Ye et al. (2021) Ye, J., Cai, Q., Shen, C., et al. 2021, ApJ, 909, 45, doi: 10.3847/1538-4357/abdeb5
  * Yuan et al. (2023) Yuan, D., Fu, L., Cao, W., et al. 2023, Nature Astronomy, 7, 856, doi: 10.1038/s41550-023-01973-3
  * Zhou et al. (2020) Zhou, G., Gao, G., Wang, J., et al. 2020, ApJ, 905, 150, doi: 10.3847/1538-4357/abc5b2
  * Zhou et al. (2021) Zhou, X., Shen, Y., Su, J., et al. 2021, arXiv e-prints, arXiv:2109.02847.  https://arxiv.org/abs/2109.02847