# Algorithms for Computing the Free Distance of Convolutional Codes

Zita Abreu, Joachim Rosenthal, Michael Schaller

###### Abstract

The free distance of a convolutional code is a reliable indicator of its
performance. However its computation is not an easy task. In this paper, we
present some algorithms to compute the free distance with good efficiency that
work for convolutional codes of all rates and over any field. Furthermore we
discuss why an algorithm which is claimed to be very efficient is incorrect.

##  1 Introduction

Since real channels are inherently noisy, error correcting codes are necessary
in all communication systems that handle digitally represented data.
Convolutional codes are one kind of error-correcting codes that is well-known
for its suitability for sequential decoding with low delay. This class of
codes has been extensively researched, see [1, 2]. The distance of a
convolutional code indicates the code’s robustness since it allows us to
analyze its capacity to protect data against errors. Codes with greater
distance are preferable because they allow us to rectify more errors. The free
distance is one of the main types of distances for convolutional codes. This
kind of distance is important for the decoding of a complete message.

In general, calculating the free distance of convolutional codes is not an
easy task. Up to now, there are not many known algorithms that calculate the
free distance. Actually, not much has been done lately to develop new
algorithms or even improve old ones. Most algorithms in the literature are
either too slow or use too much computation and storage or are even wrong, see
e.g. [3], [4], [5], [6] and [7]. Additionally, they are designed only over
$\mathbb{F}_{2}$ and for certain rates. For instance, the authors of [6] and
[7] made reference to the possibility of a generalization of the proposed
algorithms, but they did not elaborate.

In this article we present some algorithms for computing the free distance of
convolutional codes. In Section 2 some introductory concepts are given.
Section 3 provides a naive algorithm for calculating the free distance of
convolutional codes. In section 4 we prove that one of the existing algorithms
for calculating the free distance published in [3] is wrong. In section 5, we
improve and generalize the algorithm provided in [7], making it work for
convolutional codes of all rates and over any field. Finally, in Section 6, we
present a novel algorithm for determining the free distance that also supports
convolutional codes with all rates and degrees over any field. This algorithm
combines the algorithm from Section 5 with the algorithm from article [5].

##  2 Preliminaries

In this section, we provide a few key definitions and results for the
following sections. For more details, see [8] or [2].

First, denote by $\mathbb{F}_{q}[z]$ the ring of polynomials over the finite
field with $q$ elements $\mathbb{F}_{q}$.

######  Definition 1.

A convolutional code $\mathcal{C}$ of rate $k/n$ is a
$\mathbb{F}_{q}[z]$-submodule of $\mathbb{F}_{q}[z]^{n}$ of rank $k$. A matrix
$G(z)\in\mathbb{F}_{q}[z]^{k\times n}$ whose rows constitute a basis of
$\mathcal{C}$ is known as a generator matrix for $\mathcal{C}$, i.e.:

$\displaystyle\mathcal{C}$
$\displaystyle=\\{{v(z)\in\mathbb{F}_{q}[z]^{n}:v(z)=u(z)G(z)\text{ with
}u(z)\in\mathbb{F}_{q}[z]^{k}\\}.}$

######  Definition 2.

Let $G(z)=\sum_{i=0}^{\mu}G_{i}z^{i}\in\mathbb{F}_{q}[z]^{k\times n}$ with
$G_{\mu}\neq 0$ and $k\leq n$. For each $i$, $1\leq i\leq k$, the $i$-th row
degree $\nu_{i}$ of $G(z)$ is defined as the largest degree of any entry in
row $i$ of $G(z)$, in particular $\mu=\max_{i=1,\ldots,k}\nu_{i}$. The
external degree of $G(z)$ is the sum of the row degrees of $G(z)$. The
internal degree of $G(z)$ is the maximal degree of the $k\times k$ minors of
$G(z)$ and it is the same for every generator matrix of the same code.

######  Definition 3.

A matrix $G(z)\in\mathbb{F}_{q}[z]^{k\times n}$ is said to be row reduced if
its internal and external degrees are equal. In this case, $G(z)$ is called a
minimal generator matrix of the convolutional code it generates. The degree
$\delta$ of a code $\mathcal{C}$ is the internal degree of a generator matrix
of $\mathcal{C}$. A convolutional code with rate $k/n$ and degree $\delta$ is
known as a $(n,k,\delta)$ convolutional code.

######  Definition 4.

The free distance of a convolutional code $\mathcal{C}$ is given by

$d_{free}(\mathcal{C}):=\min_{v(z)\in\mathcal{C}}\left\\{wt(v(z))\ |\ v(z)\neq
0\right\\},$

where $wt(v(z))$ is the Hamming weight of
$v(z)=\sum_{t=0}^{\deg(v(z))}v_{t}z^{t}\in\mathbb{F}_{q}[z]^{n}$ that is
defined as $wt(v(z))=\sum^{\deg(v(z))}_{t=0}wt(v_{t}),$ where the weight
$wt(v)$ of $v\in\mathbb{F}_{q}^{n}$ is the number of nonzero components of
$v$.

######  Definition 5.

For $j\in\mathbb{N}_{0}$, the j-th column distance of a convolutional code
$\mathcal{C}$ is defined as

$d_{j}^{c}(\mathcal{C}):=\min\left\\{wt(v_{0},\ldots,v_{j})\ |\
{v}(z)\in\mathcal{C}\text{ and }{v}_{0}\neq 0\right\\}.$

######  Definition 6.

Let $\mathcal{C}$ be an $(n,k)$ convolutional code over $\mathbb{F}_{q}$. A
full row rank matrix $H(z)\in\mathbb{F}_{q}[z]^{(n-k)\times n}$ satisfying

$\mathcal{C}=\ker
H(z)=\\{v(z)\in\mathbb{F}_{q}[z]^{n}\,:\,H(z)v(z)^{\top}=0\\}$

is called a parity-check matrix of $\mathcal{C}$. If this matrix exists,
$\mathcal{C}$ is named non-catastrophic, otherwise it is named catastrophic.

Moreover, from [2] we know that a code is non-catastrophic if and only if
$G(z)$ is left prime which is equivalent to $G(z)$ having finite output
$u(z)G(z)\in\mathbb{F}_{q}[z]^{n}$ with $u(z)\in\mathbb{F}_{q}((z))^{k}$
implies that the input was already finite, i.e.,
$u(z)\in\mathbb{F}_{q}[z]^{k}$.

######  Definition 7.

Let $\mathcal{C}$ be a non-catastrophic convolutional code generated by the
row reduced matrix $G(z)$ with entries $g_{ij}(z)$. Let
$\overline{g}_{ij}(z)=z^{\nu_{i}}g_{ij}(z^{-1})$ where $\nu_{i}$ is the $i$-th
row degree. Then the code $\overline{\mathcal{C}}$ with generator matrix
$\overline{G}=(\overline{g}_{ij})$ is called the reverse code of
$\mathcal{C}$.

In practice, instead of moving backward via the state-transition diagram of
$\mathcal{C}$, one might go forward through the state-transition diagram of
$\overline{\mathcal{C}}$.

The next theorem provides an upper bound for the free distance.

######  Theorem 1 ( [9]).

Let $\mathcal{C}$ be an $(n,k,\delta)$ convolutional code. Then,
$d_{free}(\mathcal{C})\leq(n-k)\left(\left\lfloor\frac{\delta}{k}\right\rfloor+1\right)+\delta+1.$
This bound is called the generalized Singleton bound.

##  3 Naive computation of the Free Distance

This section provides a very simple method to calculate the free distance of a
non-catastrophic convolutional code. The main idea behind this method is to
find a non-trivial shortest path through the state-transition-diagram from the
zero-state to the zero-state, with shortness assessed in terms of the weight
of the output.

Consider $W^{\ast}$ to be an upper bound on the free distance. As example,
assign $W^{\ast}$ to the generalized Singleton bound, see Theorem 1.

Essentially, we will want to explore all paths until we either return to the
zero-state and adjust the upper bound, or until the weight of the path exceeds
$W^{\ast}$.

To ensure consistency across the algorithms provided here, we will describe
this method by subtracting weights from the upper bound and from the states
and checking if the residual weight is non-negative. Furthermore we move
backwards through the state-transition-diagram.

Let us now describe the method as follows: We initialize an empty stack. We
start at the zero-state and consider all possible extensions $S_{E}$ except
for the ones that corresponds to the all zero input. If there are extensions
$S_{E}$ from the zero-state to the zero-state which do not correspond to the
all zero input (which can only happen if a row degree is $0$) we calculate the
weights of the outputs corresponding to these extensions and set $W^{\ast}$ to
the minimum of $W^{\ast}$ and these weights. We do not put these extensions on
the stack. For all other extensions we set the weight $W_{E}$ of $S_{E}$ to be
$W^{\ast}-w_{E}$ where $w_{E}$ is the weight of the output corresponding to
the extension. Then we store the extensions $S_{E}$ together with the weights
$W_{E}$ in the stack if the weights are non-negative. Next, we select a pair
$(S,W)$ from the stack and look at all the extensions $(S_{E},W_{E})$. The
weight $W_{E}$ of an extension $S_{E}$ is now $W-w_{E}$. If the extension is
the zero-state and if the weight $W_{E}>0$, we change the upper bound
$W^{\ast}$ to $W^{\ast}-W_{E}$ and we adjust the weight of the current state
to $W-W_{E}$. Moreover we adjust the stack, i.e., delete the states with
weights less than $W_{E}$ and subtract $W_{E}$ from the weights of the
remaining ones. If the extension is not the zero-state and if $W_{E}\geq 0$ we
save $(S_{E},W_{E})$ on the stack. We continue in the same way until the stack
is empty. Finally, we return $W^{\ast}$ which is the free distance.

To demonstrate the method, consider the following example in the Figure 1.

$\begin{aligned} 00\\\ 6\hskip 2.84544pt\end{aligned}$$\begin{aligned} 01\\\
4\hskip 2.84544pt\end{aligned}$11$\begin{aligned} 10\\\
\cancel{4}{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1}}\end{aligned}$00$\begin{aligned}
\hskip 0.28453pt11\\\
\text{{\cancel{2}{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}-1}}}}\end{aligned}$11$\begin{aligned}
00\\\
\cancel{3}{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0}}\end{aligned}$10$\begin{aligned}
01\\\ \text{ 0 }\end{aligned}$01$\begin{aligned} 10\hskip 2.84544pt\\\ \text{0
\hskip 1.42271pt}\end{aligned}$00$\begin{aligned} 11\hskip 2.84544pt\\\ \text{
{-2} \hskip 2.84544pt}\end{aligned}$11$\begin{aligned} 00\hskip 2.84544pt\\\
\text{ {-1} \hskip 2.84544pt}\end{aligned}$10$\begin{aligned} 01\hskip
2.84544pt\\\ \text{{-1} \hskip 1.42271pt}\end{aligned}$01 Figure 1: Naive
computation of the free distance for a $(1,2,2)$ convolutional code with
$G(z)=\begin{bmatrix}1+z^{2}&z^{2}\end{bmatrix}\in\mathbb{F}_{2}[z]^{1\times
2}$.

######  Example 1.

$W^{\ast}$ is initialized as $6$. While we move backwards in the state-
transition-diagram, the only possible extension is the extension by $1$, i.e
$S_{1}=01$. In this state, $W_{1}=W^{\ast}-w_{1}=6-wt(11)=6-2=4$. This is not
the zero-state and $W_{1}>0$, so we add $(01,4)$ to the stack. Since it is the
only pair in the stack, we select it and look at all its extensions. Extension
$10$ will have weight $W-W_{0}=4-wt(00)=4-0=4$ and extension $11$ will
similarly have weight $2$. Again none of these states are the zero-state and
$W_{E}>0$ for both, so $(10,4)$ and $(11,2)$ are added to the stack and
$(10,4)$ will be on top of the stack. Popping from the stack we get $S=10$.
Since, one of the extensions of $S$ is $S_{0}=00$ with $W_{0}=3$, we update
the upper bound to $W^{\ast}=W^{\ast}-W_{0}=6-3=3$. Additionally, we update
$W$ to $W-W_{0}=4-3=1$. Adjusting the stack, the state $11$ with weight $2$ is
deleted from the stack. We continue with the extension of $S=10$ by $1$ to get
$S_{1}=01$ with weight $0$. Then we put it on the stack and pick it
immediately again. So $S=01$ and $W=0$. We compute the extensions $S_{0}=10$
with weight $W_{0}=0$ and $S_{1}=11$ with $W_{1}=-2$. Only $(S_{0},W_{0})$
gets put on the stack. Then we pick it from the stack. Now both extensions
have negative weight and the stack is empty, so we stop and return
$W^{\ast}=3$ as the free distance.

