# Ensemble Distillation Approaches for Grammatical Error Correction

###### Abstract

Ensemble approaches are commonly used techniques to improving a system by
combining multiple model predictions. Additionally these schemes allow the
uncertainty, as well as the source of the uncertainty, to be derived for the
prediction. Unfortunately these benefits come at a computational and memory
cost. To address this problem ensemble distillation (EnD) and more recently
ensemble distribution distillation (EnDD) have been proposed that compress the
ensemble into a single model, representing either the ensemble average
prediction or prediction distribution respectively. This paper examines the
application of both these distillation approaches to a sequence prediction
task, grammatical error correction (GEC). This is an important application
area for language learning tasks as it can yield highly useful feedback to the
learner. It is, however, more challenging than the standard tasks investigated
for distillation as the prediction of any grammatical correction to a word
will be highly dependent on both the input sequence and the generated output
history for the word. The performance of both EnD and EnDD are evaluated on
both publicly available GEC tasks as well as a spoken language task.

Index Terms—  Ensemble, Distribution Distillation, Dirichlet, Transformer,
Grammatical Error Correction

##  1 Introduction

Knowledge distillation (KD) is a general approach in machine learning where
predictions from a more complex model are modelled by a simpler model, or an
ensemble of models represented by a single model [2]. The aim is to reduce the
computational cost and memory requirements for deployment. It has been applied
successfully in many domains in deep learning, such as object detection [3],
natural language processing [4, 5], acoustic models [6, 7], and also in
adversarial defence [8]. This paper focuses on applying ensemble distillation
approaches for sequence to sequence tasks, in this case grammatical error
correction (GEC).

Deep ensembles generally perform better than single models, with the added
benefit of providing measures of prediction uncertainty [9, 10]. Unfortunately
these benefits come at the cost of higher computational power and memory. To
address this problem, ensembles can be distilled into a single model.
Normally, distillation results in information about the uncertainty in the
distilled model’s prediction being lost [11]. One approach to maintaining
these uncertainty measures after ensemble distillation is Ensemble
Distribution Distillation [12]. Here the distilled student models the
distribution of the categorical predictions from the ensemble, allowing it to
retain more information about its teacher ensemble. Usually a Dirichlet
distribution is used for the distribution over categorical distributions. This
model can be challenging to train as the student model is required to
represent significantly more information from the teacher ensemble than the
standard, cross-entropy based, ensemble distillation [12]. If the distilled
model is well trained, however, it should efficiently yield good results and
enable, for example, detecting whether inputs are in or out-of-distribution,
and yield insight into why a prediction is (un)reliable [12, 11, 13].

This paper examines distillation approaches for sequence to sequence tasks
such as machine translation and grammatical error correction. For these tasks
the inference computational cost is significantly larger than ”static” tasks,
as the decoding process must combine the sequence predictions from each
individual in the ensemble [6, 14]. Distillation schemes already exist which
tackle the problem of reducing an ensemble, to save resources while also
maintaining performance. Distribution distillation methods however, have not
been explored for sequence models [15, 16]. This distribution distillation
would enable efficient uncertainty measures for tasks such as speech-to-text
systems and neural machine translation. The application area explored here is
GEC [17]. The aim is to extend previous work on uncertainty for spoken
language assessment [13, 17] to the problem of providing feedback to the
learner, and one challenge for general GEC systems is the difference between
acceptable grammar for written and spoken English.

Initially, standard distillation approach for sequence-to-sequence models, in
particular auto-regressive sequence-to-sequence, are described. Note rather
than modelling the (approximate) sequence posterior as in [6, 14], this work
models the token-level (word) posterior as this enables a simple extension to
distribution distillation for this form of model, and can be efficiently used
to derive uncertainty measures [15]. The challenges of applying distribution
distillation are then discussed, building on top of the work done in [12]. Two
possible optimisation approaches for ensemble distribution distillation are
described, as well combining standard distillation and distribution
distillation models together. These approaches are all evaluated on the task
of performing grammatical error correction, on both written data, where the
training and test data are from the same domain, and speech data where there
is a mismatch. Uncertainty performance is assessed on these tasks, and
compared to rejection based on the manual, reference, corrections.

##  2 Ensemble Distribution Distillation

Knowledge and distribution distillation are motivated by the need for high
performance systems with low memory, power, and time, requirements, whilst
also yielding uncertainty measures [11, 12, 14, 15, 18]. There are many ways
in which this can be achieved: here we will present work on knowledge
distillation and distribution distillation. In this section, we will focus on
static models, which are used for image classification [11] and spoken
language assessment [13]. A teacher ensemble is represented by a set of
parameters $\\{\bm{\theta}^{(m)}\\}_{m=1}^{M}$, and the parameters of the
standard and distribution distilled students will be represented by
$\bm{\phi}$ and $\bm{\lambda}$, respectively (to signify a fundamental
difference). Assuming that the teacher ensemble are drawn from the posterior
${\tt P}(\bm{\theta}|\mathcal{D})$, then:

$\displaystyle{\tt P}(y|\bm{x};\bm{\theta}^{(m)})={\tt
P}(y|\bm{\pi}^{(m)}),\medspace\bm{\pi}^{(m)}=\bm{f}(\bm{x};\bm{\theta}^{(m)})$

where $\pi_{i}$ represents the probability ${\tt P}(y=\omega_{i})$ and
$\bm{x}$ is the input. The set of categoricals
$\\{\bm{\pi}^{(m)}\\}_{m=1}^{M}$ can then be used to train new students. In
standard distillation [2], a student $\bm{\phi}$ is trained to emulate the
predictive distribution by minimising:

$\displaystyle\mathcal{L}(\bm{\phi})$ $\displaystyle={\tt
KL}\bigg{(}\left(\frac{1}{M}\sum_{m=1}^{M}\bm{\pi}^{(m)}\right)\;\bigg{\|}\;{\tt
P}(y|\bm{x};\bm{\phi})\bigg{)}$

This method essentially aims to let the student predict the mean of the
teacher while information about the spread (or disagreement) between members
of the ensemble is lost. This means that although distillation can
significantly improve memory and power consumption, the uncertainty metrics
will be lost.

Instead, following [12], we let a new student model $\bm{\lambda}$ predict the
parameters $\bm{\alpha}$ of a Dirichlet distribution ${\tt
Dir}(\cdot\medspace;\bm{\alpha})$:

$\displaystyle{\tt P}(\bm{\pi}|\bm{x},\bm{\lambda})={\tt
Dir}(\bm{\pi};\bm{\alpha}),\medspace\bm{\alpha}=\bm{f}(\bm{x};\bm{\lambda}),\medspace\alpha_{0}=\bm{1}^{T}\bm{\alpha}$

The Dirichlet is a distribution over categoricals, and can encapsulate
information about the uncertainty and spread in
$\\{\bm{\pi}^{(m)}\\}_{m=1}^{M}$. The training can then be based on negative
log-likelihood (NLL):

$\displaystyle\mathcal{L}(\bm{\lambda})=-\frac{1}{M}\sum_{m=1}^{M}\ln\bigg{(}{\tt
Dir}(\bm{\pi}^{(m)};\bm{\alpha})\bigg{)}$

Now it is possible to quantify uncertainty in the prediction. Three measures
of uncertainty, and their associated source, are used in this work: total,
expected data and knowledge uncertainty. These are related by [19, 20]:

$\displaystyle\underbrace{\mathcal{I}[y,\bm{\pi}|\bm{x},\bm{\lambda}]}_{\begin{subarray}{c}Knowledge\\\
Uncertainty\end{subarray}}=\underbrace{\mathcal{H}\big{[}{\tt
P}(y|\bm{x};\bm{\lambda})\big{]}}_{\begin{subarray}{c}Total\\\
Uncertainty\end{subarray}}-\underbrace{\mathbb{E}_{{\tt
P}(\bm{\pi}|\bm{x};\bm{\lambda})}\big{[}\mathcal{H}[{\tt
P}(y|\bm{\pi})]\big{]}}_{\begin{subarray}{c}Expected\ Data\\\
Uncertainty\end{subarray}}$ (1)

These uncertainties have different properties depending on whether the input,
$\bm{x}$, is in the same or different domain (distribution) to the training
set. If $\bm{x}$ is in-domain (ID), the ensemble members should return
consistent predictions, giving low expected data uncertainty (DU). If $\bm{x}$
is out-of-domain (OOD), the members will generate inconsistent predictions
giving high total uncertainty (TU) [21, 22, 23].

##  3 Sequence Ensemble Distillation

Applying both standard and distribution distillation to sequence models adds
another layer of complexity. This section covers token-level distillation
schemes, as these allow uncertainties for individual tokens, words, as well as
being combined to yield sequence-level uncertainties [15]. Extending the
notation in the previous section to sequence-to-sequence models, the pair
$(\bm{x},\bm{y})$ denotes the input-output reference sequence pair, and when
necessary, $\bm{\hat{y}}$ represents the corresponding predicted sequence for
${\bm{x}}$. The teacher ensemble $\\{\bm{\theta}^{(m)}\\}_{m=1}^{M}$ now makes
predictions by taking expectations over the following ensemble member
predictions111For a discussion of alternative approaches to ensemble
predictions based on sequences see [15].:

$\displaystyle{\tt P}(y_{l}|\bm{y}_{<l},\bm{x};\bm{\theta}^{(m)})={\tt
P}(y_{l}|\bm{\pi}_{l}^{(m)}),\medspace\bm{\pi}_{l}^{(m)}=\bm{f}(\bm{y}_{<l},\bm{x};\bm{\theta}^{(m)})$

and students are instead trained on the set
$\\{\bm{\pi}^{(1:M)}_{l}\\}_{l=1}^{L}$.

###  3.1 Ensemble Distillation

For sequence models, distillation can be performed in multiple ways [14, 24].
The approach adopted in this work is token-level knowledge distillation, and
is one of the simplest methods. The teacher and student $\bm{\phi}$ use the
same reference back-history ${\bm{y}}_{<l}$ (teacher-forcing) and input
$\bm{x}$. The KL-divergence between the ensemble and student token-level
categorical distributions is then minimised:

$\displaystyle\mathcal{L}(\bm{\phi})=\frac{1}{L}\sum_{l=1}^{L}{\tt
KL}\bigg{(}\left(\frac{1}{M}\sum_{m=1}^{M}{\bm{\pi}}_{l}^{(m)}\right)\;\bigg{\|}\;{\tt
P}\big{(}y_{l}|\bm{y}_{<l},\bm{x};\bm{\phi}\big{)}\bigg{)}$ (2)

Extending distribution distillation described in the previous section, this
section introduces distribution distillation for sequence models. For token-
level distribution distillation, the student with parameters $\bm{\lambda}$
predicts the Dirichlet distribution with parameters $\bm{\alpha}_{l}$ for the
$l$-th token:

$\displaystyle{\tt P}(\bm{\pi}|\bm{y}_{<l},\bm{x};\bm{\lambda})={\tt
Dir}(\bm{\pi};\bm{\alpha}_{l}),\medspace\bm{\alpha}_{l}=\bm{f}(\bm{y}_{<l},\bm{x};\bm{\lambda})$

Given $\\{\bm{\pi}^{(1:M)}_{l}\\}_{l=1}^{L}$ from the sequence ensemble, the
distribution distilled model can be trained using negative log-likelihood
(NLL):

$\displaystyle\mathcal{L}(\bm{\lambda})=-\frac{1}{ML}\sum_{l=1}^{L}\sum_{m=1}^{M}\ln\bigg{(}{\tt
Dir}(\bm{\pi}_{l}^{(m)};\bm{\alpha}_{l})\bigg{)}$ (3)

Eq. (3) optimizes the parameters of the distribution distilled model directly
from the ensemble predictions for each time instance. Though this yields an
appropriate criterion to find ${\bm{\lambda}}$, predicting the distribution
over the ensemble members for all back-histories, it may be very challenging
to optimise the network parameters. To simplify this optimization, a two stage
approach may be adopted. First, the predictions of the ensemble for each back-
history is modelled by a Dirichlet distribution with parameters
${\tilde{\bm{\alpha}}}_{l}$ where

$\displaystyle\bm{\tilde{\alpha}}_{l}=\operatorname*{arg\,max}_{\bm{\alpha}}\bigg{\\{}\frac{1}{M}\sum_{m=1}^{M}\ln\Big{(}{\tt
Dir}(\bm{\pi}_{l}^{(m)};\bm{\alpha})\Big{)}\bigg{\\}}$

The distribution distillation parameters are now trained to minimize the KL-
divergence between this ensemble Dirichlet and the predicted distilled
Dirichlet:

$\displaystyle\mathcal{L}(\bm{\lambda})=\frac{1}{L}\sum_{l=1}^{L}{\tt
KL}\Big{(}{\tt Dir}(\bm{\pi};\bm{\tilde{\alpha}}_{l})\;\Big{\|}\;{\tt
Dir}(\bm{\pi};\bm{\alpha}_{l})\Big{)}$ (4)

This has the same general form as eq. (2) but is based on the KL-divergence
between distributions of distributions, rather than just the expected
posterior distribution based on the ensemble.

Once a distribution distilled sequence model has been obtained, the
probability of predicting class $c$ is then:

$\displaystyle{\tt
P}(y=\omega_{c}|\bm{y}_{<l},\bm{x};\bm{\lambda})=\mathbb{E}_{{\tt
P}(\bm{\pi}|\bm{y}_{<l},\bm{x};\bm{\lambda})}\Big{[}{\tt
P}(y=\omega_{c}|\bm{\pi})\Big{]}=\frac{\alpha_{l,c}}{\alpha_{l,0}}$

as would be expected when the student parametrises a Dirichlet. This shows
that achieving sequence based distribution distillation is a straightforward
generalization of both ensemble distribution and sequence distillation.

Both eq. (3) and (4) result in viable objectives for obtaining a distribution
distilled sequence model $\bm{\lambda}$, and from which uncertainty metrics
can be derived. Instead of calculating entropy over the classes as in eq. (1),
uncertainties for sequence models have to enumerate all possible sequences
$\bm{y}\in\mathcal{Y}$. Let $\bm{\Pi}=\\{\bm{\pi}_{l}\\}_{l=1}^{L}$ be the
sequence of categorical distributions from which the sequence $\bm{y}$ is
generated:

$\displaystyle{\tt P}(\bm{\Pi}|\bm{x};\bm{\lambda})=\prod_{l=1}^{L}{\tt
Dir}(\bm{\pi}_{l};\bm{\alpha}_{l}),\medspace{\tt
P}(\bm{y}|\bm{\Pi})=\prod_{l=1}^{L}{\tt P}(y_{l}|\bm{\pi}_{l})$

The sequence uncertainties for the student $\bm{\lambda}$ can be expressed as:

$\displaystyle\underbrace{\mathcal{I}[\bm{y},\bm{\Pi}|\bm{x};\bm{\lambda}]}_{\begin{subarray}{c}Knowledge\\\
Uncertainty\end{subarray}}=\underbrace{\mathcal{H}\big{[}{\tt
P}(\bm{y}|\bm{x};\bm{\lambda})\big{]}}_{\begin{subarray}{c}Total\\\
Uncertainty\end{subarray}}-\underbrace{\mathbb{E}_{{\tt
P}(\bm{\Pi}|\bm{x};\bm{\lambda})}\big{[}\mathcal{H}[{\tt
P}(\bm{y}|\bm{\Pi})]\big{]}}_{\begin{subarray}{c}Expected\ Data\\\
Uncertainty\end{subarray}}$ (5)

However, as noted in [15], calculating the uncertainties require the
intractable computation of enumerating all $\bm{y}\in\mathcal{Y}$. Instead the
same approximations made in [15] will be utilized here. Given $S$ sampled
sequences: $\forall\bm{y}_{<l}^{(s)}\subset\bm{y}^{(s)}\sim{\tt
P}(\bm{y}|\bm{x},\bm{\lambda})$, the following Monte-Carlo approximations can
be made:

$\displaystyle\mathcal{H}\big{[}{\tt
P}(\bm{y}|\bm{x};\bm{\lambda})\big{]}\approx$
$\displaystyle\frac{1}{S}\sum_{s,l}\mathcal{H}\big{[}{\tt
P}(y_{l}|\bm{y}_{<l}^{(s)},\bm{x};\bm{\lambda})\big{]}$ (6)
$\displaystyle\mathbb{E}_{{\tt
P}(\bm{\Pi}|\bm{x};\bm{\lambda})}\big{[}\mathcal{H}[{\tt
P}(\bm{y}|\bm{\Pi})]\big{]}\approx$
$\displaystyle\frac{1}{S}\sum_{s,l}\mathbb{E}_{{\tt
P}(\bm{\pi}|\bm{y}_{<l}^{(s)},\bm{x};\bm{\lambda})}\big{[}\mathcal{H}[{\tt
P}(y_{l}|\bm{\pi})]\big{]}$ (7)

The sequence-level uncertainties, or rates, in [15] may be obtained by
normalising the quantities in eqs. (6) and (7) by the sequence length.

###  3.2 Guided Uncertainty Approach

Training distribution distilled models can be significantly harder than
training to standard distillation. Hence, a two model approach will also be
explored. This is based on the observation that the distribution distilled
model tends to have high Spearman’s rank correlation with the ensemble
predicted uncertainties in teacher-forcing mode used in training. In contrast,
evaluating the model in free-run decoding, so predictions are based on
${\hat{\bm{y}}}_{<l}$, the same consistency between the model and the ensemble
was not observed. To address this, the distribution distilled model was fed
the back-history from the distilled model, as the distilled model was found to
be less sensitive to the teacher-forcing and free-running mismatch.

Assuming we have a distilled model $\bm{\phi}$ obtained from eq. (2), and a
distribution distilled model $\bm{\lambda}$ obtained from either eq. (3) or
(4), one can then perform free-run decoding according to:

$\displaystyle\bm{\hat{\pi}}_{l}=\bm{f}(\bm{\hat{y}}_{<l},\bm{x};\bm{\phi});\>\>\>\>\hat{y}_{l}\sim{\tt
P}(y_{l}|{\hat{\bm{y}}}_{<l},{\bm{x}};{\bm{\phi}})$
$\displaystyle\hat{\bm{\alpha}_{l}}=\bm{f}(\bm{\hat{y}}_{<l},\bm{x};\bm{\lambda})$

The distilled model $\bm{\phi}$ is used to predict the output sequence, which
then guides the second model $\bm{\lambda}$ to return higher quality
uncertainties that can be derived from
$\\{\hat{\bm{\alpha}_{l}}\\}_{l=1}^{L}$. Although this method (referred to as
guided uncertainty approach; GUA) does not yield the same efficiency as
standard distribution distillation, it ensures that the best attributes of
$\bm{\phi}$ and $\bm{\lambda}$ are maintained in testing.

##  4 Experimental Evaluation

###  4.1 Data and Experimental Setup

The data used in this paper is the same as those described in [25]. This data
has been manually annotated with grammatical errors and corrections. The
training set and FCE (a specific test set) have been taken from the Cambridge
Learner Corpus [26] and includes written exams of candidates with different L1
languages. The FCE test set is a held-out subset of the corpus and therefore,
in-domain (ID) with the training data. NICT-JLE [27] is a public speech corpus
based on non-native speech. Only manual transcriptions of these interviews
involving Japanese learners are available, along with the grammatical
corrections. No audio is available. BULATS [28, 29] is a corpus based on a
free speaking business English test, where candidates were prompted for
responses up to 1 minute. Both manual (BLT) and ASR (BLT${}_{\text{asr}}$)
transcriptions were used, the average ASR WER for this data was 19.5% (see
[30] for details). The candidates are drawn from 6 L1s and span the full range
of proficiencies in English. These sets are out-of-domain (OOD) and will be
used to test the performance of any uncertainties derived from distribution
distilled models and ensembles, see Table (1).

Set | Train | FCE | NICT | BLT |  BLT${}_{\text{asr}}$
---|---|---|---|---|---
# of Sentences | 1.8M | 2.7K | 21.1K | 3.7K | 3.6K
Ave. Length | 13.4 | 14.0 | 6.6 | 16.6 | 16.7
Domain | Ref. | ID | OOD | OOD | OOD
Table 1: Number of sentences, mean sentence length, and domain.