Observe that, the algorithm will terminate since it will encounter repeated
states for each path it can take in the tree. As the code is non-catastrophic,
the weight of the repeated state $S$ will be less than it was previously when
encountered again. Otherwise there would be a cycle of only zero output from
$S$ to $S$, resulting in infinite input and finite output, which, by
definition, does not occur if the code is non-catastrophic. Obviously this
argument provides only a rough upper bound for the number of steps that the
algorithm performs. In practice, it will terminate with considerably less
steps than the bound.

Overall, this naive computation of the free distance works over
$\mathbb{F}_{q}$ and for convolutional codes of all rates, but it performs
poorly. Therefore, we focused on developing better algorithms to calculate the
free distance.

##  4 The Heapmod algorithm is wrong

In this section we study the algorithm presented in [3], called Heapmod
algorithm and claimed to be better than the one published in [7], and we
demonstrate that it is incorrect. To be more specific, we rewrite the stop
condition of the algorithm to provide a clearer understanding of its
incorrectness. With our version of the stop condition, it is evident that the
incorrectness is related to an old conjecture by Costello [10] and a class of
counterexamples presented in [11] by Miczo and Rudolph.

The authors of the Heapmod algorithm construct the code tree in a slightly
different manner than usual. Namely, instead of having the state on the nodes
and the outputs on the edges they have the state together with the input on
the node and there are no outputs on the edges. A decimal representation of
registers is used, where the registers consist of input and state. To clarify,
see figure 2 to check how the first part of the tree looks like. Note that the
entire tree is not drawn.

$\scriptstyle{1}$
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\scriptstyle{(10,1\mod
4})}$ $\scriptstyle{101}$ $\scriptstyle{2}$
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\scriptstyle{(01,2\mod
4})}$ $\scriptstyle{011}$ $\scriptstyle{3}$
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\scriptstyle{(11,3\mod
4})}$ $\scriptstyle{110}$ $\scriptstyle{4}$
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\scriptstyle{(00,0\mod
4})}$ $\scriptstyle{000}$ $\scriptstyle{5}$
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\scriptstyle{(10,1\mod
4})}$ $\scriptstyle{101}$ $\scriptstyle{6}$
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\scriptstyle{(01,2\mod
4})}$ $\scriptstyle{011}$ $\scriptstyle{7}$
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\scriptstyle{(11,3\mod
4})}$ $\scriptstyle{110}$ • Index $i$ • Register • Output Figure 2:
Illustration of the Heapmod algorithm.

In this section we let $\mathcal{C}$ be a $(n,1,M)$ convolutional code. We
state the algorithm below, using the same notation as in [3].

  1. 1.

Calculate the weight $w(O_{i})$ of the output $O_{i}$ corresponding to
register $i$ for $i\in\\{1,\ldots,2^{M+1}-1\\}$.

  2. 2.

Compute

$l_{j}=(2j+1)\cdot 2^{M}$

for $j\in\\{0,\ldots,2^{M+1}-1\\}$.

  3. 3.

For each $j$ form a path from $l_{j}$ to the root of the tree and sum the
weights of the corresponding outputs. The minimum over $j$ of these sums is
the free distance.

Taking into account the tree with $M=1$ whose first part is represented in
Figure 2, the indices $l_{j}$ would be $2,6,10,14$. Note that coming to
$l_{j}$ corresponds to a return to the zero-state. This is easy to see, since
$l_{0}=2^{M}$ and $l_{j+1}=(2(j+1)+1)\cdot 2^{M}=(2j+1+2)\cdot
2^{M}=l_{j}+2^{M+1}$. In this way, adding $2^{M+1}$ gives us exactly the next
occurrence of the same register. Therefore, let us reformulate the algorithm
in the following way: It considers the $2^{M+1}$ shortest (in terms of steps
through the state-transition-diagram) paths from the zero-state to the zero-
state.

A simple counterexample for the correctness of the Heapmod algorithm is the
non-catastrophic code generated by

$G(z)=\begin{pmatrix}z^{6}+z^{4}+1,&z^{6}+z^{5}+z^{4}+z^{3}+z+1\end{pmatrix}.$

For this example, the Heapmod algorithm returns a free distance of $9$ instead
of $8$ (which was the correct answer). To be explicit
$wt((1+z^{4}+z^{6}+z^{8})G(z))=wt(\begin{pmatrix}z^{14}+1,&z^{14}+z^{13}+z^{9}+z^{3}+z+1\end{pmatrix})=8$.
Similar examples can easily be obtained.

The underlying reason is the following: For each $j$ the length from $l_{j}$
to the root of the tree is at most

$\log_{2}(l_{j})\leq\log_{2}(2\cdot 2^{M+1}\cdot 2^{M})=2M+2.$

That means that the algorithm explores only inputs up to length $2M+2-M=M+2$
since we also need $M$ zero inputs at the end to arrive at the registers
corresponding to the $l_{j}$. Note that in the example above the input length
is $9$ since the degree of the input polynomial is $8$. Furthermore, since the
input length is 9, it is evident that the algorithm does not explore
everything because $M=6$ and the algorithm only explores inputs up to length
$M+2=8$.

It was a conjecture of Costello in [10] that inputs of length $2M$ would be
enough to determine the free distance. In the same paper he also proved an
upper bound quadratic in $M$ for the input length for systematic codes. Miczo
and Rudolph proved in [11] that no linear upper bound on the input length is
enough to get the free distance, which also gives another proof of the
incorrectness of the Heapmod algorithm.

##  5 The optimized FAST algorithm

In this section we will introduce an optimized version of the FAST algorithm,
a stack algorithm which computes the free distance of convolutional codes,
presented in [7]. More precisely, we present a generalized version that works
for arbitrary convolutional codes of all rates and over any field, we explain
the modifications we made and we address their impact in relation to the
efficiency of the algorithm.

In the Fast algorithm, the last non-zero entry in a state $S$ tells us the
minimum number of steps $u$ we must go in the backward direction to reach the
zero-state. Additionally, moving forward from the zero-state to $S$ requires
at least $u$ inputs, which tells us that the weight of this path is at least
the $u-1$-st column distance $d_{u-1}$. This means that if the weight of a
state $S$ is less than $d_{u-1}$ we do not have to extend this state, see
Figure 3.

$00\cdots 0$ S $*0\cdots 01$ $S$ $*1\underbrace{0\cdots 0}_{m\mathrm{\
zeros}}$ $S_{m}$ $*1\underbrace{0\cdots 0}_{m-1\mathrm{\ zeros}}$ $\geq M$
steps$\geq M-m$ steps$E\neq 0$$E=0$ Figure 3: Idea behind the (optimized) Fast
algorithm for a convolutional code of rate $1/n$ over $\mathbb{F}_{2}$.

The optimized Fast algorithm aims to avoid unnecessary state extensions in the
code tree that are not considered in the Fast algorithm. If one calculates the
column distances $[d_{0},\ldots,d_{M}]$ using brute force one can improve this
by storing the actual distance of the state to the zero-state. Here distance
refers to the path with the fewest steps, which is not always the path with
lowest weight. This costs more memory, but can make the calculation
considerably faster. See Figure 4 and 5, which illustrates the difference in
efficiency between the Fast algorithm and the optimized Fast algorithm for
eighty arbitrarily chosen convolutional codes from the literature with
different rates and degrees and over distinct fields. Figure 4 correlates the
number of new nodes visited to the free distance of the chosen convolutional
codes, whereas Figure 5 shows the number of nodes visited for each
convolutional code separately. Both show a significant reduction in the number
of nodes we need to visit using the optimized Fast algorithm.

Figure 4: Efficiencies of the FAST algorithm and the optimized FAST algorithm.
Figure 5: Another perspective of Figure 4.

Let us first introduce the notation for presenting the optimized Fast
algorithm (Algorithm 1). The integer $m(S)_{i}$ represents the number of zeros
to the right of the rightmost non-zero entry in row $i$ of state $S$. The
$d(S)$ is the distance of a state to zero where we only consider paths of
length at most $M+1$, where $M$ is the memory. This distance will be
precomputed when we compute the column distances. Finally, the $j$-th column
distance is denoted by $d_{j}$.

Algorithm 1 OPTIMIZED FAST ALGORITHM

Given a generator matrix $G$ of a convolutional code $\mathcal{C}$ with row
degrees $\nu_{i}$ and $W^{*}$ an upper bound of the free distance, we search
the code tree of $\mathcal{C}$ to determine the free distance.

Calculate $d_{j}$ with $j=0,\ldots,M$ and the distance $d(S)$ of the state $S$
to the zero-state.

Calculate all extensions $S$ of the zero-state corresponding to non-zero input
together with their weights $W$. If $S$ ends in the zero-state set
$W^{\ast}=\min(W^{\ast},W)$ and do not put the extension on the stack. For the
states $S$ which do not end in the zero-state put the tuples $(S,W)$ on the
stack if $W\geq 0$.

while stack is nonempty do

Take $(S,W)$ from the stack.

Compute all extensions $S_{E}$ of $S$ and their weights $W_{E}$.

Go through all of them.

if $S_{E}$ is the zero-state and $W_{E}\geq 0$ then

Set $W^{*}=W^{*}-W_{E}$, adjust the weights of the extensions and of the
states on the stack and continue.
---

else

$\sigma=\max(\nu_{i}-m(S_{E})_{i})$

if (($W_{E}<d(S_{E})$ and $W_{E}<d_{M}$) or $W_{E}<d_{\sigma-1}$) then

Discard this extension.

end if

end if

Pick one of the extensions that are left and set it to $S$ and its weight to
$W$.
---

Put the other extensions that are left on the stack.
---

end while

Return $W^{\ast}$

For a better understanding of the algorithm, consider that the variable
$\sigma$ represents the minimum number of steps required to return to the
zero-state. Note that this variable is used in the second part of the second
if in order to discard the extension.

Regarding the first part of the second if, the condition $W_{E}<d(S_{E})$
implies that we cannot return to the zero-state from $S_{E}$ in less than
$M+1$ steps. To get to the zero-state, we need to take at least $M+1$ steps.
If $W_{E}<d_{M}$, we also cannot return to the zero-state in $M+1$ or more
steps.

Let us conclude with some observations on the algorithm’s implementation. We
recommend to pick the extension by zero if it is one of the remaining
extensions. Returning to the zero-state can lead to an adjustment of the upper
bound, which can significantly reduce the runtime. It may even make sense to
employ more advanced heuristics to swiftly return to the zero-state and,
hopefully, adjust the upper bound.

##  6 A new algorithm

In this section, we propose a new algorithm for computing the free distance of
convolutional codes that is based on a combination of two algorithms, namely
the optimized FAST algorithm from Section 5 and the algorithm proposed by
Larsen in [5]. The last one is a bidirectional algorithm that extends paths
forward and backward simultaneously. Our algorithm is also bidirectional and
it works for arbitrary convolutional codes of all rates and over any field.

This new algorithm also uses a state-transition diagram and keeps the notation
used of Algorithm 1. However, now when a path reaches a state, we stored it in
an array together with the information about the path’s type (forward or
backward) and the weight. If there are many paths to a state, the lowest
weight is stored.

In this algorithm we will also need to calculate the weight by starting from
the zero-state with weight $0$ and increasing it. We denote this weight by
$W_{L}$. Note that, for the weight $W$ which we get by starting from the zero-
state with weight $W^{\ast}$ and decreasing it we get $W+W_{L}=W^{\ast}$.

The following new notation must be fixed

  * •

$T$ \- Type of state $S$. The type consists of two parts, $T1$ and $T2$. $T1$
indicates whether the state $S$ is $D$(ead) or non$-D$ and $T2$ is the type,
$F$(orward) or $B$(ackward), of the minimum-weight path first found to $S$.
Thus we have four possible types: $F$, $B$, $DF$, and $DB$.

Note that if a variable depends on $T$ we use the same notation as in
Algorithm 1 in case $T2=B$. If $T2=F$, we use the same variables as in the
fast algorithm but for the reverse code generated by the generator matrix
$\overline{G}$. By $G(T)$ we denote $G$ if $T=B$ and $\overline{G}$ if $T=F$.
Similarly $\nu_{i}(T)$ is the $i$-th row degree of $G$ if $T=B$ and the $i$-th
row degree of $\overline{G}$ if $T=F$.

The variable $W^{\ast}$ remains the current upper bound of the free distance
and the new algorithm for a rate$-k/n$ convolutional code is stated as
Algorithm 2. In the algorithm where there is continue we skip all further
steps in the for-loop and proceed in the for-loop with the next extension $E$.

Algorithm 2 NEW ALGORITHM

Given a generator matrix $G$ of a convolutional code $\mathcal{C}$ with row
degrees $\nu_{i}$ and $W^{*}$ an upper bound of the free distance, we search
the code tree of $\mathcal{C}$ to determine the free distance.

1) Start at the zero-state and calculate all forward and backward extensions
$S$ and their weights $W_{L}$. If some of these extensions have non-zero input
and give the zero-state adjust the upper bound $W^{\ast}$. We do not save
these extensions. We store all other extensions with the following exception:
If a $B$-extension is equal to some $F$-extension save only the one with lower
weight $W_{L}$ and adjust the upper bound if the sum of the weights is less
than $W^{\ast}$. Together with a state $S$ we also store its weight $W_{L}$
and type $T$.

while not Done do

2) Search through the storage array for the lowest weight non$-D$ state $S$.
It has weight $W_{L}$ and is of type $T$.
---

for extensions $E$ do

3) Determine terminal state $S_{E}$ and weight $W_{E,L}$ of the $E$-extension
of $S$.
---

4)  If $S_{E}$ is the zero-state set $W^{\ast}=\min(W_{E,L},W^{\ast})$ and
continue.
---

5) $\sigma(T)=\max(\nu_{i}(T)-m(S_{E},T)_{i})$. If
$W_{E,L}>\frac{W^{\ast}+n-1}{2}$ or ($W_{E}<d(S_{E},T)$ and $W_{E}<d_{M}(T)$)
or $W_{E,L}<d_{\sigma-1}(T)$ continue.
---

6) Check through the array for $S_{k}=S_{E}$. If no such $S_{k}$ is found,
store $S_{E},W_{E,L},T$ in the array and continue.
---

7) The type and weight of $S_{k}$ are $T_{k}$ and $W_{k,L}$. If ${T2}_{k}=T2$,
go to 11).
---

8) Set $W^{\ast}=\min(W^{\ast},W_{E,L}+W_{k,L})$.

9) If $W_{E,L}\geq W_{k,L}$, continue.

10) Set $W_{k,L}=W_{E,L}$, ${T2}_{k}=T2$ and continue.

11) If $T1_{k}=D$, continue.

12) Set $W_{k,L}=\min(W_{k,L},W_{E,L})$.

end for

13) Set $T1=D$.

14) Set $\textit{Done}=(2W_{L}\geq W^{\ast}$ or all states are $D)$.

end while

return $d_{free}=W^{\ast}$.

Regarding paper [5], it is interesting to note that a proof of correctness of
the algorithm proposed there is incomplete. In fact, we demonstrate this in
the Appendix and we explain why this can be remedied.

The main improvements of the new algorithm compared to the Algorithm in [5]
are that this algorithm works for all rates, degrees and over any field and
that the algorithm visits fewer states.

In short, regarding the algorithm, step 4 is required in order to generalize
it. More precisely, without this step the algorithm is wrong if $k>1$.
Additionally, with this step we also circumvent the issues in the proof of the
algorithm proposed by Larsen, which we elaborate on in the Appendix. Step 5,
on the other hand, allows us similarly as in the optimized Fast algorithm to
skip certain states, namely avoiding unnecessary state extensions in the code
tree.

Finally, Figure 6 shows how the new algorithm compares to the Optimized Fast
Algorithm for the same eighty convolutional codes as before. Note that this
comparison is not entirely fair since the new algorithm has to search through
an array to visit the nodes.

Figure 6: Efficiencies of the Optimized Fast algorithm and the new algorithm.

##  7 Conclusion and Future Work

Some algorithms to compute the free distance of convolutional codes were
presented. Depending on the convolutional code the bottleneck of these
algorithms might either be the number of states visited or the precomputation
of the column distances. Therefore an efficient method for computing the
$j$-th column distance is an interesting topic for further research.

Additionally, for constructing convolutional codes with good free distance it
will be important to have algebraic criteria. The state-transition diagram can
be represented algebraically by linear system representations, see [12]. It
would be therefore interesting to review the new algorithm presented here in
the algebraic systems representation.

To see the implementations of the algorithms mentioned here, visit
https://github.com/uscpr/algorithms-for-computing-the-free-distance-of-
convolutional-codes.

##  Appendix A Comments on Larsen’s Algorithm

In this section, we will review the algorithm presented by Larsen in [5]. We
will demonstrate that the proof of the correctness of the mentioned algorithm
is incomplete. However, despite this, the algorithm works. Therefore, we
generalize it and use it as a reference for the construction of the new
Algorithm presented in Section VI.

In general, the proof the correctness of the algorithm is incomplete because
it only implicitly takes the structure of the convolutional code into account.
In fact, it is easy to construct a state-transition diagram that looks like
the state-transition diagram of a convolutional code but on which the
algorithm fails if we are allowed to choose the weights on the edges as we
wish. Consider the following example.

######  Example 2.

Take the following state-transition diagram with the respective weights on the
edges which has the same graph as a state-transition diagram for a
convolutional code of memory $2$, see Figure 7.

$01$$11$$00$$10$20102001201 Figure 7: Example of a state-transition diagram
that verifies that the proof of correctness of the algorithm proposed by
Larsen is incomplete.

Before we introduce the algorithm proposed by Larsen, let us take into account
that in this algorithm when a state is reached by a path it is stored in an
array together with information on the type of the path (forward or backward)
and the weight of the path. If there are many paths to a state, then the
lowest weight is stored. The notation used is the same as in the new algorithm
described in Section VI. For completeness we state the algorithm introduced by
Larsen, see Algorithm 3.

Algorithm 3 ALGORITHM PRESENTED BY LARSEN

Let $W^{\ast}$ denote the current upper bound on the free distance. Then the
algorithm for computing the free distance a $1/n$ convolutional code is the
following.

1) Set $W^{\ast}$ as an upper bound on the free distance. Store
$S=(1,0\ldots,0)$, together with the weight $W_{L}$ of $S$ we get from the
edge coming from $(0,\ldots,0)$ and $T=F$ and $S=(0,0\ldots,01)$, together
with its weight $W_{L}$ we get from the edge going to $(0,\ldots,0)$ and $T=B$
in the array.

2) Search through the storage array for the lowest weight non$-D$ state $S$.
It has weight $W_{L}$ and is of type $T$.

3) If $2W\geq W^{\ast}$ or all states are $D$, go to 17).

4) Set $E=0$. Determine terminal state $S_{E}$ and weight $W_{E,L}$ of the
$0$-extension of $S$. Go to 6).

5) Set $E=1$. Determine terminal state $S_{E}$ and weight $W_{E,L}$ of the
$1$-extension of $S_{m}$.

6)If $W_{E,L}>\frac{(W^{\ast}+n-1)}{2}$, go to 15).

7) Check through the array for $S_{k}=S_{E}$. If such an $S_{k}$ is found, go
to 9).

8) Find an unoccupied location in the array and store $S_{E}$,$W_{E,L}$, $T$.
Go to 15).

9) The type and weight of $S_{k}$ are $T_{k}$ and $W_{k,L}$. If
${T2}_{k}={T2}$ go to 13).

10) Set $W^{\ast}=\text{min}(W^{\ast},W_{E,L}+W_{k,L})$.

11) If $W_{E,L}\geq W_{k,L}$, go to 15).

12) Set $W_{k,L}=W_{E,L}$ and ${T2}_{k}={T2}_{m}$. Go to 15).

13) If ${T1}_{k}=D$, go to 15).

14) Set $W_{k,L}=\text{min}(W_{k,L},W_{E,L})$.

15) If $E=0$, go to 5).

16) Set ${T1}=D$. Go to 2).

17) Return $d_{free}=W^{\ast}$.

Let us consider now for Example 2 the steps of the algorithm. In step 1, we
set, for example, $W^{\ast}=100$ and we put $((10),1,F)$ and $((01),10,B)$ in
the array. In step 2, we set $S=(10),W_{L}=1,T=F$. The condition 3 is not
satisfied so we go to step 4. The $0$-extension is $S_{0}=(01)$ and
$W_{0,L}=3$. Condition 6 is not satisfied so we go to step 7. Since $(01)$ is
in the array, we go to step 9. We set $S_{k}=(01),W_{k,L}=10,T_{k}=B$. Since
$T2_{k}\neq T2$, we go to step 10 and set $W^{\ast}=13$. As $W_{0,L}<W_{k}$,
we go to 12 and fix $W_{k,L}=3$ and $T2_{k}=F$. Therefore the array now
contains $((10),1,F)$ and $((01),3,F)$. Then we consider the $1$-extension of
$S=(10)$ in the $F$ direction which is $S_{1}=(11)$. So $((11),2,F)$ is stored
in the array. After that $S=(10)$ will be set dead and we go to step 2 again.
Here pick $S=(11),W_{L}=2,T=F$. Condition 3 is not satisfied, so we go to step
4. The $0$-extension of $S$ is $S_{0}=(01)$ and $W_{0,L}=2$. Condition 6 is
not valid and we find in 7 that $S_{k}=(01)$ is in the array. In 9, we get
$W_{k,L}=3,T_{k}=F$ and go to step 13. The condition in 13 is not satisfied
and we go to 14. In step 14, we set $W_{k,L}=2$. Afterwards we continue with
the $1$-extension. Nothing changes in the array through the steps of the
$1$-extension. Afterwards we set $(11)$ dead. Then in step 2 we go look for
the lowest weight non-dead state in the array which is $S=(01)$ with $W_{L}=2$
and $T=F$. In step 3 we continue. We consider the $0$-extension of $S$ which
is $S_{0}=(00)$ with $W_{0,L}=12$. Let us assume that the condition 6 is not
satisfied (otherwise it becomes only easier). The state $(00)$ is not in the
array, so we store it with weight $12$. Again in the step for the
$1$-extension nothing changes in the array. Then $(01)$ is set dead. At this
moment only $(00)$ is the non-dead state in the array. Obviously $2\cdot
12=24>W^{\ast}=13$ and we stop and return $13$. This is not the distance of
the shortest non-trivial path from $(00)$ to $(00)$. The correct answer would
be $12$.

We will now describe the problem in the proof presented in [5] and explain why
the algorithm still works for convolutional codes.

Firstly, Larsen considers an arbitrary path in the state-transition diagram
and uses a case analysis with two possibilities. The first one is that all
non-zero states in the path are dead and the second is that there are paths
with some non-dead states. The author claims that in the first case, the case
with only dead states, there must be a part in the path where two adjacent
states are of different $T2$ type. Obviously in the example above this is not
true. In reality, it is simple to construct an example of a convolutional code
where this is not the case, such as the code generated by
$G(z)=(z^{2}+1,z^{2})$. Going through the algorithm we can observe, as
previously, that the claim is incorrect for the code generated by $G(z)$.
However now the algorithm returns the correct result, which is $3$. Let us
elaborate a bit more.

Let us assume that there are two previous states for the penultimate state
$(0\cdots 01)$, namely $(0\cdots 011)$ and $(0\cdots 010)$. The crucial point
for the counterexample is that there are $F$-paths from $(0\cdots 0)$ through
both $(0\cdots 011)$ and $(0\cdots 010)$ to $(0\cdots 01)$ which have weight
less than the weight of the edge from $(0\cdots 01)$ to $(0\cdots 0)$. More
intuitively and considering Figure 8, the weight of the green and red paths is
smaller than the weight of the blue edge.

$0\cdots 000$$10\cdots 00$$0\cdots 011$$0\cdots 010$$0\cdots 001$ Figure 8:
Illustration on the idea to fix the proof in Larsen’s paper for convolutional
codes of rate $1/n$ over $\mathbb{F}_{2}$.

However, we claim that for a convolutional code it can never happen that a
path from $(0\cdots 0)$ through $(0\cdots 011)$ to $(0\cdots 01)$ has weight
less than the edge from $(0\cdots 01)$ to $(0\cdots 0)$, i.e., we claim that
the green path has weight at least the weight of the blue edge. To prove this,
let us assume the opposite. Then at least one of the $n$ outputs (where $n$ is
the code length) must give a $1$ on the blue edge and a $0$ on the green path.
Therefore it must be some linear combination of the registers with the
following properties. We write the linear combination in polynomial form
$\sum_{i=0}^{M}a_{i}z^{i}$, where $a_{i}\in\mathbb{F}_{2}$ indicates which
registers contribute to the output. We need to have $a_{M}=1$, otherwise the
blue edge has a $0$ as output. Furthermore, $a_{M-1}=1$, otherwise the output
for the transition from $(0\cdots 011)$ to $(0\cdots 01)$ would give output
$1$ and this transition belongs to the green path. Writing it as a polynomial
we get $\sum_{i=0}^{M-2}a_{i}z^{i}+z^{M-1}+z^{M}$. But then any input will
give weight at least $2$ which means that the green path has weight at least
$1$ and so the claim follows.