All models used in this work are transformers based on [31], with the default
parameters of 6 layers, $d_{model}=512$, $d_{ff}=2048$, 8 heads, and 10%
dropout. Input words were mapped to randomly initialized embeddings of
dimension 512. Random seeds were used to generate an ensemble of 5 models
$M=5$ [10]. As 5 models are used in the ensemble, this increases the memory
requirements by a factor 5 over a single model, and approximately 5 for the
decoding cost depending on the form of beam-search being used. It is possible
to utilize more advanced ensemble generation methods, but is not key for this
work [19, 20, 32]. Additionally, sequence models can make predictions in
different ways: expectation of products, and products of expectations [15].
Since it has been found in [15] that products of expectations performs better,
it will be adopted in this work when evaluating ensembles. Beam-search with a
beam of one will be used, and uncertainty metrics will be evaluated on this
single output.

As noted in [12], when training, the target categorical is often concentrated
in one of the classes, while the Dirichlet predicted by the student often has
a distribution spread over many classes. This implies that the common support
between the teacher and student is poor. Optimizing KL-divergence when there
is a small non-zero common support between target and prediction can lead to
instability, and therefore, temperature annealing will be used to alleviate
this [2]. In this work, when training a distribution distilled system, the
targets from the ensemble are annealed with a temperature $T$, following a
schedule from $T=10.0$ to $3.0$. Further reducing the temperature down to
$T=1.0$ resulted in instability during training. To remain consistent, all
results concerning ensembles will also be based on the final temperature
$T=3.0$.

The GEC performance metric used in this work is GLEU [33], and uncertainty
metrics will be based at the token level. Furthermore, the distribution
distilled models and ensembles will be compared when a fraction of samples are
rejected (based on some metric) and passed on to an oracle. In these cases,
the rejection will be based on the highest sentence level uncertainty metrics,
and will be compared to simply rejecting sentences based on their length. An
additional metric based on knowing the true target will be ${\tt
L}\cdot(1-{\tt GLEU})$ (referred to as manual), where ${\tt L}$ is the true
sequence length; the inclusion of length is important as sentence level GLEU
does not take length into account even though it has an effect when used at a
dataset level [34]. System performance will be evaluated using relative area
under the curve ${\tt AUC_{RR}}$ as defined by:

$\displaystyle{\tt AUC_{RR}}=\frac{{\tt AUC}-{\tt AUC}_{\tt random}}{{\tt
AUC}_{\tt manual}-{\tt AUC}_{\tt random}}$ (8)

following the same reasoning as in [35] to simplify comparison between
metrics; ${\tt AUC}_{\tt random}$ refers to fully random rejection. Finally
GLEU performance will also be reported at 10% rejection.

###  4.2 Results

Table 2 shows the GLEU performance of a range of different models. As expected
the ensemble performs best, showing performance gains over the individual
ensemble members. Distilling the ensemble, using eq. (2), again yielded
performance gains over the individual ensemble members, though not to the same
level as the original ensemble but at reduced memory and computational cost.

Test set |  Ind±σ | Ens. | Dist. | NLL | KL | GUA
---|---|---|---|---|---|---
FCE |  69.5±0.11 | 70.6 | 69.9 | 68.0 | 68.8 | 69.9
NICT |  47.2±0.20 | 48.0 | 47.8 | 44.7 | 45.7 | 47.8
BLT |  49.8±0.15 | 50.9 | 50.7 | 48.2 | 48.9 | 50.7
BLT${}_{\text{asr}}$ |  31.3±0.09 | 31.6 | 31.5 | 30.9 | 31.3 | 31.5
Table 2: GLEU performance of individual, ensemble and standard and
distribution distilled models. Ind gives shows the ensemble member mean GLEU,
and standard deviation.

Secondly, the table also shows the performance of distribution distilled
models. The simplified training using KL-divergence, eq. (4), outperforms NLL,
eq. (3), though neither matches the average performance of the ensemble
members. This illustrates the challenges in training sequence ensemble
distribution distillation models. As the KL model performed better, it will be
used as the uncertainty model in GUA, though both the NLL and KL based models
yielded high Spearman rank correlation coefficient with the ensemble
uncertainty measures. The same trend can be seen for both the in-domain data
(FCE) and the out-of-domain data.

Table 3 shows the average word-level uncertainties predicted for two datasets,
FCE and the most mismatched BLT${}_{\text{asr}}$.

Test Set | Model | TU | DU | KU
---|---|---|---|---
FCE | Ensemble | 8.41 | 8.24 | 0.16
GUA | 8.37 | 8.20 | 0.18
BLT${}_{\text{asr}}$ | Ensemble | 8.85 | 8.67 | 0.18
GUA | 8.85 | 8.65 | 0.20
Table 3: Word-level uncertainties for an in-domain (FCE) and out-of-domain
(BLT${}_{\text{asr}}$) dataset. Total (TU), Data (DU) and Knowledge (KU)
Uncertainties.

The guided model approach behaves similarly to the ensemble. High knowledge
uncertainty (KU) is, according to theory, an indication of a sample being OOD
and in this case as expected KU is higher for BLT${}_{\text{asr}}$—an out of-
of-domain test set. Furthermore, it can be seen that all uncertainties for
BLT${}_{\text{asr}}$ are higher than FCE across both models.

Table 4 shows relative AUC (AUCRR) for FCE and BLT${}_{\text{asr}}$, as well
as mix of the two which assess whether the uncertainty measures can detect the
more challenging OOD speech data. It is interesting that as a baseline simply
using the length of the of the output sequence (Length) is a good baseline, as
longer sentences will tend to have more complex grammatical structure and
opportunity for mistakes. KU performs best in rejecting challenging inputs.
for both the ensemble and GUA, outperforming simple length selection for all
conditions.

Test Set Model Length TU DU KU FCE Ensemble 0.701 0.734 0.734 0.740 GUA 0.701
0.736 0.735 0.750 BLT${}_{\text{asr}}$ Ensemble 0.895 0.914 0.914 0.909 GUA
0.895 0.902 0.901 0.917 FCE + Ensemble 0.810 0.840 0.840 0.843
BLT${}_{\text{asr}}$ GUA 0.813 0.837 0.837 0.856

Table 4: AUCRR performance. The combined FCE + BLT${}_{\text{asr}}$ refers to
the case when the system is run on both ID and OOD data.

Table 5 shows the performance when a fixed percentage, 10%, is ”rejected” and
manually corrected. These results can be compared to the baseline numbers in
Table 2, show significant gains from rejecting 10% of the data. Similar trends
to that shown in AUCRR can be seen, KU generally yields the best performance,
outperforming the baseline length approach in all conditions. GUA yields
similar performance to the original ensemble, but at reduced memory and
computational costs.

Test Set Model Length TU DU KU Manual FCE Ens. 80.5 80.8 80.8 81.0 83.4 GUA
80.1 80.5 80.5 80.5 82.9 BLT${}_{\text{asr}}$ Ens. 52.4 54.4 54.2 53.2 56.5
GUA 51.9 52.6 52.5 53.4 56.3 FCE + Ens. 65.9 66.9 66.8 66.9 68.7
BLT${}_{\text{asr}}$ GUA 65.3 65.8 65.8 66.2 68.3

Table 5: GLEU performance when 10% of inputs are rejected (based on a metric),
and passed on to manual correction.

##  5 Conclusion

This work describes the application of ensemble distillation approaches to
sequence data. Two forms of distillation are discussed: standard distillation
where the distilled model predicts the ensemble mean; and distribution
distillation where the distribution over the ensemble predictions is modelled.
Though more challenging to train, ensemble distribution distillation yields
both an ensemble prediction and uncertainty measures associated with
prediction which allows, for example, more challenging predictions to be
manually corrected. The approaches were evaluated on grammatical error
correction tasks using data either matched to the written training data, or
mismatched speech data. Standard distillation was found to work well for
sequence distillation, but distribution distillation acting alone did not
yield good performance. By the combining the predictions from standard
distillation with the uncertainty predictions of distribution distillation,
both good performance and uncertainty measures could be obtained, with reduced
computational memory costs compared to the original ensemble.

Future work will examine: improved optimisation approaches and approximations
for sequence ensemble distribution distillation; and alternative uncertainty
criteria described in [15].