In sum, we can conclude the correctness of the algorithm as follows. Either
the red path has at least the same weight as the blue edge. In this situation,
Larsen’s claim that in every path in which all non-zero states are dead there
exist two adjacent states of different $T2$ type is true and the correctness
remains as it is in the rest of Larsen’s proof. If the red path has weight
less than the blue edge, then in particular it has weight less than the green
path as we have shown above. This means that the adjustment of the upper bound
through the red path and the blue edge that the algorithm does is valid, and
there is no further need for an adjustment due to the green path.

Taking into account the return to the zero-state, which is important for
$k>1$, one can adapt the proof of Larsen also for the new algorithm proposed
in Section VI.

## Acknowledgment

This work is supported by the SNSF grant n. 212865, by CIDMA through FCT,
https://doi.org/10.54499/UIDB/04106/2020,
https://doi.org/10.54499/UIDP/04106/2020 and by FCT grant UI/BD/151186/2021.

## References

  * [1] R. Johannesson and K. S. Zigangirov, Fundamentals of Convolutional Coding, ser. Digital and Mobile Communication. New Jersey: Wiley-IEEE Press, 1999.
  * [2] J. Lieb, R. Pinto, and J. Rosenthal, Convolutional Codes. CRC Press, 2021.
  * [3] O. David and V. Lyandres, “Heapmod algorithm for computing the minimum free distance of convolutional codes,” in 21st IEEE Convention of the Electrical and Electronic Engineers in Israel. Proceedings. IEEE, 2000, pp. 435–438.
  * [4] M. Rouanne and J. Daniel J. Costello, “An algorithm for computing the distance spectrum of trellis codes,” IEEE Journal on Selected Areas in Communications, vol. 7, no. 6, p. 929, August 1989.
  * [5] K. Larsen, “Comments on "an efficient algorithm for computing free distance",” IEEE Transactions on Information Theory, vol. 19, no. 4, pp. 577–579, 1973.
  * [6] L. Bahl, C. Cullum, W. Frazer, and F. Jelinek, “An efficient algorithm for computing free distance,” IEEE Transactions on Information Theory, vol. 18, no. 3, pp. 437–439, 1972.
  * [7] M. L. Cedervall and R. Johannesson, “A fast algorithm for computing distance spectrum of convolutional codes,” IEEE Transactions on Information Theory, vol. 35, no. 6, pp. 1146–1159, 1989.
  * [8] S. Lin and D. J. Costello, Error Control Coding: Fundamentals and Applications. Prentice Hall Professional Technical Reference, 1994.
  * [9] J. Rosenthal and R. Smarandache, “Maximum distance separable convolutional codes,” Appl. Algebra Engrg. Comm. Comput., vol. 10, no. 1, pp. 15–32, 1999.
  * [10] D. J. Costello, “A construction technique for random-error-correcting convolutional codes,” IEEE Transactions on Information Theory, vol. 15, no. 5, pp. 631–636, 1969.
  * [11] A. Miczo and L. D. Rudolph, “A note on the free distance of a convolutional code,” IEEE Trans. Inf. Theory, vol. 16, no. 5, pp. 646–648, 1970.
  * [12] J. Lieb and J. Rosenthal, “Erasure decoding of convolutional codes using first-order representations,” Math. Control Signals Syst., vol. 33, pp. 499–513, 2021


Transverse-momentum-dependent parton distribution functions
for spin-1 hadrons

S. Kumano1,2$\star$ and Qin-Tao Song3

1 KEK Theory Center, Institute of Particle and Nuclear Studies, KEK,

Oho 1-1, Tsukuba, Ibaraki, 305-0801, Japan

2 J-PARC Branch, KEK Theory Center, Institute of Particle and Nuclear Studies,
KEK,

and Theory Group, Particle and Nuclear Physics Division, J-PARC Center,

Shirakata 203-1, Tokai, Ibaraki, 319-1106, Japan

3 School of Physics and Microelectronics, Zhengzhou University,

Zhengzhou, Henan 450001, China

*<EMAIL_ADDRESS><EMAIL_ADDRESS>
Proceedings for the XXVIII International Workshop
on Deep-Inelastic Scattering and Related Subjects,
Stony Brook University, New York, USA, 12-16 April 2021
10.21468/SciPostPhysProc.?

## Abstract

We show transverse-momentum-dependent parton distribution functions (TMDs) for
spin-1 hadrons including twist-3 and 4 functions by taking the decomposition
of a quark correlation function in the Lorentz-invariant way with the
conditions of Hermiticity and parity invariance. We found 30 new TMDs in the
tensor-polarized spin-1 hadron at twists 3 and 4 in addition to 10 TMDs at
twist 2. Since time-reversal-odd terms of the collinear correlation function
should vanish after integrals over the partonic transverse momentum, we
obtained new sum rules for the time-reversal-odd structure functions,
$\bm{\int d^{2}k_{T}g_{LT}=\int d^{2}k_{T}h_{LL}=\int d^{2}k_{T}h_{3LL}}$
$\bm{=0}$, at twists 3 and 4. We also indicated that transverse-momentum-
dependent fragmentation functions exist in tensor-polarized spin-1 hadrons.
The TMDs can probe color degrees of freedom, so that they are valuable in
providing unique opportunities for creating interdisciplinary physics fields
such as gluon condensate, color Aharonov-Bohm effect, and color entanglement.
We also found three new collinear PDFs at twists 3 and 4, and a twist-2
relation and a sum rule were derived in analogy to the Wandzura-Wilczek
relation and the Burkhardt-Cottingham sum rule on the structure function
$\bm{g_{2}}$.

##  1 Introduction

The field of transverse-momentum-dependent parton distribution functions
(TMDs) is one of rapidly-developing fields in hadron physics. It is
interesting because the explicit color degrees of freedom can be probed by the
TMDs. Depending on the color flow in hadrons, the TMDs have opposite signs,
for example, in semi-inclusive deep inelastic scattering and Drell-Yan
processes. The TMD physics is related to fundamental aspects of quantum
physics, such as color Aharonov-Bohm effect and color entanglement, and the
TMDs are also valuable for understanding the color glass condensate. The color
is confined in hadrons and it does not appear easily in experimental
observables. The TMDs provide a unique opportunity to shed light on new hadron
phenomena involving the color.

The TMDs have been investigated so far for the spin-1/2 nucleons. On the
collinear structure function $b_{1}$ of the spin-1 deuteron, there was a
measurement by the HERMES collaboration in 2005 [1]. In future, there are
experimental projects to investigate structure functions of spin-1 hadrons,
especially the spin-1 deuteron, at the Jefferson laboratory (JLab) [2],
Fermilab [3], and Nuclotron-based Ion Collider fAcility (NICA) [4], and
electron-ion colliders in US and China. We also have been investigating
structure functions of spin-1 hadrons theoretically by considering these
experimental projects [5, 6, 7, 8]. Therefore, it is a good opportunity and
timely to investigate TMDs also for the spin-1 hadrons. However, only the
twist-2 theoretical formalism was available [9]. The twist-3 and 4 parts were
investigated recently in Ref. [7], where new structure functions and sum rules
were proposed at twists 3 and 4. In addition, a twist-2 relation and a sum
rule were derived for the twist-3 collinear structure function $f_{LT}$ [8].
We explain these works in this paper.

##  2 Correlation functions and polarizations for spin-1 hadrons

The TMDs and collinear parton distribution functions (PDFs) are defined from
the correlation function given by the matrix element as

$\displaystyle\Phi_{ij}^{[c]}(k,P,T\,|\,n)$
$\displaystyle=\int\\!\frac{d^{4}\xi}{(2\pi)^{4}}\,e^{ik\cdot\xi}\langle\,P,T\left|\,\bar{\psi}_{j}(0)\,W^{[c]}(0,\xi)\psi_{i}(\xi)\,\right|P,\,T\,\rangle.$
(1)

Here, $k$ and $P$ are quark and hadron momenta, $T$ indicates the tensor
polarization of the hadron, $\xi$ is the space-time coordinate of the quark,
$\psi$ is the quark field, $W^{[c]}(0,\xi)$ is the gauge link, $c$ indicates
the integral path, and $n$ is the lighcone vector
$n^{\mu}=(1,0,0,-1)/\sqrt{2}$. Since only the tensor polarization is
considered in this work, the vector polarization $S$ is not explicitly
written. The TMD and collinear correlation functions are defined by
integrating the quark transverse momenta as

$\displaystyle\Phi^{[c]}(x,k_{T},P,T)$
$\displaystyle=\\!\int\\!dk^{+}dk^{-}\,\Phi^{[c]}(k,P,T\,|n)\,\delta(k^{+}\\!-xP^{+}),$
(2) $\displaystyle\Phi(x,P,T)$
$\displaystyle=\\!\int\\!d^{2}k_{T}\,\Phi^{[c]}(x,k_{T},P,T).$ (3)

The covariant form of the tensor polarization $T^{\mu\nu}$ of a spin-1 hadron
is expressed by the tensor-polarization parameters $S_{T}^{x}$, $S_{T}^{y}$,
$S_{L}$, $S_{LL}$, $S_{TT}^{xx}$, $S_{TT}^{xy}$, $S_{LT}^{x}$, and
$S_{LT}^{y}$ as

$\displaystyle T^{\mu\nu}=\frac{1}{2}$
$\displaystyle\left[\frac{4}{3}S_{LL}\frac{(P^{+})^{2}}{M^{2}}\bar{n}^{\mu}\bar{n}^{\nu}-\frac{2}{3}S_{LL}(\bar{n}^{\\{\mu}n^{\nu\\}}-g_{T}^{\mu\nu})\right.$
$\displaystyle\left.+\frac{1}{3}S_{LL}\frac{M^{2}}{(P^{+})^{2}}n^{\mu}n^{\nu}+\frac{P^{+}}{M}\bar{n}^{\\{\mu}S_{LT}^{\nu\\}}-\frac{M}{2P^{+}}n^{\\{\mu}S_{LT}^{\nu\\}}+S_{TT}^{\mu\nu}\right].$
(4)

Here, $a^{\\{\mu}b^{\nu\\}}$ indicates the symmetrized combination
$a^{\\{\mu}b^{\nu\\}}=a^{\mu}b^{\nu}+a^{\nu}b^{\mu}$, $M$ is the hadron mass,
and $P^{+}$ is the lightcone momentum given by $P^{+}=(P^{0}+P^{3})/\sqrt{2}$.
Using this tensor together with available Lorentz vectors $P$, $k$, and $n$,
we wrote the general form of the correlation function by considering the
condition of Hermiticity and parity invariance as [7]

$\displaystyle\Phi(k,P,T\,|n)$
$\displaystyle=\frac{A_{13}}{M}T_{kk}+\frac{A_{14}}{M^{2}}T_{kk}\not{P}+\cdots+\frac{A_{20}}{M^{2}}\varepsilon^{\mu\nu
Pk}\gamma_{\mu}\gamma_{5}T_{\nu k}$ $\displaystyle+\frac{B_{21}M}{P\cdot
n}T_{kn}+\frac{B_{22}M^{3}}{(P\cdot n)^{2}}T_{nn}+\cdots+\frac{B_{52}M}{P\cdot
n}\sigma_{\mu k}T^{\mu n},$ (5)

where the notation $X_{\mu k}\equiv X_{\mu\nu}k^{\nu}$ is used for the
contraction of the tensor $X_{\mu\nu}$ with $k^{\nu}$. The important point of
this expression is that the terms with the lightcone vector $n$ are included
for finding the twist-3 and 4 structure functions, as investigated in finding
twist-3 and 4 structure functions in the spin-1/2 nucleons [10]. Integrating
this expression as given in Eq. (2), we obtained possible tensor-polarized
structure functions up to twist 4.

##  3 TMDs and PDFs of spin-1 hadrons up to twist 4

TMDs for tensor-polarized spin-1 hadrons

Table 1:  Twist-2 TMDs.

Table 2:  Twist-3 TMDs.

Table 3:  Twist-4 TMDs.

The TMDs are defined from the correlation functions of Eqs. (2) and (5) by
taking traces with various $\gamma$ matrices as
$\Phi^{\left[\Gamma\right]}(x,k_{T},T)\equiv\frac{1}{2}\,\text{Tr}\left[\,\Phi(x,k_{T},T)\Gamma\,\right],$
and possible TMDs are listed in Tables 3, 3, and 3. The upper half sections
(unpolarized U, longitudinally polarized L, transversely polarized T) of these
tables are the same as the ones of the spin-1/2 nucleons. In this work, the
tensor-polarized sections (LL, LT, TT) are studied. The square brackets $[\ ]$
indicate chiral-odd distributions and the others are chiral-even ones. The
T-even (T: time reversal) and T-odd distributions are separately written. The
twist-2 part was studied in Ref. [9] by calculating $\Phi^{[\gamma^{+}]}$,
$\Phi^{[\gamma^{+}\gamma_{5}]}$, and $\Phi^{[i\sigma^{i+}\gamma_{5}]}$
($\Phi^{[\sigma^{i+}]}$). Our new results were for the twist-3 and twist-4
functions in Tables 3 and 3 [7].

Since the full expressions are lengthy, we explain the outline by taking an
example on the correlation function $\Phi^{[\gamma^{i}]}$ and related twist-3
TMDs. The twist-3 TMDs were obtained by calculating $\Phi^{[\gamma^{i}]}$,
$\Phi^{\left[{\bf 1}\right]}$, $\Phi^{\left[i\gamma_{5}\right]}$
$\Phi^{[\gamma^{i}\gamma_{5}]}$ $\Phi^{[\sigma^{ij}]}$, and
$\Phi^{[\sigma^{-+}]}$. Among them, $\Phi^{[\gamma^{i}]}$ is expressed by the
polarization parameters and TMDs as

$\displaystyle\Phi^{[\gamma^{i}]}(x,k_{T},T)=\frac{M}{P^{+}}\bigg{[}f^{\perp}_{LL}(x,k_{T}^{\,2})S_{LL}\frac{k_{T}^{i}}{M}\\!+\\!f^{\,\prime}_{LT}(x,k_{T}^{\,2})S_{LT}^{i}$
$\displaystyle\ \hskip 28.45274pt-
f_{LT}^{\perp}(x,k_{T}^{\,2})\frac{k_{T}^{i}S_{LT}\cdot
k_{T}}{M^{2}}-f_{TT}^{\,\prime}(x,k_{T}^{\,2})\frac{S_{TT}^{ij}k_{T\,j}}{M}$
$\displaystyle\ \hskip
28.45274pt+f_{TT}^{\perp}(x,k_{T}^{\,2})\frac{k_{T}\cdot S_{TT}\cdot
k_{T}}{M^{2}}\frac{k_{T}^{i}}{M}\bigg{]},$

where the TMDs $f^{\perp}_{LL}$, $\cdots$, $f_{TT}^{\perp}$ are expressed by
the expansion coefficients $A_{i}$ and $B_{i}$ of Eq. (5). The guideline for
assigning the TMD names $f^{\perp}_{LL}$, $\cdots$, $f_{TT}^{\perp}$ is
explained in Ref. [7] depending on polarizations and transverse-momentum
factors. In this way, we obtained five TMDs in $\Phi^{[\gamma^{i}]}$ and they
are listed in Table 3. In listing the TMDs in the tables, we redefined the
TMDs by $F(x,k_{T}^{\,2})\equiv
F^{\,\prime}(x,k_{T}^{\,2})-(k_{T}^{\,2}/(2M^{2}))\,F^{\perp}(x,k^{\,2}_{T})$
where $k_{T}^{\,2}=-\vec{k}_{T}^{\,2}$, and the functions without ′ are shown
in the tables. The other twist-3 TMDs are also listed in Table 3. In the same
way, the twist-4 TMDs are obtained by calculating $\Phi^{[\gamma^{-}]}$,
$\Phi^{[\gamma^{-}\gamma_{5}]}$, and $\Phi^{[\sigma^{i-}]}$ and they are
listed in Table 3. There are 40 TMDs in the tensor-polarized spin-1 hadron as
shown in Tables 3, 3, and 3. Among them, we found that 30 new TMDs exist at
twist-3 and 4 as

$\displaystyle\text{Twist-3 TMD:}\ \ f_{LL}^{\perp},\ e_{LL},\ f_{LT},\
f_{LT}^{\perp},\ e_{1T},\ e_{1T}^{\perp},\ f_{TT},\ f_{TT}^{\perp},\ e_{TT},\
e_{TT}^{\perp},\ g_{LL}^{\perp},\ g_{LT},\ g_{LT}^{\perp},$ $\displaystyle\
\hskip 65.44142ptg_{TT},\ g_{TT}^{\perp},\ h_{1L},\ h_{LT},\ h_{LT}^{\perp},\
h_{TT},\ h_{TT}^{\perp},$ $\displaystyle\text{Twist-4 TMD:}\ \ f_{3LL},\
f_{3LT},\ f_{3TT},\ g_{3LT},\ g_{3TT},\ h_{3LL}^{\perp},\ h_{3LT},\
h_{3LT}^{\perp},\ h_{3TT},\ h_{3TT}^{\perp}.$ (7)

They are classified by the time-reversal and chirality properties. Since the
T-odd collinear PDFs should vanish $f(x)_{\text{T-odd}}=0$ due to the time-
reversal invariance, we have the following sum rules at twist-3 and 4

$\displaystyle\\!\int\\!d^{2}k_{T}\,g_{LT}(x,k_{T}^{\,2})=\\!\int\\!d^{2}k_{T}\,h_{LL}(x,k_{T}^{\,2})=\\!\int\\!d^{2}k_{T}\,h_{3LT}(x,k_{T}^{\,2})=0.$
(8)

TMD fragmentation functions of spin-1 hadrons

New TMD fragmentation functions also exist for the spin-1 hadrons [7], and
they are obtained simply by changing the kinematical variables and function
names as

$\displaystyle\ \text{Kinematical variables:}\ \
x,k_{T},S,T,M,n,\gamma^{+},\sigma^{i+}\Rightarrow\
z,k_{T},S_{h},T_{h},M_{h},\bar{n},\gamma^{-},\sigma^{i-},$ $\displaystyle\
\text{Distribution functions:}\ \ f,g,h,e\hskip 65.44142pt\ \Rightarrow\
\text{Fragmentation functions:}\ D,G,H,E.$ (9)

Table 4:  Twist-2 PDFs.

Table 5:  Twist-3 PDFs.

Table 6:  Twist-4 PDFs.

PDFs for tensor-polarized spin-1 hadrons

The collinear PDFs of the tensor-polarized hadrons [11] are obtained by
integrating the TMDs over the transverse momentum as $f(x)=\int
d^{2}k_{T}f(x,k_{T}^{\,2})$. Many functions vanish by this integral and the
remaining PDFs are shown in Tables 6, 6, and 6. Therefore, the new twist-3 and
4 functions, which we found, are [7]

$\displaystyle\text{Twist-3 PDF:}\ e_{LL},\ f_{LT},\ \ \ \text{Twist-4 PDF:}\
f_{3LL}.$ (10)

The asterisks $*1$, $*2$, $*3$, $*4$ indicate the collinear PDFs $h_{1LT}(x)$,
$g_{LT}(x)$, $h_{LL}(x)$, $h_{3LT}(x)$ vanish, respectively, because of the
time-reversal invariance. However, since the time-reversal invariance cannot
be imposed in the fragmentation functions, the corresponding fragmentation
functions $H_{1LT}(z)$, $G_{LT}(z)$, $H_{LL}(z)$, $H_{3LT}(z)$, as indicated
by the replacements of Eq. (9), should exist as collinear fragmentation
functions [7, 12].

Twist-2 relation and sum rule for $\bm{f}_{LT}$

Since we obtained the new PDFs for spin-1 hadrons, it is possible to
investigate useful relations for them in the similar way to the Wandzura-
Wilczek (WW) relation and the Burkhardt-Cottingham (BC) sum rule. For finding
twist-2 relations, it is important to specify higher-twist effects. Twist-3
effects are investigated by the nonlocal operator
$\bar{\psi}(0)\big{(}\partial^{\mu}\gamma^{\alpha}-\partial^{\alpha}\gamma^{\mu}\big{)}\psi(\xi)$,
whose matrix element between tensor-polarized spin-1 hadron states is
expressed by $f_{1LL}$ and $f_{LT}$ up to twist 3. For the tensor-polarized
spin-1 hadron, twist-3 multiparton distribution functions

$\displaystyle F_{LT}(x,y),\ \ G_{LT}(x,y),\ \ H_{LL}^{\perp}(x,y),\ \
H_{TT}(x,y),$ (11)

exist in general [8]. The twist-3 matrix element is expressed by $F_{G,LT}$
and $G_{G,LT}$ defined with the gluon field tensor $G^{\mu\nu}$. Specifying
the twist-3 terms, we obtained [8]

$\displaystyle f_{LT}(x)$
$\displaystyle=\frac{3}{2}\int^{\epsilon(x)}_{x}dy\frac{f_{1LL}(y)}{y}+\int^{\epsilon(x)}_{x}dy\frac{f_{LT}^{(HT)}(y)}{y},$
(12)

where the higher-twist function $f_{LT}^{(HT)}(x)$ is expressed by the
integral of $F_{G,LT}$ and $G_{G,LT}$. Here, $\epsilon(x)=1\ (-1)$ for $x>0\
(x<0)$. In the positive (negative) $x$ region, $f(x)$ is a quark (antiquark)
distribution. We define $f^{+}$ distribution by $f^{+}(x)\equiv
f(x)+\bar{f}(x)=f(x)-f(-x)$ for $f=f_{1LL},\ f_{LT},\ f_{LT}^{(HT)}$ at $x>0$,
the relation becomes
$f_{LT}^{+}(x)=\frac{3}{2}\int^{1}_{x}dy\,f_{1LL}^{+}(y)/y$ if the higher-
twist terms are neglected. Furthermore, if the function $f_{2LT}$ is defined
by $f_{2LT}(x)\equiv(2/3)f_{LT}(x)\\\ -f_{1LL}(x)$, this integral relation
becomes

$\displaystyle
f_{2LT}^{+}(x)=-f_{1LL}^{+}(x)+\int^{1}_{x}\frac{dy}{y}f_{1LL}^{+}(y).$ (13)

We obtain a sum rule by integrating this twist-2 relation over $x$ as

$\displaystyle\int_{0}^{1}dx\,f_{2LT}^{+}(x)=0.$ (14)

These are analogous equations to the WW relation and BC sum rule. The function
$f_{1LL}$ is given by the tensor-polarized structure function $b_{1}$ as
$-(3/2)f_{1LL}^{+}=b_{1}^{q}+b_{1}^{\bar{q}}$. If the tensor-polarized
antiquark distributions vanish, the parton-model sum rule $\int dxb_{1}=0$
exists [5]. Then, there is a sum rule for the $f_{LT}$ itself as
$\int_{0}^{1}dx\,f_{LT}^{+}(x)=0$.

##  4 Conclusion

We proposed new TMD structure functions at twist 3 and 4 by decomposing the
TMD correlation function in the Lorentz invariant way with constraints of
Hermiticity and parity invariance. There are 40 TMDs for the tensor-polarized
spin-1 hadron up to twist 4, and we found new 30 TMDs at twist 3 and 4. Then,
we showed sum rules for the time-reversal-odd TMDs. There exist also
corrensponing fragmentation functions simply changing kinematical variables
and function names. Integrating the TMDs over the transverse momentum, we
obtained new collinear PDFs $e_{LL}$, $f_{LT}$, and $f_{3LL}$ at twist 3 and
4, in addition to the twist-2 function $f_{1LL}$ ($b_{1}$). We also indicated
the twist-2 relation and the sum rule for $f_{LT}$ in analogy to the Wandzura-
Wilczek relation and the Burkhardt-Cottingham sum rule.

## Acknowledgements

S. Kumano was partially supported by Japan Society for the Promotion of
Science (JSPS) Grants-in-Aid for Scientific Research (KAKENHI) Grant Number
19K03830. Qin-Tao Song was supported by the National Natural Science
Foundation of China under Grant Number 12005191 and the Academic Improvement
Project of Zhengzhou University.