## References

  * [1]
  * [2] G. Hinton, O. Vinyals, & J. Dean, ”Distilling the Knowledge in a Neural Network,” Proc. NIPS, Montreal, Canada, 2014.
  * [3] G. Chen, W. Choi; X. Yu, T. Han, & M. Chandraker, ”Learning Efficient Object Detection Models with Knowledge Distillation, ” Proc. NeurIPS, Long Beach, CA, United States, 2017.
  * [4] J. Cui, B. Kingsbury, B. Ramabhadran, G. Saon, T. Sercu, K. Audhkhasi, A. Sethy, M. Nussbaum-Thom, & A. Rosenberg, ”Knowledge Distillation Across Ensembles of Multilingual Models for Low-resource Languages,” Proc. ICASSP, New Orleans, LA, United States, 2017.
  * [5] R. Yu, A. Li, V. I. Morariu, & L. S. Davis, ”Visual Relationship Detection with Internal and External Linguistic Knowledge Distillation,” Proc. ICCV, Venice, Italy, 2017.
  * [6] J. H. M. Wong & M. J. F. Gales, ”Sequence Student-teacher Training of Deep Neural Networks,” Proc. Interspeech, San Francisco, CA, United States, 2016.
  * [7] T. Asami, R. Masumura, Y. Yamaguchi, H. Masataki, & Y. Aono, ”Domain Adaptation of DNN Acoustic Models using Knowledge Distillation,” Proc. ICASSP, New Orleans, LA, United States, 2017.
  * [8] N. Papernot, P. McDaniel, X. Wu, S. Jha, & A. Swami, “Distillation as a Defense to Adversarial Perturbations Against Deep Neural Networks,” Proc. SP, San José, Costa Rica, 2016.
  * [9] Kevin P. Murphy. Machine Learning. The MIT Press, 2012.
  * [10] B. Lakshminarayanan, A. Pritzel, & C. Blundell, ”Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles,” Proc. NeurIPS, Long Beach, CA, United States, 2017.
  * [11] A. Malinin & M. J. F. Gales, ”Predictive Uncertainty Estimation via Prior Networks,” Proc. NeurIPS, Montreal, Canada, 2018.
  * [12] A. Malininm, B. Mlodozeniec, & M. J. F. Gales, ”Ensemble Distribution Distillation,” Proc. ICLR, 2020.
  * [13] X. Wu , K. M. Knill, M. J. F. Gales, & A. Malinin, “Ensemble Approaches for Uncertainty in Spoken Languagege Assessment,” Proc. Interspeech, Shanghai, China, 2020.
  * [14] Y. Kim & A. M. Rush, “Sequence-Level Knowledge Distillation,” Proc. EMNLP, Austin, TX, United States, 2016.
  * [15] A. Malinin & M. J. F. Gales, ”Uncertainty in Structured Prediction,” arXiv.org, 2002.07650v3, 2020.
  * [16] C. Zhou, J. Gu, & G. Neubig, ”Understanding Knowledge Distillation in Non-autoregressive Machine Translation,” Proc. ICLR, 2020.
  * [17] Z. Yuan & T. Briscoe, “Grammatical Error Correction Using Neural Machine Translation,” Proc. ACL, San Diego, CA, United States, 2016.
  * [18] R. Pang, T. N. Sainath, R. Prabhavalkar, S. Gupta, Y. Wu, S. Zhang, & C. Chiu, ”Compression of End-to-End Models”, Proc. Interspeech, Hyderabad, India, 2018.
  * [19] J. M. Hernandez-Lobato & R. Adams, ”Probabilistic Backpropagation for Scalable Learning of Bayesian Neural Networks,” Proc. ICML, Lille, France, 2015.

  * [20] Y. Gal & Z. Ghahramani, ”Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning,” Proc. ICML, New York City, NY, United States, 2016.
  * [21] L. Smith & Y. Gal, “Understanding Measures of Uncertainty for Adversarial Example Detection,” Proc. UAI, Monterey, CA, United States, 2018.
  * [22] Y. Gal, ”Uncertainty in Deep Learning”, Ph.D. thesis, University of Cambridge, 2016.
  * [23] J. Quiñonero-Candela, ”Dataset Shift in Machine Learning”, The MIT Press, 2009.
  * [24] M. Huang, Y. You, Z. Chen, Y. Qian, & K. Yu, ”Knowledge Distillation for Sequence Model,” Proc. Interspeech, Hyderabad, India, 2018.
  * [25] Y. Lu, M. J. F. Gales, & Y. Wang, “Spoken Language ‘Grammatical Error Correction’,” Proc. Interspeech, Shanghai, China, 2020.
  * [26] D. Nicholls, “The Cambridge Learner Corpus: Error Coding and Analysis for Lexicography and ELT,” Proc. Corpus Linguistics, Lancaster, United Kingdom, 2003.
  * [27] E. Izumi, K. Uchimoto, & H. Isahara, “The NICT JLE Corpus Exploiting the Language Learners’ Speech Database for Research and Education,” Proc. IJCIM, Bangkok, Thailand, 2004.
  * [28] L. Chambers & K. Ingham, “The BULATS Online Speaking Test,” Research Notes, vol. 43, pp. 21–25, 2011. [Online]. Available: http://www.cambridgeenglish.org/images/23161-research-notes-43.pdf
  * [29] K. M. Knill, M. J. F. Gales, P. P. Manakul, & A. Caines, ”Automatic Grammatical Error Detection of Non-native Spoken Learner English,” Proc. ICASSP, Brighton, United Kingdom, 2019.
  * [30] Y. Lu, M. J. F. Gales, K. M. Knill, P. P. Manakul, L. Wang, & Y. Wang, ”Impact of ASR performance on spoken grammatical error detection,” Proc. Interspeech, Graz, Austria, 2019.
  * [31] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, & Ł. Kaiser, ”Attention is All You Need,” Proc. NeurIPS, Long Beach, CA, United States, 2017.
  * [32] W. Maddox, T. Garipov, P. Izmailov, D. P. Vetrov, & A. G. Wilson, ”A Simple Baseline for Bayesian Uncertainty in Deep Learning,” Proc. NIPS, Vancouver, Canada, 2019.
  * [33] C. Napoles, K. Sakaguchi, M. Post, & J. Tetreault, ”Ground truth for grammatical error correction metrics,” Proc. ACL, Beijing, China, 2015.
  * [34] Y. Wu, M. Schuster, Z. Chen, Q. V. Le, et al. ”Google’s Neural Machine Translation System: Bridging the Gap between Human and Machine Translation,” arXiv.org, 1609.08144v2, 2016.
  * [35] A. Malinin, A. Ragni, K. M. Knill, & M. J. F. Gales, ”Incorporating Uncertainty into Deep Learning for Spoken Language Assessment, ” Proc. ACL, Vancouver, Canada, 2017.

novel approach to distribution distilling sequence models, and additionally
presented a second improved objective function. Since distribution
distillation can be cumbersome, a less restricted approach, where one
distribution distills onto two models has also been introduced. The use of a
separate uncertainty model significantly improves upon the single model
distribution distillation approach, and additionally has the added benefit
that it can perform as well as any other available model. Furthermore, this
distribution distillation approach was also shown to perform similarly to an
ensemble when used in a hybrid Oracle-DNN system, showing that distribution
distilled models can perform as well as their teacher ensembles and replace
them to save resources. Future studies can build on top of this by focusing on
obtaining well performing single model distribution distilled models, and
performing sequence-level distribution distillation as opposed to word- level
performed here.


# Harmonic ECME Excited by Energetic Electrons Travelling Inside A Coronal
Loop

Mehdi Yousefzadeh  Institute of Space Sciences, Institute of Frontier and
Interdisciplinary Science, Shandong University, Shandong, China.
Hao Ning  Institute of Space Sciences, Institute of Frontier and
Interdisciplinary Science, Shandong University, Shandong, China.
Yao Chen  Institute of Space Sciences, Institute of Frontier and
Interdisciplinary Science, Shandong University, Shandong, China.
Yao Chen<EMAIL_ADDRESS>
(Received ***; Revised ***; Accepted by APJ - 04-Jan-2021)

###### Abstract

A complete understanding of solar radio bursts requires developing numerical
techniques which can connect large-scale activities with kinetic plasma
processes. As a starting point, this study presents a numerical scheme
combining three different techniques: (1) extrapolation of magnetic field
overlying a specific active region in order to derive the background field,
(2) guiding-center simulation of dynamics of millions of particles within a
selected loop to reveal the integral velocity distribution function (VDF)
around certain sections of the loop, and (3) particle-in-cell (PIC) simulation
of kinetic instabilities driven by energetic electrons initiated by the
obtained distributions. Scattering effects at various levels (weak, moderate,
and strong) due to wave$\backslash$turbulence-particle interaction are
considered using prescribed time scales of scattering. It was found that the
obtained VDFs contain strip-like and loss-cone features with positive
gradient, and both features are capable of driving electron cyclotron maser
emission (ECME), which is a viable radiation mechanism for some solar radio
bursts, in particular, solar radio spikes. The strip-like feature is important
in driving the harmonic X mode, while the loss-cone feature can be important
in driving the fundamental X mode. In the weak-scattering case, the rate of
energy conversion from energetic electrons to X2 can reach up to $\sim
2.9\times 10^{-3}$$E_{k_{0}}$, where $E_{k_{0}}$ is the initial kinetic energy
of energetic electrons. The study demonstrates a novel way of exciting X2 mode
in the corona during solar flares, and provides new sight into how escaping
radiation can be generated within a coronal loop during solar flares.

Solar corona (1483); Solar activity (1475); Radio bursts (1339); Solar coronal
radio emission (1993); Plasma astrophysics (1261)

††thanks: E-mail<EMAIL_ADDRESS>
##  1 Introduction

Solar radio bursts represent sudden enhancement of emission at radio
wavelengths released from the solar atmosphere. Various types of bursts have
been classified according to their manifestation on the dynamic spectrum, such
as bursts of type I to V, millisecond solar spikes, flare decimetric
continuum, etc (Wild et al., 1963; McLean & Labrum, 1985; see, e.g., Feng et
al., 2012, Li et al., 2017, Chen et al., 2014, and Vasanth et al., 2019 for
latest studies). Among them, solar spikes are radio bursts closely associated
with the impulsive stage of solar flares. Spikes are characterized by their
extremely high brightness temperature of up to 1015 K, the highest among all
types of solar radio bursts. They appear in metric-decimetric wavelength,
exhibiting very narrow-band spectral width with rapidly-rising intensity
profile and extremely short characteristic time scales ($\sim$ a few to tens
of ms for individual burst). Thousands of spikes can exist in an event. This
leads some authors to suggest that solar spikes represent elementary energy
release events during solar flares (see, e.g., Benz, 1986; Benz & Kane, 1986;
Fleishman & Mel’nikov, 1998).

In earlier studies, ECME is considered to be the emission mechanism of solar
spikes (see, e.g., Holman et al., 1980; Melrose & Dulk, 1982; Robinson, 1991).
The ECME mechanism was initially proposed by Twiss (1958), Gaponov (1959) and
Schneider (1959). Wu & Lee (1979) realized the importance of relativistic
effect in the resonance condition, and applied ECME to Decametric Radiation
(DAM) of Jupiter. Later, the theory was used to explain solar radio bursts,
such as solar spikes (Melrose & Dulk, 1982; Sharma & Vlahos, 1984; Aschwanden,
1990) and type IIIs (Wu et al., 2005). According to Wu & Lee (1979) and
follow-up studies, ECME can be excited by energetic electrons with a positive
gradient along the perpendicular direction in the velocity space, i.e.,
$\partial f/\partial{v_{\perp}}>0$, in plasmas with
$\omega_{pe}/\Omega_{ce}<1$. In addition, ECME favors the growth of X mode. In
plasmas with $\omega_{pe}/\Omega_{ce}>1$ the growth rates of X and O modes
decline rapidly while Z (or the electrostatic upper hybrid mode at large wave
number) and W become the dominant modes (Sharma & Vlahos, 1984; Lee et al.,
2009; Yi et al., 2013; Ni et al., 2020).

Most earlier studies suggest the solar spikes are excited by ECME of energetic
electrons with the loss-cone distribution, i.e., via the loss-cone maser. Such
distribution can form when electrons are trapped and bounced within a coronal
loop. Nevertheless, unresolved issues exist that hinder us from a complete
understanding of the exact emission mechanism and how one can use radio data
to infer coronal properties. One major difficulty lies in the fact that any
solar radio bursts are a consequence of multi-scale process. An event would
start from large-scale magnetic eruptive activities in the astronomical to MHD
scale, followed by particle acceleration via magnetic reconnection or shocks
in the particle kinetic scale, and excitation and propagation of radiations.
Present computational resources are far from enough to simulate such cross-
scale physics.

Another difficulty is that the loss-cone maser for solar spikes mainly excite
the fundamental X mode (X1), and it is difficult for X1 to pass through the
second harmonic absorption layer during its escape from the source (Melrose &
Dulk, 1982). This gives the escaping difficulty of ECME when being applied to
solar radio bursts. To resolve this, Melrose & Wheatland (2016) proposed that
a low-density cavity exists in the source region, which can duct the X mode
out. Existence of such low-density cavity in the corona has not been verified
yet. Another option is to have direct excitation of X at higher harmonics such
as X2 and X3, since the absorption rates at the third or fourth harmonic layer
are considerably less in comparison to that at the second harmonic. Yet,
direct and efficient excitation of X2 and X3 are rarely reported.

Most, if not all, existing studies on emission mechanism of solar radio
bursts, including solar spikes, started from analytically-prescribed
distribution function of energetic electrons. Such function is not determined
by the particle dynamics in a self-consistent way. In this study, we present a
preliminary effort to connect macro-scale dynamics of energetic electrons with
micro-scale kinetic instabilities energized by these electrons to explain the
possible emission mechanisms behind solar spikes. This is done by combining
three sets of numerical techniques working at different scales. It was found
that the VDFs deduced from the simulations contain certain interesting
features that are capable of exciting emission efficiently at the second
harmonic (X2) via ECME, thus providing a novel approach to resolve the
escaping difficulty of ECME. The following section introduces numerical
techniques used here, the third section presents analysis and major results,
summary and discussion are given in the last section.

##  2 Numerical Schemes

To bridge the physics of different spatial scales and simplify the effort as
much as possible, we combine three sets of numerical techniques, including (1)
the nonlinear force-free field (NLFFF) extrapolation method (Wiegelmann, 2004
and Wiegelmann et al., 2006) to describe large-scale topology of corona
magnetic field within which electrons are traveling and trapped, (2) the
guiding-center (GC) simplification method (e.g., Northrop, 1963; Gordovskyy et
al., 2010) to simulate dynamical motions of a large-number of electrons so as
to infer evolution of their VDF, and (3) the particle-in-cell (PIC) simulation
fed by the obtained VDFs to explore the excited kinetic instabilities and
radiations. This represents a significantly-simplified approach for the
proposed study.

###  2.1 NLFFF Extrapolation of the Coronal Magnetic Field

To start the simulation, we selected one specific active region (NOAA
AR11283), which released a group of eruptions including three big flares
(M5.3, X2.1, and X1.8) on September 6 and 7, 2011. It has been investigated by
several authors (e.g., Ruan et al., 2014; Liu et al., 2014; Romano et al.,
2015). Here we start from the magnetic condition before the X2.1 flare. The
NLFFF extrapolation was based on the magnetogram data observed at 21:30 UT by
the Helioseismic and Magnetic Imager (HMI; Schou et al., 2012).

The HMI data and the result of extrapolation are shown in Figure 1(a), from
which we observe a clear twisted flux rope structure whose eruption releases
the X2.1 flare. Three representative field lines with different amount of
twist are plotted. The loop shown in panel (b) contains no significant twist,
the field strengths at the two foot-points are close to 1800 G and 1200 G,
respectively. The location with the minimum field strength is close to the
loop top. The loop in panel (c) presents the M-like topology with significant
twist, the field strengths at the two foot-points are nearly 300 G and 1800 G,
respectively. The location with the minimum field strength is located at the
right foot-point. The field line in panel (d) has more twists, the field
strength at its left foot is 1400 G, decreasing to a local minimum at 800 G
and then increasing to the maximum at nearly 1300 G around the top part, then
decreasing to 400 G and further increasing to 1200 G at the right foot.

As a starting point, we conduct simulation of particle dynamics within the
simple untwisted loop configuration, i.e., that shown in panel (b). Note that
selection of the specific AR and the loop is just to provide a realistic
magnetic condition of the corona so as to initiate the simulation. The loop
forms a simple magnetic trap (or mirror) and represents a typical field line
in the background corona overlying the pre-eruption flux rope structure. Using
other twisted configuration would be interesting yet more complex, and the
essential physics should be revealed by the study using the simple untwisted
loop structure. From the extrapolation technique we can obtain three-
dimensional distribution of field strength and direction along the whole loop.
The information will be fed to the following GC simulation of particle
dynamics as input. The time scale of the bouncing motion can be estimated with
the loop length and the average energy of electrons. Here, the height of the
selected loop structure is about 100 arcsec, and the length is nearly 300
arcsec. Thus, it takes $\sim 3$s for electrons with an energy of 15 keV to
travel such distance.

###  2.2 Guiding-center simulation of dynamics of energetic electrons

VDF of energetic electrons accounting for radio emission remains unknown. It
is related to the acceleration and heating mechanism and further wave-particle
interaction during solar flares. In any case, strong heating seems to be an
intrinsic part of solar flares, therefore one may obtain heated electrons with
high-temperature Maxwellian-like distribution. The heated electrons will be
injected into a loop which might be the post-reconnection loop or one of the
loops directly connected to the reconnection region. If reconnection takes
place high in the corona, as indicated by X-ray data of flares (Masuda et al.,
1994; Sui & Holman, 2003; Chen et al., 2017), top of loops could be a probable
place of injection. To simulate the following dynamic motion, we inject two
million of electrons with a Maxwellian distribution around the loop top in an
impulsive way (labelled with “I” in Figure 1(b)). The length of the injection
section is $\sim$ 5 arcsec. The VDF is expressed by

$f_{0}=\frac{1}{(2\pi)^{3/2}v_{0}^{3}}\exp\left(-\frac{u^{2}}{2v_{0}^{2}}\right)$
(1)

where $v_{0}=0.24$ c and c is the speed of light. The GC simplification is
described with standard relativistic equations (Northrop, 1963; Gordovskyy et
al., 2010). To apply this method, the characteristic time and spatial scales
of particle gyro-motion should be much less than those of the magnetic field
variation. For the simple untwisted large-scale loop structure, such
conditions are well satisfied. The classical Bulirsch-Stoer approach is
employed to solve the GC equations, which is characterized by a high level of
accuracy in energy conservation.

Pre-existing waves and those excited after injection can scatter electrons and
modify their VDFs. To simulate such effect, we introduce a parameter
$\tau_{s}$ to represent the time scale of pitch-angle scattering: within each
interval of $\tau_{s}$, pitch angle of an electron will be changed to a random
value in a range of [0, $\pi$]. A larger $\tau_{s}$ represents weaker
scattering and vice versa. We refer to Piana et al. (2007) and Chen &
Petrosian (2013) for characteristic values of $\tau_{s}$ during solar flares.
For electrons of few tens of Kev, $\tau_{s}$ could be around 0.1s to 1s. Thus,
we set $\tau_{s}$ to be 5s, 1s, and 0.5s, to represent weak (Case W), moderate
(Case M), and strong (Case S) levels of scattering, respectively.

###  2.3 VPIC simulation of ECME driven by energetic electrons

We use the Vector-PIC (VPIC) open-source code released by the Los Alamos
National Labs for PIC simulation. The code employs a second-order, explicit,
leapfrog algorithm to update particle positions and velocities, along with a
second-order finite-difference time-domain solver for evolution of electric
and magnetic fields which are described by the full Maxwell equations (Bowers
et al., 2008a,2008b_Bowers, 2008b, 2009).

To initiate the PIC simulation, we adopt Maxwellian distribution for
background plasmas consisting of electrons and protons. The background
magnetic field is set to be along $\vec{e}_{z}$
($\vec{B}_{0}=B_{0}\vec{e}_{z}$) and the wave vector ($\vec{k}$) is within the
$xOz$ plane. VDFs of energetic electrons are taken from the above GC
simulation, and only those at the loop top are used. According to the NLFFF
extrapolation the field strength around the loop top is $\sim$400 G, the
number density of background plasmas $n_{0}$ is assumed to be $10^{9}$ cm-3,
thus the ratio of plasma oscillation frequency and electron gyro-frequency,
i.e., $\omega_{pe}/\Omega_{ce}$, is 0.25. This means that ECME should be
important. The density ratio of energetic and background electrons
($n_{e}/n_{0}$) is taken to be 0.05. The grid spacing is
$\Delta=2.7\,\lambda_{D}$ where $\lambda_{D}$ is the Debye length of
background electrons and the grid number is [512, 512]. The simulation time
and domain are 2000 $\omega^{-1}_{pe}$ for Cases W and M, and 2500
$\omega^{-1}_{pe}$ for Case S, and the simulation domain in space is given by
$Lx=Lz=25\,c/\omega_{pe}$. The resolved ranges of wave number and frequency
are then [-16, 16] $\Omega_{ce}/c$ and [0, 3.2] $\Omega_{ce}$, respectively.
In each cell for each species, we employ 1000 macroparticles. The charge
neutrality condition is maintained.

Note that VDFs of energetic electrons within the loop change with many
factors, such as time, location, form of initial distribution, location of
injection, scattering effect, etc. This results in lots of variations. Keeping
this in mind, the present study is taken as a starting point along the
proposed line of research, with two main aims: (1) how VDFs evolve along the
loop structure while considering various levels of scattering, and (2) whether
and how VDFs are capable of efficient excitation of harmonic X modes (X2, as
well as other wave modes) so as to partially address the escaping difficulty
associated with ECME. In the following, we first present evolution of VDFs
within four different sections of the loop, then compare outputs of PIC
simulations that are initiated by three sets of VDF obtained at the same
location (loop top) and time (1s after injection) with different levels of
scattering. Results with VDFs at other location and time, as well as using
different forms of initial distribution within a different loop structure will
be examined in future.