## References

  * [1] A. Airapetian et al. (HERMES Collaboration), Phys. Rev. Lett. 95 (2005) 242001.
  * [2] J.-P. Chen et al., Proposal to Jefferson Lab PAC-38, PR12-11-110 (2011); M. Jones et al., A Letter of Intent to Jefferson Lab PAC 44, LOI12-16-006 (2016).
  * [3] D. Geesaman et al., Proposal to Fermilab PAC, P-1039 (2013).
  * [4] A. Arbuzov et al., Prog. Nucl. Part. Phys. 119, 103858 (2021) 1.
  * [5] F. E. Close and S. Kumano, Phys. Rev. D 42 (1990) 2377.
  * [6] S. Hino and S. Kumano, Phys. Rev. D 59 (1999) 094026; D 60 (1999) 054018; T.-Y. Kimura and S. Kumano, Phys. Rev. D 78 (2008) 117505; S. Kumano, Phys. Rev. D 82 (2010) 017501; W. Cosyn, Yu-Bing Dong, S. Kumano, and M. Sargsian, Phys. Rev. D 95 (2017) 074036; S. Kumano and Qin-Tao Song, Phys. Rev. D 94 (2016) 054022; D 101 (2020) 054011 & 094013\.
  * [7] S. Kumano and Qin-Tao Song, Phys. Rev. D 103 (2021) 014025.
  * [8] S. Kumano and Qin-Tao Song, arXiv:2106.15849, submitted for publication.
  * [9] A. Bacchetta and P. Mulders Phys. Rev. D 62 (2000) 114004.
  * [10] K. Goeke, A. Metz and M. Schlegel, Phys. Lett. B 618 (2005) 90; A. Metz, P. Schweitzer and T. Teckentrup, Phys. Lett. B 680 (2009) 141.
  * [11] P. Hoodbhoy, R. L. Jaffe, and A. Manohar, Nucl. Phys. B 312 (1989) 571; J. P. Ma, C. Wang, and G. P. Zhang, arXiv:1306.6693 (2013).
  * [12] X. Ji, Phys. Rev. D 49 (1994) 114.


# Doppler Signature of a Possible Termination Shock in an Off-Limb Solar Flare

Ryan J. French,1 Sijie Yu,2 Bin Chen,2 Chengcai Shen3 and Sarah A. Matthews4
1National Solar Observatory, 3665 Discovery Drive, Boulder, CO 80303, USA
2Center for Solar-Terrestrial Research, New Jersey Institute of Technology, 3L
King Jr. Blvd., Newark, NJ 07102-1982, USA
3Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA
02138, USA
4University College London, Mullard Space Science Laboratory, Holmbury St.
Mary, Dorking, Surrey, RH5 6NT, UK

(Accepted to MNRAS, Feb 6th 2024)

###### Abstract

We report striking Doppler velocity gradients observed during the well-
observed September 10th 2017 solar flare, and argue that they are consistent
with the presence of an above-the-looptop termination shock beneath the flare
current sheet. Observations from the Hinode Extreme-ultraviolet Imaging
Spectrometer (EIS) measure plasma sheet Doppler shifts up to 35 km/s during
the late-phase of the event. By comparing these line-of-sight flows with
plane-of-sky measurements, we calculate total velocity downflows of 200+ km/s,
orientated $\approx 6-10^{\circ}$ out of the plane of sky. The observed
velocities drop rapidly at the base of the hot plasma sheet seen in extreme
ultraviolet, consistent with simulated velocity profiles predicted by our
2.5-D magnetohydrodynamics model that features a termination shock at the same
location. Finally, the striking velocity deceleration aligns spatially with
the suppression of Fe XXIV non-thermal velocities, and a 35–50 keV hard X-ray
looptop source observed by the Reuven Ramaty High Energy Solar Spectroscopic
Imager (RHESSI). Together, these observations are consistent with the presence
of a possible termination shock within the X8.2-class solar flare.

###### keywords:

Solar flare – Termination Shock – Current Sheet

††pubyear: 2024††pagerange: Doppler Signature of a Possible Termination Shock
in an Off-Limb Solar Flare–References

##  1 Introduction

The standard model of eruptive solar flares, also known as the CSHKP model
after the original works from Carmichael (1964); Sturrock (1968); Hirayama
(1974); Kopp & Pneuman (1976), is able to explain many observable phenomena in
solar flares. In this model, magnetic reconnection is driven by the inflow of
oppositely oriented magnetic fields beneath an erupting magnetic flux rope.
Within this reconnection site, magnetic free energy is converted into plasma
heating, particle acceleration, and electromagnetic radiation. Newly formed
magnetic fields are ejected from either end of the current sheet site,
contributing towards the erupting flux rope (and associated coronal mass
ejection) above, and forming hot flare loops, below.

However, the standard model is unable to account for some of the higher-energy
features observed in flares. In particular, it does not address the specific
physical mechanisms through which magnetic energy is converted into the
kinetic energy of accelerated particles observed in solar flares. In recent
years, proposed additions to the simplified 2-dimensional standard flare model
have been developed to explain observations of high energy particles. The
cartoon in Figure 3C of Chen et al. (2015) presents such a scenario after
Masuda et al. (1994), incorporating hard X-ray (HXR) and radio observations
into the standard model configuration, explained by the presence of a
termination shock. In solar flares, a termination shock will persist as super-
magnetosonic reconnection outflows from the reconnection site, collide with
the dense post-reconnected flare loops. Termination shocks of similar magnetic
geometry are also believed to occur in other astrophysical plasmas, such as in
supernova remnants (e.g. Miles, 2009).

However, convincing evidence of solar flare termination shocks are
comparatively rare. Chen et al. (2015) and Chen et al. (2019) use high-cadence
radio imaging spectroscopy to reveal the morphology and dynamics of a
termination shock accelerating solar flare energetic electrons, using
observations of a C1.9 class eruptive flare from the Karl G. Jansky Very Large
Array (VLA). Work from Luo et al. (2021) finds similar results, with VLA
observations of a much larger M8.4-class flare. The synergy of ultraviolet and
X-ray observations have also been used to deduct the likely presence of a
termination shock, in work from Polito et al. (2018). Polito et al. (2018)
detect large (200+ km/s) Doppler shifts in the high temperature Fe XXI line
with the Interface Region Imaging Spectrograph (IRIS, De Pontieu et al.,
2014). These Doppler shifts, observed in an on-disk X1-class flare, were found
to diverge above the flare loop arcade, aligning with a high energy 35–70 keV
HXR source observed by the Reuven Ramaty High Energy Solar Spectroscopic
Imager (RHESSI, Lin et al., 2002). The authors attributed these observations
to deflecting flows consistent with the existence of a termination shock above
the looptop. These studies all provide observational evidence that strongly
supports the existence of a termination shock, through measurements of ion
dynamics and electron populations. Simulations have also been able to produce
termination shocks within solar flares, first in 2D (Forbes, 1986; Takasao et
al., 2015; Shen et al., 2018) and more recently 3D (Shen et al., 2022; Shibata
et al., 2023).

According to predicted termination shock geometry, solar flare termination
shocks are sat within hot cusp-like features, above the bright flare loops,
but below the current sheet base. These hot cusp shapes are observed
frequently in solar flares at the limb of the Sun, easy to spot in the hot 131
and 94 Å channels of the Solar Dynamics Observatory’s Atmospheric Imaging
Assembly (SDO/AIA, Lemen et al., 2012). The tip of the cusp meets where we
expect the base of the current sheet to be, but observations of visible heated
plasma around the current sheet (known as a plasma sheet), are far more
infrequent. Because of their infrequency, spectroscopic measurements of the
plasma sheet and cusp region in (E)UV are even rarer. Only a handful of
spectroscopic observations of plasma sheets or cusp regions exist for solar
flares, including the famous September 10th 2017 flare. Many facets of this
flare have been studied, including dynamics of the hot, long-lived plasma
sheet (e.g, Warren et al., 2018; Cheng et al., 2018; Longcope et al., 2018;
French et al., 2019; Chen et al., 2020a), flare loop-top region (Chen et al.,
2020a; Fleishman et al., 2020; Yu et al., 2020), and flare arcade to the south
(Reeves et al., 2020) that is associated with the 3D geometry of the erupting
flux rope (Chen et al., 2020b). In this study, we revisit spectral
observations from the Hinode EUV Imaging Spectrometer (EIS, Culhane et al.,
2007) of the September 10th 2017 flare, in particular the late phase
observations presented by French et al. (2020). We combine EIS Doppler
measurements along the plasma sheet with plane-of-sky (POS) EUV downflows
detected within the same region and time period by Yu et al. (2020). Together,
the Doppler (line-of-sight, LOS) and POS downflows provide a 3D measurement of
the velocity profile, which we compare with velocities predicted by
termination shock simulations.

##  2 Summary of Observations

Figure 1: Left column: AIA 193 Å images at 16:09, 16:44 and 18:39 UT. Times
are consistent across each row. Middle column: Hinode-EIS Fe XXIV 192.04 Å
intensity. Right column: Hinode-EIS Fe XXIV 192.04 Å Doppler velocity. Green
and cyan contours provide spatial locations of highest non-thermal velocity
and Doppler velocity respectively. Doppler velocity contours are set at $-$16
km/s, and non-thermal velocity contours 100, 70, & 60 km/s for 16:09, 16:44
and 18:39 UT respectively. Doppler contours are not included in the 16:09 UT
frame, due to PSF effects. For the 16:44 UT maps, RHESSI contours of 6–12 keV
and 35–50 keV HXR emission are also included as blue and red contours
respectively. The four gray horizontal-dashed lines at 16:44 and 18:39 UT mark
the Narrow (short dashes) and Wide (long dashes) plasma sheet cross-sections
analyzed in Figure 2.

The X8.2-class September 10th 2017 flare erupted over the western limb from AR
12673. Due to its position on the Sun and within the solar cycle, the
September 10th 2017 flare was observed across the electromagnetic spectrum by
a fleet of space and ground-based instruments, making it the most published
solar flare of all time (French, 2022).

The first column of Figure 1 shows AIA 193 Å emission of the flare over the
west limb, at snapshots coinciding with the start of three Hinode EIS rasters.
The AIA 193 Å broadband filter contains two strong spectral lines, the cool Fe
XII 193 line (peaking at 1.6 MK) and hot Fe XXIV 192.04 Å line (peaking at 18
MK). The hot Fe XXIV line is only produced in the strongest of flares, so not
common in AIA 193 Å images. The three rows of Figure 1 show the evolution of
AIA 193 Å emission, with growth of the flare loops and hot cusp region, and
persistence of the hot, horizontal plasma sheet.

Hinode-EIS observed the September 10th 2017 flare for long periods of its
duration, including from flare onset at 15:44 UT, through the peak at 16:07
UT, up to 16:52 UT. Observations continued later in the flare’s evolution,
from 18:39 UT to 19:31 UT. EIS observed a 2″ slit-width raster scan, rastering
from west to east with a cadence of 8 minutes 52 seconds. EIS measured spectra
of multiple spectral lines, including the hot Fe XXIV 192.04 Å line included
in the broadband AIA 193 Å window. The hot formation temperature of Fe XXIV
renders it useful for studies of hot flaring plasma. The second column of
Figure 1 presents intensity measurements of EIS Fe XXIV 192.04 Å aligned with
AIA, obtained by fitting a single Gaussian curve to 1x3 binned data (to create
square 3x3″pixels). The three EIS intensity maps presented in Figure 1 were
chosen to show the first EIS raster of the plasma sheet at 16:09 UT, the final
EIS raster of the first observing sequence (16:44 UT), and first EIS raster of
the continued observations at 18:39 UT. That is, EIS was collecting 9 minute
rasters between the first two frames of Figure 1, but with a data gap between
the second and third frame. In these images, the solar limb is just off the
left edge of the image field-of-view (FOV). There are clear similarities
between the EIS 192.04 Å and AIA 193 Å emission, specifically in the presence
of the hot plasma sheet and above-looptop cusp region. All EIS maps have been
cropped from their initial FOV. Due to a shift in the EIS pointing from 18:39
UT onwards, the right edge of the 18:39 UT FOV sits within the plotted region.