##  3 Numerical Results and Analysis

###  3.1 VDFs given by the GC simulation

In our GC simulation, motions of two millions of electrons are resolved. To
evaluate evolutions of VDFs along the loop, we select four sections, referred
to as A, B, C, and D, respectively (see Figure 1b). Their lengths are taken to
be large enough ($\sim$12 arcsec) so as to get sufficient number of electrons.
The VDFs at certain moment are obtained by collecting the information of all
electrons within the specific section at the moment of interest.

Note that electrons with small initial pitch angles may get lost from the
loop. Electrons that are bounced backwards give rise to the well-known loss-
cone distribution. If electrons are also mirrored from the other side of the
loop, a double-sided loss-cone distribution is obtained. For the selected
loop, the ratio of maximum to minimum $B$, i.e., the mirror ratio, is $\sim$
3, this gives a loss cone angle of $\sim$ 35∘. When scattering effect is
included, all electrons will get lost eventually with stronger scattering
leading to faster loss. These deductions are verified by our GC simulations.

Figure 2(a-d) present temporal evolution of the obtained VDFs at different
loop sections (A-D) for Case W with $\tau=5s$. The double-sided loss cone
appears as an overall morphology of VDFs within the first few seconds. Angle
of the loss cone and timing of its development are consistent with the above
estimate. The loss-cone angle increases from Section A to Section D, as a
result of the decrease of mirror ratio. As seen from the accompanying movie
for section A (around the top), the initial VDF is Maxwellian since this
section partially overlaps with the injection region while lower sections
include no electrons initially.

In addition to the loss-cone feature, there appear several strip-like features
that exhibit positive gradient. The lower strip feature appears shortly after
the injection, corresponding to particles with relatively large pitch angles
and therefore small bouncing distance; the upper one appears with large
$v_{\perp}$ and large-positive $v_{\parallel}$, then mitigating to regions
with smaller $v_{\perp}$ and negative $v_{\parallel}$. The strips appear
earlier for sections closer to the top. They are given by energetic electrons
that have experienced bouncing by at least one time and arrive at the section
window at the specific time. Similar multi-strip features have been
demonstrated earlier by White et al. (1983) in their analytical evaluation of
VDFs for a one-dimensional magnetic field line. We have repeated their
calculation and obtain similar results, confirming the validity of the
simulation result presented here. The features are of major interest here
since they emerge along the whole structure with significant positive
gradients of VDFs, this makes them to be candidates of efficient driver of
ECME.

We also conducted GC simulations for cases with moderate (Case M) and strong
(Case S) levels of scattering. The solutions within Section A at the same time
after the injection ($t=1s$) are presented in Figures 2e and 2f (and the
accompanying movie). They are basically similar to that for Case W. With
stronger scattering effect, electrons are lost from the loop at a faster pace,
the strip feature becomes less sharp, more blurred, and number of observable
strips may decrease. For example, only very weak signature of strips remains
in Case S (see Figure 2f).

###  3.2 PIC simulations of ECME radiation

The PIC simulations are initiated by the three sets of VDFs of energetic
electrons that correspond to different levels of scattering, referred to as
Cases W, M, and S, respectively. The VDFs are obtained at 1s after injection
and have been shown in Figures 2(a, e, f). The parameter setup of the PIC
simulation has been presented above. In the following texts, we first
introduce results for Case W, and then compare them with other cases.

####  3.2.1 PIC simulation for Case W

We start from Case W which corresponds to weak level of scattering with $\tau$
= 5s. The PIC-evolved distribution has been presented in Figure 3a ($t\sim
1000\,\omega_{pe}^{-1}$) and the accompanying movie. After $t\sim
500\,\omega_{pe}^{-1}$, energetic electrons start to diffuse significantly and
the strip-feature gets blurred, indicating strong wave-particle interaction.

In Figure 4a we present the $\vec{k}$-space intensity distribution of the
strongest wave mode at the corresponding wave vector $\vec{k}$, the
$\omega$-$k$ dispersion analysis is presented in Figure 5a and the online
animation. We see that the strongest mode is the quasi-perpendicular
($90^{\circ}<\theta_{B}<110^{\circ}$) Z mode with frequency slightly less than
$\Omega_{ce}$, while the harmonic X mode (X2) is also at a high intensity
along the quasi-perpendicular direction ($95^{\circ}<\theta_{B}<105^{\circ}$)
with frequency slightly less than $2\,\Omega_{ce}$. The X1 emission is very
weak and mainly around $30^{\circ}$. Energy curves of these modes have been
plotted in Figure 3d. From the energy curves, X2 and Z modes first grow
linearly at basically the same growth rates which are fitted to be 0.033 and
0.026 $\omega_{pe}$, respectively (see Figure 3d for fitting lines); before
$t\sim 500\,\omega^{-1}_{pe}$, the two modes reach the maximum intensities of
$\sim 2.9\times 10^{-3}$ and $\sim 4.8\times 10^{-3}$ $E_{k_{0}}$,
respectively, where $E_{k_{0}}$ represents the total kinetic energy of
energetic electrons.

The VDFs exhibit two types of features with positive gradients, i.e., the
loss-cone and the strip-like feature. To figure out the driving agency of each
mode, we have plotted resonance curves onto the corresponding VDF map as shown
in Figure 2a. The curves are for the specific parameter set of
($\omega,k,\theta$) at which the strongest intensities for X2 and Z modes are
achieved, while for X1 we select a representative propagation angle of 35∘. It
is clear that the resonance curves of X2 and Z modes almost overlap with each
other, and both modes are excited by the same upper strip-like feature through
the electron cyclotron maser instability (Wu & Lee, 1979). The lower strip-
like feature passes through the background Maxwellian distribution, therefore
unable to excite any wave mode. On the other hand, the resonance curve of X1
crosses the loss-cone feature, results in a weak excitation.

In summary, for Case W we get efficient excitation of Z and X2, yet without
significant excitation of X1, both Z and X2 are associated with the upper
strip-like feature of the VDF, while the loss-cone feature only leads to weak
excitation of X1.

####  3.2.2 Comparison of PIC simulations with different levels of scattering

In this subsection, results of PIC simulations for Cases M ($\tau=1s$) and S
($\tau=0.5s$) are presented, and differences of the three cases are
highlighted. The initial and PIC-evolved VDFs have been presented in Figures
2(e-f) and 3(b-c). The corresponding $\vec{k}$-space intensity distribution of
waves within time range of $1500-2000\,\omega_{pe}^{-1}$ for Cases S $\&$ M
have been presented in Figure 4(b-c), results of $\omega$-$k$ dispersion
analysis have been presented in Figure 5(b-c) and the online animation, and
energy curves for various wave modes have been presented in Figure 3(d-f).

In general, it takes longer time for modes to saturate in cases with stronger
scattering. As seen from the energy curves (Figure 3), in Case M the duration
of the linear stages of X2 and Z are both $\sim 800\,\omega_{pe}^{-1}$, longer
than the corresponding duration of the two modes in Case W ($\sim
500\,\omega_{pe}^{-1}$), and shorter than the duration of the linear growth of
X1 in Case S which is $\sim 1800\,\omega_{pe}^{-1}$.

In Case M, intensities of both X2 ($\sim 3.0\times 10^{-5}\,E_{k_{0}}$) and Z
($\sim 3.3\times 10^{-4}\,E_{k_{0}}$) are much smaller than those of Case W,
and X1 ($\sim 2.1\times 10^{-5}\,E_{k_{0}}$) becomes considerably enhanced
over the noise. Ranges of frequency and wave number of X2 ($\sim
1.92-1.95\,\Omega_{ce}$) and Z ($\sim 0.96-0.98\,\Omega_{ce}$) modes are quite
close to those of Case W. In Case S, X1 becomes the strongest, reaching a
maximum energy of $\sim 3.8\times 10^{-4}\,E_{k_{0}}$, while both Z and X1
have energies of $\sim 2\times 10^{-7}\,E_{k_{0}}$) much weaker than those in
other cases. X1 has a symmetric arc-like emission pattern along quasi-parallel
to oblique direction ($0^{\circ}<\theta_{B}<55^{\circ}$), its intensity is
close to the thermal noise level for larger $\theta_{B}$. As read from the
VDFs shown in Figures 2 and 3, densities of energetic electrons along the edge
of the loss cone in Case S are much higher that those in Case W. This is due
to different level of scattering. With stronger scattering the loss cone
feature is better-developed within shorter period and the strip-like features
are less significant, while for weaker scattering, a significant part of
energetic electrons is still bouncing and the VDF is characterized by strip-
like features at the time of interest (1s after injection). It is apparent
that the upper strip-like feature is more efficient (than the loss-cone
feature) in converting electron energy into wave modes. This explains why the
loss-cone feature in Case W and Case S plays very different role and why the
duration of linear stage in weaker-scattering case lasts shorter.

Resonance curves shown in Figures 2e and 2f are plotted with parameters
provided in the figure caption. Again, Z and X2 modes (if exist) are excited
by the upper strip-feature in Cases M $\&$ W, and the very-strong emission of
X1 in Case S is driven by the loss-cone feature. The linear growth rate of X1
in Case S is much smaller than that of X2 and Z in Cases W and M. Note that
similar emission pattern of X1 via the loss-cone maser has been reported by
earlier studies (see, e.g., Wagner et al., 1983,1984_Wagner, 1984; Yoon &
Ziebell, 1995).

##  4 Summary and Discussion

As a starting point to bridge the large-scale dynamics and small-scale plasma
kinetic maser instabilities associated with solar radio bursts, in particular,
solar spikes, we developed a numerical scheme combining techniques including
the magnetic field extrapolation to describe the magnetic configuration of a
normal loop, the guiding-center method to infer the temporal evolution of VDFs
along various sections of the loop while taking the effect of pitch-angle
scattering into account, and PIC simulations to further explore the kinetic
instabilities driven by electrons with the obtained VDFs. Consistent with
earlier studies, VDFs of energetic electrons that are released from the loop
top manifest interesting strip-like feature with significant positive velocity
gradient, together with the well-known loss-cone feature. According to further
PIC simulations, the strip-like feature is essential for efficient excitation
of ECME at harmonic X mode (in the weak-scattering case), while the loss-cone
feature can be efficient in exciting the fundamental X mode (in the strong-
scattering case). Efficient amplification of X2 favors the escape of ECME
radiation from the corona, this effectively reduces the limitation of applying
ECME to solar spikes. The study provides new insight into how escaping
radiation, in particular, the harmonic X mode, can be generated within a
coronal loop.

In most earlier studies relevant to ECME and solar radio spikes (e.g., see
Aschwanden, 1990), loss-cone distributions have been applied. This
distribution may lead to efficient excitation of X1 mode, while being unable
to drive efficient excitation of X2 or higher harmonics. During its outward
escape, it is difficult for the X1 mode to pass the second-harmonic absorption
layer where the absorption coefficient is notoriously large. This gives rise
to the escaping difficulty of fundamental emission from the corona. To resolve
this issue, efficient excitation of radio bursts at higher harmonics is
required. Thus, the most significant point of this study is the efficient
excitation of X2 under certain conditions. Here we demonstrate that strip-like
features of VDFs can form due to the initial bouncing motion of energetic
electrons that are released at the loop top. These features are formed during
the early stage of the VDF relaxation towards a well-developed loss-cone
distribution, and they are critical to efficient excitation of harmonic X
mode.

It should be pointed out that excitation of ECME radiation depends sensitively
on many conditions. The effect of pitch-angle scattering on both VDFs and
further excitation of wave modes has been demonstrated. In addition, initial
distribution of electrons, injection location and the way of injection (e.g.,
impulsively or continuously, see White et al., 1983), ratio of the
characteristic frequencies ($\omega_{pe}/\Omega_{ce}$), all these factors
shall affect the obtained form of VDFs and further excitation of various wave
modes. Future study shall expand the parameter regime for a more-complete
understanding of the multi-scale radiation processes.

## Acknowledgments