Gaussian fitting of the Fe XXIV emission line also provides diagnostics of
Doppler velocity and non-thermal velocity (excess line broadening beyond the
thermal width) for the flaring-temperature plasma. The third column of Figure
1 show maps of EIS Fe XXIV Doppler velocity, aligned to the other intensity
panels in the figure. Due to the shape of the point-spread function (PSF),
care is needed when interpreting off-limb Doppler measurements, which can
generate artifacts along the slit direction. These artifacts are most
prominent at sharp boundaries between high and low intensity regions, such as
that between the bright off-limb features and the dark background of space.
This PSF has a stark effect in the early EIS Doppler velocity maps of the
flare examined here, within the bright plasma sheet. In Figure 1, the slit
orientation is top to bottom, rastering from right to left. In the bright
plasma sheet region, photons from bright pixels can contaminate the wings of
adjacent pixels, creating an artificial Doppler signature. This PSF effect is
discussed in detail by Warren et al. (2018), and easily seen in the 16:09 UT
Doppler velocity frame in Figure 1. In this frame, horizontal streaks of
strong parallel red and blue shifts are seen along the plasma sheets. These
Doppler measurements are not real, and an artifact of the EIS PSF. Later in
the flares evolution, growth of the cusp region around the plasma sheet
removes the sharp intensity transition, and therefore reduces influence of the
PSF on the Doppler maps. For these later Doppler maps (as presented in French
et al., 2020), we can reliably extract LOS velocity information from the
Doppler profiles. Examining the Doppler maps at 16:44 and 18:39 in Figure 1,
this study takes particular interest at the strong blue shift signal within
the plasma sheet region towards the right of the frame, and strong velocity
gradient at the base of the plasma sheet. For these two later EIS rasters,
blue-shift contours (at $-$16 km/s) are overlaid onto the intensity maps in
left and center panels (in cyan). We use the standard definition for blue
shift in this study, defined as velocities flowing towards the observer. Blue
shift velocities at the base of the plasma sheet, where plasma is flowing
towards the Sun, immediately indicate that the plasma sheet is tilted slightly
away from us, out of the POS above the limb. This is consistent with the known
geometry of the September 10th 2017 flare (e.g. Chen et al., 2020b).

Also included within every panel in Figure 1 are green contours of EIS 192.04
Å non-thermal velocity. Values of 100, 70 and 60 km/s are chosen for the three
presented EIS times respectively. EIS non-thermal velocity measurements for
this flare have been examined previously, e.g. Warren et al. (2018); Longcope
et al. (2018) for the earlier EIS rasters, and French et al. (2020) for the
later EIS rasters. The non-thermal velocities have been interpreted in these
studies as macroscopic turbulence along the LOS, as expected from
reconnection-induced turbulence within the plasma sheet.

Intermittent measurements of HXRs were also observed by the RHESSI for this
flare, which had four operational detectors at the time (1, 3, 6, 8). For the
16:44 UT panels in Figure 1, RHESSI 6-12 and 35-50 keV emission are
overplotted with blue and red contours respectively. The RHESSI HXR images are
the same values presented in Figure 2e of (Gary et al., 2018), and were
reconstructed using the clean algorithm with data from detectors 3, 6 and 8.
The 6–12 keV signal, consisting primarily of thermal emission, is located at
the top of the flare loops. The 35–50 keV signal, however, containing mainly
non-thermal emission, peaks higher above the limb at the intersection between
the hot flare cusp and plasma sheet. Unfortunately, RHESSI was not observing
during the later 18:39 UT period.

##  3 velocity gradients along the plasma sheet

###  3.1 Line-of-Sight Velocities

Figure 2: Plane-of-sky gradients of EIS 192.04 Å Doppler and non-thermal
velocity, relative to 35-50 keV RHESSI emission. The RHESSI emission curve is
normalised to the maximum/minimum y axis values. The two columns present the
cross-section of these values along the narrow and wide cross-sections
presented in Figure 1, respectively. The top and bottom row each show the same
cross-sections at 16:44 and 18:39 UT.

Dashed horizontal gray lines in Figure 1 mark the locations of two cross-
sections along the plasma sheet. The outer dashed lines mark the Wide cross-
section, which encompasses all EUV emission from the current sheet. The inner
dashed lines mark the location of the Narrow cross-section, centered within
the strong blue-shift signal in the plasma sheet. Along these regions, we take
the cross-section of EIS Fe XXIV Doppler and non-thermal velocities, with
RHESSI 35-50 keV emission. Values are averaged in the y-direction, across the
cross-section region width.

Figure 2 plots both the Narrow and Wide cross-sections for the 16:44 and 18:39
UT time frames. Doppler velocity values are plotted on the left y-axis, and
non-thermal velocity values on the right y-axis. The Doppler velocity shows
blue shift, but we plot the magnitude of the blue shift for ease of comparison
to the non-thermal velocity trends.

All four panels demonstrate a strong gradient in both Doppler and non-thermal
velocities, suggesting a rapid deceleration of plasma as it approaches around
$x=1020\arcsec$. Both EIS Fe XXIV parameters then experience another gradual
or rapid decline below $x=1020\arcsec$ at lower altitudes, but it is this
first drop and local minimum we are most interested in here. The gradient is
most pronounced in the narrow cross-section, but still visible along the wider
region cross-section. Both Doppler and non-thermal velocity demonstrate
similar trends, especially in the wide cross-section, where the rate of
velocity decrease match for both time periods, (the lines are close to
parallel with the margin of error). The gradients between each parameter
differ more-so along the narrow-cross-section, resulting from the fact that
maximum Doppler and non-thermal velocities are positioned slightly differently
at a given width across the plasma sheet (Doppler velocities are centered off-
center across the plasma-sheet, whilst the non-thermal velocity shows a wider,
bifurcated pattern). Within the narrow cross-section, blue shift Doppler
velocities decrease from 35 to 0 km/s in under 60″ down the plasma sheet,
during both raster times.

At 16:44 UT, RHESSI observations are available for comparison. At this time,
the strong EIS Fe XXIV velocity gradients disappear at the same location as
peak high energy 35-50 keV RHESSI emission. A sudden velocity drop and a HXR
source located above the flare arcade are independently two key signatures
expected from the presence of a termination shock, which we see here at both
$\approx$ 1 $\&$ 3 hours into the flare evolution respectively.

###  3.2 Plane-of-Sky velocities

Figure 3: Plasma sheet downflows and upflows along the September 10th 2017
flare, detected in difference imaging of AIA 131 Å observations, (modified
from Figure 8C in Yu et al., 2020). Red lines mark upflows, red downflows, and
green selected downflows near EIS raster times. X’s mark the X-points
identified between upflow/downflow pairs. Plane-of-sky velocities are
distributed between 100-900 km/s, with an average of 250 km/s.

In Yu et al. (2020), the authors use running difference analysis of 131 Å
images to detect faint downflows and upflows along the September 10th 2017
plasma sheet. In doing this, a number of reconnection X-points are located
along the plasma sheet structure during this later-phase of the flare. Figure
3 shows a modified version of Figure 8c from Yu et al. (2020), displaying the
timings and height above the limb of all detected downflows, upflows and
X-points. Upflows and downflows are plotted in blue and red respectively, with
the exception of the four downflows closest to each EIS raster time, which are
plotted in green. EIS raster times are marked with gray dashed lines.

By tracking these plasma flows, we can measure the plasma velocity within the
POS, complementary to the LOS velocities measured by EIS Fe XXIV Doppler
values. The AIA 131 Å channel samples the hot Fe XXI line, which has a peak
temperature sensitivity of 10 MK. This is cooler than the Fe XXIV line
measured by EIS, which has a peak temperature sensitivity of 18 MK. Despite
this temperature difference, both instruments are observing hot flaring plasma
within the plasma sheet, and thus an empirical comparison of POS and LOS
velocities are valid.

The downflows plotted in Figure 3 have POS velocities ranging from 100-900
km/s, with an average of 250 km/s. Using basic right-angle trigonometry, we
can compare the observed 250 km/s POS velocity and 35 km/s Doppler velocity to
provide an initial estimate for the magnitude of the total plasma downflow
velocity vector. A trigonometric comparison with these values yields a plasma
sheet orientation 8∘ away from the POS, with 99% of the full velocity vector
captured within the observed POS downflows. To take a deeper investigation
into the velocity vector and how it changes along the plasma sheet, we can
compare observations to simulated plasma velocities from MHD simulations.

##  4 Comparison to Simulations

Figure 4: Left: Contribution function for Fe XXIV 192.02 Å. Center:
Temperature map of the MHD model frame presented in Figure 5(top). Right: Fe
XXIV temperature mask for temperature map in center frame.

###  4.1 Model Scaling and Temperature Mask

We compare observed velocities in the LOS and POS with the simulated velocity
profile of a termination shock, utilizing the 2.5D Athena flare model (Shen et
al., 2018). The model is dimensionless, and so scaling of the simulation
output is needed to compare with observations. Scales for four primary
characteristics are defined, and used to calculate the remaining dependent
characteristics. The first three primary characteristics are plasma beta
($=0.1$), ambient temperature ($=2$ MK) and magnetic field scaling (40 G), set
as standard values for the flaring corona. The fourth primary scaling
characteristic is length scale, which we set to match the limb to looptop
distance between the simulation and observations, independently for the 16:44
and 18:39 UT time periods of the September 10th 2017 flare. The dependent
characteristics, including pressure, density, temperature, current, time and
velocity, are scaled relative to the primary ones.

As a 2.5D flare simulation, parameters are all symmetric along the LOS axis.
To more accurately compare simulation LOS velocities with those measured in
the hot Fe XXIV line by EIS, we apply a temperature scaling mask to the
simulation output. Figure 4 (left) shows the contribution function of the Fe
XXIV 192.04 Å line, normalized to the maximum sensitivity of the line. By
cross-referencing the temperature sensitivity to the temperature map output by
the simulation (as an example shows in the center panel of Figure 4), we can
create a scaling mask to determine where in the simulation map contains
temperatures analogous to those observed by Fe XXIV. This is a quick empirical
method to remove flows of cooler plasma from the velocity comparison to
observations. The right panel of Figure 4 shows an example temperature mask
for one simulation frame, essentially removing the contribution of pixels
outside of the plasma sheet, cusp, or flare loops.

###  4.2 Line-of-Sight Comparison

Given the geometry of the flare (see Shen et al., 2018, for simualtion
information), the flare is positioned in the POS, with plasma downflows
primary in this plane. To replicate the LOS velocities observed by EIS, we
tilt the simulation map out of the POS, simultaneously scaling the LOS and POS
velocity components by a small angle. Given the small angles involved, it is
not required to alter the POS length scales during this process, given that an
angle of $10^{\circ}$ changes the length scale by $<1\%$.

Panel A of Figure 5 shows the LOS velocity map of a simulation frame late into
the impulsive phase, with an implemented tilt of 10∘ away from the POS. The Fe
XXIV temperature mask has been applied to this map, to empirically replicate
the Fe XXIV Doppler velocity observed by EIS. For this simulation frame, we
apply a length scaling to create a flare size close the September 10th 2017
flare at the time of the 16:44 UT EIS raster. As described earlier in this
section, dependent characteristics (including the characteristic velocity) are
then set relative to this length scale and other primary characteristics. This
model frame has a termination shock located at the horizontal gray dashed line
($\approx$95 Mm above the limb), determined by a minimum in the divergence of
velocity field (Shen et al., 2018).

Qualitatively comparing the simulated LOS velocities in Figure 5A with the EIS
Fe XXIV 192.04 Å Doppler maps in Figure 1, we immediately notice the strong
blue shift along the plasma sheet. To compare the two quantitatively, Figure
5B plots cross-sections of LOS velocity for both the simulation and EIS
observations. The solid black lines show the cross-section of simulated LOS
velocity, taken between the vertical gray dashed lined in panel A. We plot
this velocity cross-section for three different tilts, of angles 6 to 10∘ away
from the POS. The cyan curve plots the same EIS Doppler narrow cross-section
plotted in Figure 2. Please note that the axes have been flipped from this
earlier Figure 2, with distance plotted along the y-axis, and units of
arcseconds converted to Mm.

Figure 5B shows that the observed EIS velocity profile aligns nicely with the
simulated Doppler velocities (vdoppler) of tilts -6 to -10∘, matching the
gradient along the plasma sheet expected from the presence of a termination
shock. This comparison was completed for multiple frames from the simulation
evolution, with similar results. All of these frames included a termination
shock at the base of the plasma sheet.

Figure 5: Comparison of line-of-sight (LOS) and plane-of-sky (POS) velocities
between a simulated solar flare (2.5 ATHENA model) and observations. The top
and bottom row present a comparison to the 16:44 and 18:39 UT time periods of
the September 10th 2017 flare respectively. The first column presents the LOS
velocity map for the flare simulation, tilted 10∘ out of the POS and with a
temperature mask applied (to replicate Fe XXIV emission). The vertical dashed
lines mark the location of the plasma sheet cross-section, and horizontal
dashed line the termination shock location, determined by the local minimum in
the divergence of the velocity field. The second column plots cross-sections
along the plasma sheet, comparing simulation (vdoppler) and observed (EIS Fe
XXIV Doppler) LOS velocities. Simulation velocities are plotted at different
tilt angles from the LOS. The third column presents the cross-section of
velocity down the plasma sheet (in the POS), alongside the total simulation
velocity and Alfvén speed. The Alfvén speed is measured on the outside of the
current sheet. The height and speeds of observed POS downflows (from AIA 131 Å
observations) are also plotted for comparison. The fourth column plots the
cross-section of the Mach number along the model current sheet.