The present study is supported by the National Natural Science Foundation of
China (11790303 (11790300), 11750110424, and 11873036). The authors
acknowledge the open-source Vector Particle In Cell (VPIC) code provided by
Los Alamos National Labs (LANL), Super Cloud Computing Center (BSCC, URL:
http://www.blsc.cn/) for providing HPC resources, and Dr. Alexander William
Degeling and Xiangliang Kong (SDU) $\&$ Dr. Mahboub Hosseinpour (TU) for
helpful discussion.

## References

  * Aschwanden et al. (1990) Aschwanden, M. J., Benz, A. O., Schwartz, R. A., Lin, R., Pelling, R. M., & Stehling, W. 1990, Sol. Phys., 130, 39
  * Aschwanden (1990) Aschwanden, M. J. 1990, A&A, 237, 512
  * Benz (1986) Benz, A. O. 1986, Sol. Phys., 104, 99
  * Benz & Kane (1986) Benz, A. O., & Kane, S. R. 1986, Sol. Phys., 104, 179
  * Bowers et al. (2008a) Bowers, K. J., Albright, B. J., Bergen, B., et al. 2008a, IEEE Press, 63, 11
  * (6) Bowers, K. J., Albright, B. J., Yin, L., et al. 2008b, Physics of Plasmas, 15, 055703
  * (7) Bowers, K. J., Albright, B. J., Yin, L., et al. 2009, Journal of Physics Conference Series, 012055
  * Chen & Petrosian (2013) Chen, Q., & Petrosian, V. 2013, ApJ, 777, 33
  * Chen et al. (2014) Chen, Y., Du, G., Feng, L., et al. 2014, ApJ, 787, 59
  * Chen et al. (2017) Chen Y., Wu Z., Liu W. et al. 2017, ApJ, 843, 8
  * Feng et al. (2012) Feng, S. W., Chen, Y., Kong, X. L., et al. 2012, ApJ, 753, 21
  * Feng et al. (2013) Feng, S. W., Chen, Y., Kong, X. L., et al. 2013, ApJ, 767, 29
  * Fleishman & Mel’nikov (1998) Fleishman G. D. & Mel’nikov V. F. 1998, PhyU, 41, 1157
  * Gaponov (1959) Gaponov, A. V. 1959, Izv VUZ Radiofizika , 2, 450-462
  * Gordovskyy et al. (2010) Gordovskyy, M., Browning, P. K., & Vekstein, G. E. 2010a, A&A, 519,A21
  * Holman et al. (1980) Holman, G. D., Eichler, D., & Kundu, M. R. 1980, IAU Symp. 86, Radio Physics of the Sun, ed. M. R. Kundu & T. E. Gergely (Cambridge: Cambridge Univ. Press), 457
  * Lee et al. (2009) Lee, K. H., Omura, Y., Lee, L. C., et al. 2009, Phys. Rev. Lett., 103, 105101
  * Li et al. (2017) Li, C. Y., Chen, Y., Wang, B., et al. 2017, Sol. Phys., 292, 82
  * Liu et al. (2014) Liu, C., Deng, N., Lee, J., et al. 2014, ApJ, 795, 128
  * Masuda et al. (1994) Masuda, S., Kosugi, T., Hara, H., Tsuneta, S., & Ogawara, Y. 1994, Natur, 371, 495
  * McLean & Labrum (1985) McLean, D. J. & Labrum, N. R. 1985, Solar Radiophysics: Studies of Emission from the Sun at Metre Wavelengths
  * Melrose & Dulk (1982) Melrose, D. B. & Dulk, G. A. 1982, ApJ, 259, 844
  * Melrose et al. (1986) Melrose, D. B., Dulk, G. A., Cairns, I. H. 1986, A&A, 163, 229
  * Melrose & Wheatland (2016) Melrose, D. B. & Wheatland, M. S. 2016, Sol. Phys., 291, pp. 3637, 658
  * Northrop (1963) Northrop, T. 1963, The Adiabatic Motion of Charged Particles (New York: Interscience)
  * Ni et al. (2020) Ni, S., Chen, Y. et al. 2020, ApJLett., V891,1
  * Piana et al. (2007) Piana, M., Massone, A. M., Hurford, G. J., et al. 2007, ApJ, 665, 846
  * Robinson (1991) Robinson, P. A. 1991, Sol. Phys., 134, 299
  * Romano et al. (2015) Romano P., Zuccarello F., Guglielmino S. L., et al. 2015, A&A, 582, 55
  * Ruan et al. (2014) Ruan, G., Chen, Y., Wang, S., et al. 2014, ApJ, 784, 165
  * Schneider (1959) Schneider, J. 1959, Phys. Rev. Lett., 2, 504
  * Schou et al. (2012) Schou, J., et al. 2012, Sol. Phys., 275, 229
  * Sharma & Vlahos (1984) Sharma, R. R. & Vlahos, L. 1984, BAAS16, 536
  * Sui & Holman (2003) Sui, L., & Holman, G. D. 2003, ApJ, 596, L251
  * Twiss (1958) Twiss, R. Q. 1958, Australian Journal of Physics, 11, 564
  * Vasanth et al. (2019) Vasanth, V., Chen, Y., Lv, M., et al. 2019, ApJ, 870, 30
  * Wagner et al. (1983) Wagner, J. S., Lee, L. C., Wu, C. S., et al. 1983, J. Geophys. Res. Lett., V10(6), 483
  * (38) Wagner, J. S., Lee, L. C., Wu, C. S., et al. 1984, Radio Science, 19, 509
  * Wheatland et al. (2000) Wheatland, M. S., Sturrock, P. A., & Roumeliotis, G. 2000, ApJ, 540, 1150
  * White et al. (1983) White, S. M., Melrose, D. B., $\&$ Dulk, G. A., Proc. Astron. Soc. Australia. 1983, 5 (2), 188
  * Wild et al. (1963) Wild, J. P., Smerd, S. F., & Weiss, A. A. 1963, ARA&A, 1, 291
  * Wiegelmann (2004) Wiegelmann, T. 2004, Sol. Phys., 219, 87
  * Wiegelmann et al. (2006) Wiegelmann, T., Inhester, B., & Sakurai, T. 2006, Sol. Phys., 233, 215
  * Wu & Lee (1979) Wu, C. S. & Lee, L. C. 1979, ApJ, 230, 621
  * Wu et al. (2005) Wu, C. S., Wang, C. B., Zhou, G. C., Wang, S., & Yoon, P. H. 2005, ApJ, 621, 1129
  * Yi et al. (2013) Yi, S., Lee, S.-Y., Kim, H.-E., et al. 2013, J. Geophys. Res., 118, 7584
  * Yoon & Ziebell (1995) Yoon, P. H., & Ziebell L.F. 1995, Phys. Rev. E 51, 4908-4916

Figure 1: (a) The magnetogram of HMI for AR 11283 before the eruption of the
major X2.1 flare on September 6, 2011, over-plotted with field lines given by
NLFFF extrapolation; (b-d) representative lines with different amount of
twist. The color in panels (b) represents the field strength of lines
presented in (b-d). The letters “A-D” represent the sections within which VDFs
will be examined, and “I” represents the region of injection. Figure 2: (a-d)
The obtained VDFs at a specific time (1s after injection) within different
loop sections (see Figure 1b) for the case with weak level of scattering,
i.e., Case W with $\tau=5s$. Numbers of particles ($N$) used to get the
integral VDFs are written on top of each panel. (e-f) VDFs corresponding to
moderate (M), and strong (S) levels of scattering for Section A, also obtained
at 1s after injection. Resonance curves corresponding to the X2, X1, and Z
modes are over-plotted. For X2, parameters used to plot these curves are (the
wave frequency $\omega=1.98\,\Omega_{ce}$, the total wave number
$k=1.92\,\Omega_{ce}/c$, the propagation angle $\theta=100^{\circ}$, and the
harmonic number $n=2$); for X1 they are ($1.08\,\Omega_{ce}$,
$0.65\,\Omega_{ce}/c$, $35^{\circ}$, and 1), for Z they are
($0.96\,\Omega_{ce}$, $1.15\,\Omega_{ce}/c$, $100^{\circ}$, and 1). The
animation of this figure (panels (a), (e) and (f)) for the main three cases
(W, M and S) begins at t = 0s and advances 0.1s at a time up to t = 3s then
continues with 2s interval till the end, at t $\sim$ 19s. The realtime
duration of the video is 7s.
(An animation of this figure is available). Figure 3: Upper panels: VDFs at
$t=1000\,\omega_{pe}^{-1}$ obtained by the PIC simulation for Cases W, M, and
S; lower panels: Temporal profiles of energies of various wave modes (X2, X1,
and Z), normalized to the total energy of energetic electrons ($E_{k_{0}}$).
Dashed lines represent exponential fitting of linear growth rates, which are
(0.033, 0.0016, 0.026) $\omega_{pe}^{-1}$ for X2, X1, and Z in Case W, and
(0.01, 0.0031, 0.015) $\omega_{pe}^{-1}$ for X2, X1, and Z in Case M, and
0.0053 $\omega_{pe}^{-1}$ for X1 in Case S, respectively. Energy profiles are
calculated within squares plotted in Figure 4b. The video begins at
$t=0\,\omega_{pe}^{-1}$ and advances $100\,\omega_{pe}^{-1}$ at a time up to
$t=2000\,\omega_{pe}^{-1}$. The realtime duration of the video is 4s.
(An animation of the upper panel of this figure is available). Figure 4:
Distribution of maximum intensity of electric field components in the
$k_{\parallel}$-$k_{\perp}$ space as shown by the color-map of 20 $log_{10}$
$[(Ex,Ey,Ez)/(cB_{0})]$ for (a) Case W, (b) Case M, and (c) Case S, with left
panels for $E_{x}$, middle panels for $E_{y}$, and right panels for $E_{y}$.
The time interval of analysis is 1500-2000 $\omega_{pe}^{-1}$, “X2”, “X1”, and
“Z” stand for the harmonic X, fundamental X, and Z mode, respectively. Squares
in panel b represent the areas used to calculate energy of respective wave
modes. See Figure 3 for the obtained energy profiles. Figure 5: Wave
dispersion diagrams for (a) Case W, (b) Case M, and (c) Case S at three
propagation angles ($\theta_{B}=35^{\circ},95^{\circ}$, and $100^{\circ}$).
The time interval of analysis is taken to be 1500 $<\omega_{pe}^{-1}<$ 2000,
“X2”, “X1”, and “Z” stand for the harmonic X, fundamental X, and Z mode,
respectively. The video begins at $\theta_{kB}=0^{\circ}$ and advances
$5^{\circ}$ at a time up to $\theta_{kB}=180^{\circ}$. The realtime duration
of the video is 7s.
(An animation of this figure is available).


$G$ is a pseudoforest, which contradicts our hypothesis. Hence, there exists a
vertex $r\in V(G)$ with $\bar{x}_{r}>0$.

Let $T_{1}$, …, $T_{k}$ denote the components of $H(\bar{y})$. If some $T_{i}$
which is not a dyad has $\mathrm{defect}(T_{i})=0$, then $T_{i}$ is
interesting. Hence, we may assume that $\mathrm{defect}(T_{i})\geq 1$ for all
components that are not dyads. By Lemma 19, this implies
$\mathrm{defect}(T_{i})=0$ for all dyads and $\mathrm{defect}(T_{i})=1$ for
all the other components. Then, the unique component containing $r$ is
interesting. ∎

We now restate and prove Theorem 5. See 5

###### Proof.

Aiming toward a contradiction, suppose that $\bar{x}_{v}<1/3$ for all vertices
$v\in V(G)$. Hence, each one of the above lemmas apply. Let $T$ denote an
interesting component of $H(\bar{y})$, which exists by Lemma 20. Since it is
interesting, $T$ has at least two vertices and each leaf of $T$ is an edge of
$G$.

Let $r\in V(T)\cap V(G)$ denote any vertex such that $\bar{x}_{v}=0$ for all
vertices $v\in V(T)\cap V(G)$ distinct from $r$. We call $r$ the _root_ of
$T$. For $w\in V(G)$, let $d_{T}(w)$ denote the number of edges $e\in V(T)\cap
E(G)$ that are incident to $w$. (This is a slight abuse of notation since
$d_{T}(w)$ is also defined for vertices $w\in V(G)$ outside of $T$.) We let
$N(T)$ denote the set of vertices $w\in V(G)$ such that $d_{T}(w)>0$ and
$w\notin V(T)$. Finally, we let $\mathcal{L}(T)\subseteq E(G)$ denote the leaf
set of $T$.

By Lemma 17 and since $T$ is an interesting component of the support graph
$H(\bar{y})$, we have

$\displaystyle|V(T)\cap E(G)|-|V(T)\cap V(G)|$ $\displaystyle=\sum_{e=vw\in
V(T)\cap
E(G)}\left(\bar{x}_{v}+\bar{x}_{w}+\bar{y}_{e,v}+\bar{y}_{e,w}\right)-\sum_{v\in
V(T)\cap V(G)}\left(\bar{x}_{v}+\sum_{e\in\delta(v)}\bar{y}_{e,v}\right)$
$\displaystyle=(d_{T}(r)-1)\bar{x}_{r}+\sum_{w\in N(T)}d_{T}(w)\bar{x}_{w}\,.$

We notice that $|V(T)\cap E(G)|-|V(T)\cap V(G)|=|\mathcal{L}(T)|-1$. Hence,
from the equations above, we get

$(d_{T}(r)-1)\bar{x}_{r}+\sum_{w\in
N(T)}d_{T}(w)\bar{x}_{w}=|\mathcal{L}(T)|-1\,.$ (13)

Let $M\mathrel{\mathop{\mathchar 58\relax}}=(d_{T}(r)-1)+\sum_{w\in
N(T)}d_{T}(w)=(d_{T}(r)-1)+|\mathcal{L}(T)|$ denote the sum of the
coefficients in the left-hand side of (13). We note that
$d_{T}(r)\leq|\mathcal{L}(T)|$. This implies $M\leq 2|\mathcal{L}(T)|-1$.
Thus,

$\frac{1}{M}\left((d_{T}(r)-1)\bar{x}_{r}+\sum_{w\in
N(T)}d_{T}(w)\bar{x}_{w}\right)\geq\frac{|\mathcal{L}(T)|-1}{2|\mathcal{L}(T)|-1}\geq\frac{1}{3}\,.$

This is a contradiction since the left-hand side is an average of values that
are all strictly less than $1/3$. ∎

##### Acknowledgements.

Samuel Fiorini and Stefan Weltge would like to thank Gwenaël Joret and Yelena
Yuditsky for preliminary discussions. Karthekeyan Chandrasekaran and Samuel
Fiorini would like to acknowledge the support of the trimester program on
Discrete Optimization at the Hausdorff Institute for Mathematics (2021) for
facilitating initial discussions on the problem. Chandra Chekuri, Samuel
Fiorini and Stefan Weltge would like to thank the Oberwolfach Research
Institute for Mathematics and the organizers of the Combinatorial Optimization
workshop (2021) during which this research started.

## References

  * [BBF95] Vineet Bafna, Piotr Berman, and Toshihiro Fujito.  Constant ratio approximations of the weighted feedback vertex set problem for undirected graphs.  In Algorithms and Computations, pages 142–151, 1995.
  * [BG96] Ann Becker and Dan Geiger.  Optimization of Pearl’s method of conditioning and greedy-like approximation algorithms for the vertex feedback set problem.  Articial Intelligence, 83:167–188, 1996.
  * [BWB18] Austin Buchanan, Yiming Wang, and Sergiy Butenko.  Algorithms for node-weighted steiner tree and maximum-weight connected subgraph.  Networks, 72(2):238–248, 2018.
  * [BYGNR98] Reuven Bar-Yehuda, Dan Geiger, Joseph Naor, and Ron M Roth.  Approximation algorithms for the feedback vertex set problem with applications to constraint satisfaction and bayesian inference.  SIAM Journal on Computing, 27(4):942–959, 1998.
  * [CGHW98] Fabián A. Chudak, Michel X. Goemans, Dorit S. Hochbaum, and David P. Williamson.  A primal–dual interpretation of two 2-approximation algorithms for the feedback vertex set problem in undirected graphs.  Operations Research Letters, 22(4):111–118, 1998.
  * [Cha00] Moses Charikar.  Greedy Approximation Algorithms for Finding Dense Components in a Graph.  In Approximation Algorithms for Combinatorial Optimization, pages 84–95, 2000.
  * [CM16] Chandra Chekuri and Vivek Madan.  Constant factor approximation for subset feedback set problems via a new LP relaxation.  In Proceedings of the Twenty-Seventh Annual ACM-SIAM Symposium on Discrete Algorithms, SODA, pages 808–820, 2016.
  * [DS05] Irit Dinur and Samuel Safra.  On the hardness of approximating minimum vertex cover.  Annals of mathematics, pages 439–485, 2005.
  * [ENSZ00] Guy Even, Joseph Naor, Baruch Schieber, and Leonid Zosin.  Approximating minimum subset feedback sets in undirected graphs with applications.  SIAM Journal on Discrete Mathematics, 13(2):255–267, 2000.
  * [ENZ00] Guy Even, Joseph Naor, and Leonid Zosin.  An 8-approximation algorithm for the subset feedback vertex set problem.  SIAM Journal on Computing, 30(4):1231–1252, 2000.
  * [EP62] Paul Erdös and Lajos Pósa.  On the maximal number of disjoint circuits of a graph.  Publ. Math. Debrecen, 9:3–12, 1962.
  * [Fuj05] Satoru Fujishige.  Submodular functions and optimization.  Elsevier, 2005.
  * [KR08] Subhash Khot and Oded Regev.  Vertex cover might be hard to approximate to within 2- $\varepsilon$.  Journal of Computer and System Sciences, 74(3):335–349, 2008.
  * [LFFW19] Mugang Lin, Qilong Feng, Bin Fu, and Jianxin Wang.  An approximation algorithm for the $\ell$-pseudoforest deletion problem.  Theoretical Computer Science, 806, 08 2019.
  * [MS81] Burkhard Monien and Reinald Schulz.  Four approximation algorithms for the feedback vertex set problem.  In WG, pages 315–326, 1981.
  * [PRS18] Geevarghese Philip, Ashutosh Rai, and Saket Saurabh.  Generalized pseudoforest deletion: Algorithms and uniform kernel.  SIAM Journal on Discrete Mathematics, 32(2):882–901, 2018.
  * [Wes00] Douglas B. West.  Introduction to Graph Theory.  Prentice Hall, 2000.