The second row of Figure 5 repeats analysis of the first column, instead
comparing the 18:39 UT EIS velocity observations to another simulation frame,
later in the simulated flare evolution. To match up the length scales between
the two, the characteristic length variable is increased by 25% from the
simulation frame analyzed in the first row, to account for the growth in the
flare arcade observed in reality (which is not present over the shorter
timescales of the simulation).

Given the shift in the FOV between the 16:44 and 18:39 UT, we do not capture
the velocity profile at such high altitudes in the plasma sheet as at 18.39 UT
as we did at 16:44 UT. Without observing the height of maximum velocity at
this time, a comparison to simulations contains slightly more uncertainty than
at earlier times – with just the rising velocity at lower altitudes captured.
Nonetheless, the velocity gradient at the observed section of the plasma sheet
aligns nicely with the simulation within the same angle tilt range. Subsequent
EIS rasters beyond 18:39 UT contain less and less of the velocity profile, as
the loop system continues to rise out of the EIS FOV.

###  4.3 Plane-of-Sky Comparison

In order to fully compare the observed and simulated velocity field, we must
also consider the POS velocity. Panel Figure 5C and G compare velocities in
the POS. The solid black line shows the simulated velocity along the same
plasma sheet cross-section (vertical gray lines in panel A), plotting velocity
in the Y direction instead of LOS direction (shown in panel B). A tilt of 10∘
is used to calculate this $v_{y}$ profile. The solid red line shows the total
velocity magnitude along the same slit, containing both the $v_{y}$ and LOS
components of the simulation plasma sheet velocity. Due to the small angle
involved ($\cos{[10^{\circ}]}\approx 98.5\%$), the total velocity magnitude is
only marginally larger than the $v_{y}$ component. Similarly to the LOS
velocity, the strong velocity gradient, reaching a minimum at the termination
shock in the simulation, is also present in the POS velocity.

Also plotted for reference is the Alvén velocity ($v_{a}$) adjacent to the
current sheet cross-section (colored magenta). In the adjacent panel, we plot
the cross-section of the Mach number along the current sheet,
$M=v_{t}/(c_{s}^{2}+v_{a}^{2})^{0.5}$, where $v_{t}$ is total velocity and
$c_{s}$ the local sound speed. Comparing current sheet downflow velocities
with $v_{a}$ and $M$, we see plasma velocities start to decrease at the
termination shock location, at an altitude of decreasing $v_{a}$ and $M>1$.
$v_{a}$ then reaches a local minimum as $M$ falls below unity. Although the
model provides a lower $v_{a}$ than we might initially expect, it is important
to note that $v_{a}$ is larger outside of the current sheet than within it.
Crucially, despite these lower values at the direct edge of the current sheet,
the Mach number remains greater than unity (as required to create the presence
of a shock).

To compare simulated POS velocities with observations, we overplot in Figure
5C and G the velocity and heights of plasma sheet downflows presented in
Figure 3. We use only the four downflows with times closest to each EIS
raster. These downflows are plotted in Figure 5 as vertical dashed lines, with
a color corresponding to the downflow time. For 5C, where comparison is made
to the 16:44 UT time period, all four observed downflows intersect the
simulated velocity profile. The same goes for the 18:39 UT comparison in 5G,
with the exception of a single, much faster downflow, which is off the scale.
The intersection of the observed downflows lines and simulated velocities
demonstrate that the simulation, with a velocity profile determined by the
presence of a termination shock, is successfully able to replicate both LOS
and POS observed velocities. It is important to note that the absence of
visible downflows at certain altitudes in AIA 131 Å observations does not
necessitate that downflows are not present at these locations and times. As
determined by the standard model, downflows are present continuously, with
only narrow windows of flowing plasma visible. With measured downflows of
100-170 and 150-210 km/s present around our two time windows (16:44 and 18:39
UT), the simulation comparison would suggest that faster downflows are
present, but not observed. The simulations indicate these could be as high as
210 and 240 km/s for each time frame respectively.

##  5 Conclusions

In this work, we have analyzed LOS and POS velocities observed along the
plasma sheet of an off-limb flare, at two time periods (16:44 and 18:39 UT)
during the flare’s late phase. By comparing LOS velocities of 35 km/s (from
EIS Fe XXIV 192.04 ÅDoppler observations) with POS velocities of 100-210 km/s
(from AIA 131 Å difference imaging), we deduce total velocities marginally
higher than the POS measurements, and a flare configuration oriented 6-10∘ out
of the POS. Although plasma downflows are continuous, our detections of them
are not. It is therefore probable that further downflows, both faster and
slower than those observed, are also present, but do not create the
temperature or density contrast necessary to detect in EUV images.

2.5D ATHENA simulations of an eruptive flare are able to replicate velocity
magnitudes and gradients in the LOS and POS, through the presence of a
termination shock at the base of the plasma sheet. Along the plasma sheet, the
simulation also suggests that downflow velocities are higher than we are able
to directly measure, possibly as high as 240 km/s during the 18:39 UT window.
With a termination shock needed to replicate the observed velocity gradients
within the simulations, this provides strong empirical evidence for the
existence of a termination shock.

In addition to the Fe XXIV Doppler velocities observed by EIS, the non-thermal
velocities also experience a sharp decline to a local minimum at the plasma-
sheet base. Non-thermal velocities in this flare have been previously
interpreted as a signal of macroscopic turbulence along the LOS, a key
prediction of ongoing reconnection (Warren et al., 2018; Longcope et al.,
2018). The simultaneous drop in the non-thermal velocities alongside the
Doppler velocities provides further evidence of a termination shock front
suppressing these flows. The non-thermal velocity variation along the current
sheet has been reported from 3D flare models in the past (e.g., Shen et al.,
2023). Models predict an increase in non-thermal velocity beneath the
termination shock due to the enhanced plasma turbulence in the below-shock
region, which is also consistent with observations (Figure 2, at heights below
$\sim$1015 and $\sim$1025″ at 16:44 and 18:39 UT respectively). However, the
deduced non-thermal velocity from observed Fe XXIV might be affected by
multiple-scale hot plasma flows along the LOS. Therefore, it is difficult to
compare quantitatively in the current study. Further investigation of non-
thermal velocities around the termination shock, and how the shock front
affects turbulent flows, are left as an area of future study.

Finally, the spatial location of RHESSI 35-50 keV X-rays, at the location of
minimum Doppler/non-thermal velocities within the flare cusp (high above the
thermal 6–12 keV source located at the flare looptop), indicate trapped high-
energy electrons above the flare looptop. This is another key prediction of a
flare termination shock. Together, these observational constraints provide
evidence for the presence of a termination shock within the September 10th
2017 solar flare.

Termination shocks can occur as long as super-magnetosonic reconnection
outflows are present during the flare energy release processes. Given the
magnitude and longevity of the September 10th 2017 flare, it is not
necessarily surprising to see evidence of a termination shock this late into
the flare’s evolution. Previous works on this event have found other examples
of high energy flare behavior several hours into the flare, including
prolonged high plasma-sheet velocities/temperatures and growth of the hot
loop-top cusp (French et al., 2020), prolonged plasma sheet downflows and
microwave emission (Yu et al., 2020), and extended (12+ hours) gamma ray
emission (Omodei et al., 2018).

To our knowledge, the results in this paper present the first reported
evidence of reconnection outflows driving a termination shock using direct
Doppler shift measurements.

## Data Availability

Hinode/EIS, RHESSI and AIA data are publicly available to all. Data cubes of
the two featured simulation frames will be provided via reasonable request to
the corresponding author.

## Acknowledgements

Hinode is a Japanese mission developed and launched by ISAS/JAXA, with NAOJ as
domestic partner and NASA and UKSA as international partners. It is operated
by these agencies in co-operation with ESA and NSC (Norway). R.F. thanks
support from the Brinson Prize Fellowship. S.J. and B.C. are supported by US
National Science Foundation grant AGS-2108853 and NASA grant 80NSSC20K1318 to
New Jersey Institute of Technology. C.S. is supported by NSF grant AST-2108438
and NASA grants 80NSSC21K2044, 80NSSC20K1318 to the Smithsonian Astrophysical
Observatory. S.A.M. is supported by UKSA grant (ST/X002063/1) and STFC grant
(ST/W001004/1).

## References

  * Carmichael (1964) Carmichael H., 1964, NASA Special Publication, 50, 451
  * Chen et al. (2015) Chen B., Bastian T. S., Shen C., Gary D. E., Krucker S., Glesener L., 2015, Science, 350, 1238
  * Chen et al. (2019) Chen B., Shen C., Reeves K. K., Guo F., Yu S., 2019, ApJ, 884, 63
  * Chen et al. (2020a) Chen B., et al., 2020a, Nature Astronomy, 4, 1140
  * Chen et al. (2020b) Chen B., Yu S., Reeves K. K., Gary D. E., 2020b, ApJ, 895, L50
  * Cheng et al. (2018) Cheng X., Li Y., Wan L. F., Ding M. D., Chen P. F., Zhang J., Liu J. J., 2018, ApJ, 866, 64
  * Culhane et al. (2007) Culhane J. L., et al., 2007, Sol. Phys., 243, 19
  * De Pontieu et al. (2014) De Pontieu B., et al., 2014, Sol. Phys., 289, 2733
  * Fleishman et al. (2020) Fleishman G. D., Gary D. E., Chen B., Kuroda N., Yu S., Nita G. M., 2020, Science, 367, 278
  * Forbes (1986) Forbes T. G., 1986, ApJ, 305, 553
  * French (2022) French R., 2022, Astronomy and Geophysics, 63, 5.38
  * French et al. (2019) French R. J., Judge P. G., Matthews S. A., van Driel-Gesztelyi L., 2019, ApJ, 887, L34
  * French et al. (2020) French R. J., Matthews S. A., van Driel-Gesztelyi L., Long D. M., Judge P. G., 2020, ApJ, 900, 192
  * Gary et al. (2018) Gary D. E., et al., 2018, ApJ, 863, 83
  * Hirayama (1974) Hirayama T., 1974, Sol. Phys., 34, 323
  * Kopp & Pneuman (1976) Kopp R. A., Pneuman G. W., 1976, Sol. Phys., 50, 85
  * Lemen et al. (2012) Lemen J. R., et al., 2012, Sol. Phys., 275, 17
  * Lin et al. (2002) Lin R. P., et al., 2002, Sol. Phys., 210, 3
  * Longcope et al. (2018) Longcope D., Unverferth J., Klein C., McCarthy M., Priest E., 2018, ApJ, 868, 148
  * Luo et al. (2021) Luo Y., Chen B., Yu S., Bastian T. S., Krucker S., 2021, ApJ, 911, 4
  * Masuda et al. (1994) Masuda S., Kosugi T., Hara H., Tsuneta S., Ogawara Y., 1994, Nature, 371, 495
  * Miles (2009) Miles A. R., 2009, ApJ, 696, 498
  * Omodei et al. (2018) Omodei N., Pesce-Rollins M., Longo F., Allafort A., Krucker S., 2018, ApJ, 865, L7
  * Polito et al. (2018) Polito V., Galan G., Reeves K. K., Musset S., 2018, ApJ, 865, 161
  * Reeves et al. (2020) Reeves K. K., Polito V., Chen B., Galan G., Yu S., Liu W., Li G., 2020, ApJ, 905, 165
  * Shen et al. (2018) Shen C., Kong X., Guo F., Raymond J. C., Chen B., 2018, ApJ, 869, 116
  * Shen et al. (2022) Shen C., Chen B., Reeves K. K., Yu S., Polito V., Xie X., 2022, Nature Astronomy, 6, 317
  * Shen et al. (2023) Shen C., Polito V., Reeves K. K., Chen B., Yu S., Xie X., 2023, Frontiers in Astronomy and Space Sciences, 10, 19
  * Shibata et al. (2023) Shibata K., Takasao S., Reeves K. K., 2023, ApJ, 943, 106
  * Sturrock (1968) Sturrock P. A., 1968, in Structure and Development of Solar Active Regions. p. 471
  * Takasao et al. (2015) Takasao S., Matsumoto T., Nakamura N., Shibata K., 2015, ApJ, 805, 135
  * Warren et al. (2018) Warren H. P., Brooks D. H., Ugarte-Urra I., Reep J. W., Crump N. A., Doschek G. A., 2018, ApJ, 854, 122
  * Yu et al. (2020) Yu S., Chen B., Reeves K. K., Gary D. E., Musset S., Fleishman G. D., Nita G. M., Glesener L., 2020, ApJ, 900, 17