Self-similar sequence transformation for critical exponents

V.I. Yukalov1,2 and E.P. Yukalova3

1Bogolubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, Dubna 141980, Russia

2Instituto de Fisica de São Carlos, Universidade de São Paulo,
CP 369, São Carlos 13560-970, São Paulo, Brazil

3Laboratory of Information Technologies,
Joint Institute for Nuclear Research, Dubna 141980, Russia

E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

Self-similar sequence transformation is an original type of nonlinear sequence
transformations allowing for defining effective limits of asymptotic
sequences. The method of self-similar factor transformations is shown to be
regular. This method is applied for calculating the critical exponents of the
$O(N)$-symmetric $\varphi^{4}$ theory in three dimensions by summing
asymptotic $\varepsilon$ expansions. It is shown that this method is
straightforward and essentially simpler than other summation techniques
involving complicated numerical calculations, while enjoying comparable
accuracy.

Keywords: Self-similar sequence transformation; asymptotic series; summation
methods; critical exponents

##  1 Introduction

Asymptotic expansions in powers of some parameter are widely used in physics
and applied mathematics [1, 2]. Since the parameters of interest very rarely
are really small, one needs to employ some kind of effective summation of
divergent series. The most popular are the method of Padé approximants [3],
Borel summation [4] and its variants, such as Padé-Borel summation and Borel
summation with conformal mapping. One also uses the methods of renormalization
group, conformal bootstrap, Monte Carlo simulations, and other methods (see
review [5]) requiring quite heavy numerical calculations.

In the present paper, we advocate another approach based on self-similar
approximation theory [6, 7, 8, 9]. The idea of this theory is to represent the
transition from one approximation term to another as the motion of a dynamical
system, or a renormalization-group equation, with the approximation order
playing the role of discrete time. Then the sequence of approximation terms
becomes bijective to the dynamical-system trajectory and the limit of the
sequence is bijective to the fixed point of the trajectory. In the vicinity of
a fixed point, the evolution equation (renormalization-group equation)
acquires the form of a self-similar relation, which explains the name of the
self-similar approximation theory. Mathematical details can be found in the
review articles [10, 11].

The approach has been applied to numerous problems, providing good agreement
with the exact results, when these are available, as well as with numerical
calculations and other elaborate methods, thus being compatible with other
methods of summation of divergent series. Accurate results can be obtained
even when just a few terms of an expansion are given, and when other summation
methods are not applicable at all (see reviews [10, 11]).

The method demonstrates good numerical convergence, which becomes especially
evident for the cases, where a large number (of order or larger than ten) of
perturbative terms are available. This concerns, e.g., the so-called zero-
dimensional model [12], one-dimensional anharmonic oscillator [12], and spin
glass [13].

All one needs for the application of the method is an asymptotic expansion
centered at a point on the real axis. Its analytical behavior on the whole
complex plane is not required. Thus, in the case of the zero-dimensional model
and the one-dimensional oscillator, we meet expressions that have on the
complex plane a singular point at zero [14]. This does not hinder the use of
the method resulting in numerically convergent sequence of approximants for
these models [12].

The method allows for the summation of a large class of functions, rational,
irrational and transcendental. Moreover, there exists a class of functions
that are exactly reproducible by this method [12]. This is the class of
functions having the form

$F_{k_{M}}(x)=\prod_{i=1}^{M}P_{m_{i}}^{\alpha_{i}}(x)$

of the product of polynomials

$P_{m_{i}}(x)=c_{i0}+c_{i1}x+c_{i2}x^{2}+\ldots+c_{i_{m_{i}}}x^{m_{i}}\,$

where $\alpha_{i}$ and $c_{ij}$ are either real or the powers $\alpha_{i}$ and
coefficients $c_{ij}$ are complex-valued numbers entering in complex conjugate
pairs so that $F_{k_{M}}$ is real, and

$k_{M}=\sum_{i=1}^{M}m_{i}+M\;.$

The exponential function also is shown to be reconstructed exactly starting
from the second-order approximation [12, 15].

Nonlinear differential equations, including singular equations, can be solved
by this method, first, by deriving a solution in terms of a series in powers
of a variable, and then by summing this series using the self-similar factor
transformation. For some nonlinear equations, exact soliton solutions have
been obtained [12, 16].

In the present paper, we consider two important points, one technical and the
other of great interest in physics. The first point is the proof of the method
regularity. We consider the method of self-similar factor transformation and
show that this method is regular, which implies that it sums every convergent
series to the same sum as that to which the series converges. This point is of
high importance for the justification of the approach. It is the standard way
of dealing with summation methods, when one, first, shows the regularity of
the method and then extrapolates it to divergent sequences, by demonstrating
its compatibility with other reliable methods and observing numerical
convergence for some test problems [17, 18]. Then we apply the method to the
transformation of $\varepsilon$ expansions for critical exponents. As a
concrete example, we study the $O(N)$ -symmetric $\varphi^{4}$ theory in three
dimensions. We show that the method of self-similar factor transformation
provides the accuracy comparable with other elaborate techniques involving
heavy numerical calculations, while being essentially simpler.

In Sec. 2, we formulate the method of self-similar factor transformation. We
do not plunge into the foundations of the whole theory, but will just give the
receipt of the method usage. In Sec. 3, we prove the regularity of the method.
The application to defining the critical exponents is given in Sec. 4, where
we compare our results with the most accurate values obtained by other methods
summarized in Refs. [5, 19, 20, 21]. We compare our results with Monte Carlo
simulations, conformal bootstrap, hypergeometric Meijer summation, Borel
summation, Borel summation with conformal mapping, and the method of
nonperturbative renormalization group. This comparison shows good agreement of
the self-similar factor approximants with the calculations by other methods.
The last Sec. 5 concludes.

##  2 Self-similar factor transformation

In this section, we describe the method that can be used for the summation of
arbitrary asymptotic series. Suppose we have got an asymptotic expansion for a
real function

$f_{k}(x)=f_{0}(x)\left(1+\sum_{n=1}^{k}a_{n}x^{n}\right)$ (1)

in powers of a real parameter $x$ assumed to be asymptotically small. However,
we need to find the value of the function at a finite value of the parameter.
The extrapolation of the expansion $f_{k}(x)$ to arbitrary values of the
parameter $x$ can be done by means of the self-similar factor transformation
[10, 22, 23]. This method transforms the truncated expansion (1) to the factor
form

$f_{k}^{*}(x)=f_{0}(x)\prod_{j=1}^{N_{k}}(1+A_{j}x)^{n_{j}}\;,$ (2)

where the number of factors is

$\displaystyle N_{k}=\left\\{\begin{array}[]{rl}k/2,{}&~{}k=2,4,6,\ldots\\\
(k+1)/2,{}&~{}k=3,5,7,\ldots\end{array}\right..$ (5)

The parameters $A_{j}$ and $n_{j}$ are uniquely determined by the accuracy-
through-order procedure, by equating the like-order terms in the expansions at
small $x$,

$f_{k}^{*}(x)\simeq f_{k}(x)\qquad(x\rightarrow 0)\;.$ (6)

This procedure gives the equations

$\sum_{j=1}^{N_{k}}n_{j}A_{j}^{n}=D_{n}\qquad(n=1,2,\ldots,k)\;,$ (7)

in which

$D_{n}\equiv\frac{(-1)^{n-1}}{(n-1)!}\;\lim_{x\rightarrow
0}\;\frac{d^{n}}{dx^{n}}\;\ln\left(1+\sum_{m=1}^{n}a_{m}x^{m}\right)\;.$ (8)

In the case of an even order $k$, Eq. (7) consists of $k$ equations uniquely
defining the $k/2$ parameters $A_{j}$ and $k/2$ parameters $n_{j}$. For odd
orders $k$, the system of $k$ equations (7) contains $k+1$ unknowns, where one
of the parameters $A_{j}$, say $A_{1}$, is arbitrary. Normalizing $A_{j}$ in
units of $A_{1}$ implies $A_{1}=1$, which makes the system of equations (7)
self-consistent and all parameters uniquely defined [10, 15]. If the found
parameters lead to a complex-valued approximant, it is replaced by the nearest
real-valued approximant. The final result is given by the average between the
last two approximants $[f_{k}^{*}(x)+f_{k-1}^{*}(x)]/2$ and the error bar is
defined as the half-difference between the last two different approximants
$[f_{k}^{*}(x)-f_{k-1}^{*}(x)]/2$. As is seen, the scheme is very simple and
straightforward.

##  3 Method regularity

In this section, we show that the self-similar factor transformation is a
regular method of summation.

Theorem. Let us consider the sequence of the terms

$f_{k}(x)=\sum_{n=0}^{k}a_{n}x^{n}\qquad(x\in\mathbb{D})$ (9)

that are defined on a domain $\mathbb{D}\subset\mathbb{R}$ including the point
$x=0$. The self-similar factor transformation reduces the term $f_{k}(x)$ to
the form

$f_{k}^{*}(x)\equiv a_{0}\prod_{j=1}^{N_{k}}(1+A_{j}x)^{n_{j}}$ (10)

that is a smooth function on $\mathbb{D}$, hence being infinitely
differentiable, whose parameters $A_{j}$ and $n_{j}$ are prescribed by the
accuracy-through-order procedure, such that

$\frac{f_{k}^{*(n)}(0)}{n!}=a_{n}\qquad(n=0,1,2,\ldots)\;,$ (11)

where $f_{k}^{*(n)}(x)$ is the $n$-th derivative of $f_{k}^{*}(x)$ over $x$.

If the sequence of terms (9) converges to a smooth function

$f(x)=\sum_{n=0}^{\infty}a_{n}x^{n}\;,$ (12)

then the sequence of the factor approximants (10) converges to the function

$f^{*}(x)\equiv a_{0}\prod_{j=1}^{\infty}(1+A_{j}x)^{n_{j}}$ (13)

coinciding with $f(x)$:

$f^{*}(x)=f(x)\;.$

Proof. The convergence of the sequence of terms (9) to a smooth function
$f(x)$ implies that the latter can be represented as the Taylor series

$f(x)=\sum_{n=0}^{\infty}\frac{f^{(n)}(0)}{n!}\;x^{n}\;,$ (14)

with

$\frac{f^{(n)}(0)}{n!}=a_{n}\qquad(n=0,1,2,\ldots)$ (15)

and with the remainder

$R_{k}[\;f(x)\;]=\frac{f^{(k+1)}(c_{1})}{(k+1)!}\;x^{k+1}\;,$ (16)

where $c_{1}\in(0,x)$, this remainder tending to zero when
$k\rightarrow\infty$,

$\lim_{k\rightarrow\infty}\;R_{k}[\;f(x)\;]=0\;.$ (17)

Comparing expressions (11) and (15), we see that

$f^{(n)}(0)-f_{n}^{*(n)}(0)=0\qquad(n=0,1,2,\ldots)\;.$ (18)

Because of this, and since the derivatives are continuous, there exists a
finite value $\varepsilon$ such that

$|\;f^{(n)}(c_{1})-f_{n}^{*(n)}(c_{2})\;|<\varepsilon<\infty\;,$ (19)

where $0\leq c_{2}\leq x$ and $n=0,1,2,\ldots$. Considering the difference
between the remainder (16) and the remainder

$R_{k}[\;f_{k}^{*}(x)\;]=\frac{f_{k}^{*(k+1)}(c_{2})}{(k+1)!}\;x^{k+1}\;,$
(20)

we have

$|\;R_{k}[\;f_{k}^{*}(x)\;]-R_{k}[\;f(x)\;]\;|<\frac{\varepsilon}{(k+1)!}\;|\;x\;|^{k+1}\;.$
(21)

Using the Stirling formula $n!\simeq\sqrt{2\pi n}(n/e)^{n}$, we find that

$\lim_{n\rightarrow\infty}\;\frac{\varepsilon}{(n+1)!}\;|\;x\;|^{n+1}=\frac{\varepsilon}{\sqrt{2\pi}}\;\lim_{n\rightarrow\infty}\;\left(\frac{|\;x\;|\;e}{n}\right)^{n}=0$

for any fixed $x$. From here, and taking into account the limit (17), we get

$\lim_{k\rightarrow\infty}\;R_{k}[\;f_{k}^{*}(x)\;]=0\;.$ (22)

Then, by the Taylor theorem, the factor approximant (13) can be represented in
the form of the Taylor series

$f^{*}(x)=\sum_{n=0}^{\infty}\frac{f^{*(n)}(0)}{n!}\;x^{n}\;,$ (23)

with

$\frac{f^{*(n)}(0)}{n!}=a_{n}\qquad(n=0,1,2,\ldots)\;.$ (24)

Comparing series (14) and (23), under equalities (15) and (24), we come to the
conclusion that the functions $f^{*}(x)$ and $f(x)$ coincide. $\square$

##  4 Critical exponents

Critical exponents can be presented in the form of $\varepsilon$-expansions in
powers of $\varepsilon=4-d$, where $d$ is space dimensionality. For
$O(N)$-symmetric $\varphi^{4}$ theory in three dimensions, the five-loop
expansions can be found in the book [24]. The summation of the five-loop
expansions by means of self-similar approximants was considered in Ref. [13]
for all $N$. It was shown that for $N=-2$ and $N\rightarrow\infty$ self-
similar approximations yield the exact values for the exponents. The results
are very accurate for large $N\gg 1$, with the errors decreasing as $1/N$ with
increasing $N$. However the accuracy for the lower $N$ was not sufficient.

Our aim in the present paper is to demonstrate that the accuracy of self-
similar factor approximants for the $O(N)$ -symmetric $\varphi^{4}$ theory in
three dimensions can be drastically improved by employing the available seven-
loop $\varepsilon$-expansions that are known for $N=1$ [25] and are derived,
using the seven-loop coupling parameter expansions [26], in Ref. [27] for
$N=0,1,2,3,4$. These expansions are as follows.

(i) $N=0$.

For $N=0$, we have

$\nu^{-1}=2-0.25\;\varepsilon-0.08594\;\varepsilon^{2}+0.11443\;\varepsilon^{3}-0.28751\;\varepsilon^{4}+0.95613\;\varepsilon^{5}-3.8558\;\varepsilon^{6}+17.784\;\varepsilon^{7}\;,$
(25)
$\eta=0.015625\;\varepsilon^{2}+0.016602\;\varepsilon^{3}-0.0083675\;\varepsilon^{4}+0.026505\;\varepsilon^{5}-0.09073\;\varepsilon^{6}+0.37851\;\varepsilon^{7}\;,$
(26)

and

$\omega=\varepsilon-0.65625\;\varepsilon^{2}+1.8236\;\varepsilon^{3}-6.2854\;\varepsilon^{4}+26.873\;\varepsilon^{5}-130.01\;\varepsilon^{6}+692.1\;\varepsilon^{7}\;.$
(27)

We calculate the corresponding self-similar factor approximants
$f_{k}^{*}(\varepsilon)$, as is explained in the previous section, and set
$\varepsilon=1$. The results of this calculation for the exponent $\nu$ are
illustrated in Table 1. Following the same procedure for the exponents $\eta$
and $\omega$, we obtain the values shown in Table 2, where we compare the
results obtained by means of factor approximants (FA) with those of other
methods: Monte Carlo simulations (MC) [28, 29, 30, 31, 32, 33], Conformal
bootstrap (CB) [34, 35, 36, 37, 38], Hypergeometric Meijer summation (HGM)
[27], Borel summation complimented by additional conjectures on the behavior
of coefficients (BAC) [39], Borel summation with conformal mapping (BCM) [40],
and Nonperturbative renormalization group (NPRG) [5, 41, 42, 43].

(ii) $N=1$.

The $\varepsilon$ expansions for the critical exponents read as

$\nu^{-1}=2-0.333333\;\varepsilon-0.11728\;\varepsilon^{2}+0.12453\;\varepsilon^{3}-0.30685\;\varepsilon^{4}+0.95124\;\varepsilon^{5}-3.5726\;\varepsilon^{6}+15.287\;\varepsilon^{7}\;,$
(28)
$\eta=0.018519\;\varepsilon^{2}+0.01869\;\varepsilon^{3}-0.0083288\;\varepsilon^{4}+0.025656\;\varepsilon^{5}-0.081273\;\varepsilon^{6}+0.31475\;\varepsilon^{7}\;,$
(29)

and

$\omega=\varepsilon-0.62963\;\varepsilon^{2}+1.6182\;\varepsilon^{3}-5.2351\;\varepsilon^{4}+20.75\;\varepsilon^{5}-93.111\;\varepsilon^{6}+458.74\;\varepsilon^{7}\;.$
(30)

The results for the factor approximants are summarized in Table 3, where they
are compared with the values obtained by other methods listed above.

(iii) $N=2$.

The $\varepsilon$ expansions are

$\nu^{-1}=2-0.4\;\varepsilon-0.14\;\varepsilon^{2}+0.12244\;\varepsilon^{3}-0.30473\;\varepsilon^{4}+0.87924\;\varepsilon^{5}-3.103\;\varepsilon^{6}+12.419\;\varepsilon^{7}\;,$
(31)
$\eta=0.02\;\varepsilon^{2}+0.019\;\varepsilon^{3}-0.0078936\;\varepsilon^{4}+0.023209\;\varepsilon^{5}-0.068627\;\varepsilon^{6}+0.24861\;\varepsilon^{7}\;,$
(32)

and

$\omega=\varepsilon-0.6\;\varepsilon^{2}+1.4372\;\varepsilon^{3}-4.4203\;\varepsilon^{4}+16.374\;\varepsilon^{5}-68.777\;\varepsilon^{6}+316.48\;\varepsilon^{7}\;.$
(33)

The calculated factor approximants for the critical exponents are presented in
Table 4, where they are compared with the exponents found by other methods
listed above.

(iv) $N=3$.

The $\varepsilon$ expansions read as

$\nu^{-1}=2-0.45455\;\varepsilon-0.1559\;\varepsilon^{2}+0.11507\;\varepsilon^{3}-0.2936\;\varepsilon^{4}+0.78994\;\varepsilon^{5}-2.6392\;\varepsilon^{6}+9.9452\;\varepsilon^{7}\;,$
(34)
$\eta=0.020661\;\varepsilon^{2}+0.018399\;\varepsilon^{3}-0.0074495\;\varepsilon^{4}+0.020383\;\varepsilon^{5}-0.057024\;\varepsilon^{6}+0.19422\;\varepsilon^{7}\;,$
(35)

and

$\omega=\varepsilon-0.57025\;\varepsilon^{2}+1.2829\;\varepsilon^{3}-3.7811\;\varepsilon^{4}+13.182\;\varepsilon^{5}-52.204\;\varepsilon^{6}+226.02\;\varepsilon^{7}\;.$
(36)

The corresponding exponents calculated by means of the factor approximants are
shown in Table 5, together with the results of calculation by other methods
listed above.

(v) $N=4$.

The $\varepsilon$ expansions take the form

$\nu^{-1}=2-0.5\;\varepsilon-0.16667\;\varepsilon^{2}+0.10586\;\varepsilon^{3}-0.27866\;\varepsilon^{4}+0.70217\;\varepsilon^{5}-2.2337\;\varepsilon^{6}+7.9701\;\varepsilon^{7}\;,$
(37)
$\eta=0.020833\;\varepsilon^{2}+0.017361\;\varepsilon^{3}-0.0070852\;\varepsilon^{4}+0.017631\;\varepsilon^{5}-0.047363\;\varepsilon^{6}+0.15219\;\varepsilon^{7}\;,$
(38)

and

$\omega=\varepsilon-0.54167\;\varepsilon^{2}+1.1526\;\varepsilon^{3}-3.2719\;\varepsilon^{4}+10.802\;\varepsilon^{5}-40.567\;\varepsilon^{6}+166.26\;\varepsilon^{7}\;.$
(39)

The resulting values of the factor approximants and the values of the
exponents found by other methods listed above are given in Table 6.

As has been explained above, the construction of the factor approximants for
each given expansion is straightforward. In order that the reader would grasp
the overall structure of these approximants, we adduce below, as an example,
the explicit expressions for $\nu^{-1}$ in the case of $N=0$. The self-similar
approximant of order $k=2$ is

$f_{2}^{*}(\varepsilon)=2(1-0.81252\varepsilon)^{0.153842}\;.$

For the approximant of order $k=3$, we have

$f_{3}^{*}(\varepsilon)=\frac{2(1+\varepsilon)^{0.169077}}{(1+0.229573\;\varepsilon)^{1.28098}}\;.$

In the fourth order $(k=4)$, we get

$f_{4}^{*}(\varepsilon)=2(1-0.343467\;\varepsilon)^{0.411916}(1+3.21435\;\varepsilon)^{0.0051269}\;.$

The fifth order $(k=5)$ yields

$f_{5}^{*}(\varepsilon)=2(1+\varepsilon)^{0.0459476}(1-0.172432\;\varepsilon)^{1.02395}(1+4.48424\;\varepsilon)^{0.00125184}\;.$

The sixth order $k=6$ gives

$f_{6}^{*}(\varepsilon)=2(1-0.261332\;\varepsilon)^{0.58267}(1+1.98816\;\varepsilon)^{0.0126183}(1+5.44898\;\varepsilon)^{0.000400602}\;.$

And the approximant of seventh order $(k=7)$ is given by the expression

$f_{7}^{*}(\varepsilon)=2(1+\varepsilon)^{0.0275211}(1-0.222915\;\varepsilon)^{0.734195}(1+3.13825\;\varepsilon)^{0.0032879}(1+6.2879\;\varepsilon)^{0.000131004}\;.$

Substituting here $\varepsilon=1$, we come to the corresponding values of the
exponent $\nu$.

All other factor approximants for the related expansions are obtained in the
same way. In addition to the exponents $\nu$, $\eta$, and $\omega$, we find
other exponents $\alpha$, $\beta$, $\gamma$, and $\delta$ through the
relations

$\alpha=2-3\nu\;,\qquad\beta=\frac{\nu}{2}\;(1+\eta)\;\,\qquad\gamma=\nu(2-\eta)\;,\qquad\delta=\frac{5-\eta}{1+\eta}\;.$
(40)

Table 7 summarizes these results.

As is seen, the self-similar factor approximants are in good agreement with
the results of other methods, while the approach using factor approximants is
very simple and allowing for the explicit analytical construction of these
approximants.

It is worth mentioning one point, where the situation remains not completely
settled. This is the case of the exponent $\alpha$ for $N=2$. The case of
$N=2$ is of special interest representing the superfluid helium and magnetic
$XY$ models [44, 45]. The situation concerns the fact that practically all
calculational methods yield the values of $\alpha$ that are in good accordance
with each other as well as with the experimental values, but that are a bit
lower than the result of Monte Carlo simulations and the value following from
the conformal bootstrap conjecture.

Table 8 illustrates the situation for $N=2$, comparing the exponent $\alpha$
obtained by different methods listed in Table 2 with experimental data,
including the extremely precise results of the specific heat measurements for
liquid helium in zero gravity [46] and the series of measurements for several
magnetic materials of the $XY$ class [47, 48, 49].

##  5 Discussion

The method of self-similar transformations is a very simple and convenient
tool for the summation of asymptotic series. The basis of the method is the
consideration of the transfer from one approximation to another as a motion in
the space of approximations, with the approximation order playing the role of
time. The motion in the vicinity of a fixed point is described by an equation
having the form of a self-similar relation, which is equivalent to a
renormalization-group equation. The fixed point of the evolution equation
defines the sought effective limit of the transformed sequence. A
representation for the effective limit acquires the form of self-similar
factor approximants. Following the usual way of dealing with sequence
transformations, we show that the method is regular and then extrapolate it to
the case of divergent sequences.

We apply the method of self-similar factor transformations for the summation
of $\varepsilon$ expansions for the $O(N)$\- symmetric theory in three
dimensions. The series of seventh order in $\varepsilon$ are used. The method
is shown to provide accurate approximations, at the same time being very
simple and allowing for the construction of explicit analytical expressions.
The results are compatible with other known methods of summation.

Employing this method, it is even possible to get reasonable estimates for
two-dimensional systems, for which one has to set $\varepsilon=2$. For the
two-dimensional systems, the symmetry can be broken only for $N=0$ and $N=1$
[50]. The estimates for the corresponding critical exponents can be compared
with the values conjectured in Ref. [51] for $N=0$ and with the known exact
values [52] for $N=1$. Thus the self-similar factor approximants for $N=0$
give $\nu=0.748$, $\eta=0.18$, and $\omega=1.62$, as compared with the exact
values $\nu=0.75$, $\eta=5/24$, and $\omega=2$. And for $N=1$, we obtain
$\nu=0.997$, $\eta=0.21$, and $\omega=1.6$, as compared with the exact values
$\nu=1$, $\eta=0.25$, and $\omega=1.75$.

Author Contributions

Both the authors, V.I. Yukalov and E.P. Yukalova, equally contributed to this
work.

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

## References

  * [1] G.E. Giacaglia, Perturbation Methods in Non-Linear Systems, Springer, New York, 1972.
  * [2] A.H. Nayfeh, Perturbation Methods, Wiley, New York, 1973.
  * [3] G.A. Baker, P. Graves-Moris, Padé Approximants, Cambridge University, Cambridge, 1996.
  * [4] H. Kleinert, Path Integrals, World Scientific, Singapore, 2004.
  * [5] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J.M. Pawlowski, M. Tissier, N. Wschebor, Phys. Rep. 910 (2021) 1.
  * [6] V.I. Yukalov, Phys. Rev. A 42 (1990) 3324.
  * [7] V.I. Yukalov, Physica A 167 (1990) 833.
  * [8] V.I. Yukalov, J. Math. Phys. 32 (1991) 1235.
  * [9] V.I. Yukalov, J. Math. Phys. 33 (1992) 3994.
  * [10] V.I. Yukalov, Phys. Part. Nucl. 50 (2019) 141.
  * [11] V.I. Yukalov, E.P. Yukalova, Physics 3 (2021) 829.
  * [12] V.I. Yukalov, E.P. Yukalova, Phys. Rev. D 103 (2021) 076019.
  * [13] V.I. Yukalov, E.P. Yukalova, Eur. Phys. J. B 55 (2007) 93.
  * [14] C.M. Bender, T.T. Wu, Phys. Rev. 184 (1969) 1231.
  * [15] V.I. Yukalov, E.P. Yukalova, Phys. Lett. A 368 (2007) 341.
  * [16] E.P. Yukalova, V.I. Yukalov, S. Gluzman, Ann. Phys. (N.Y.) 323 (2008) 3074.
  * [17] G.H. Hardy, Divergent Series, Clarendon, Oxford, 1973.
  * [18] E.J. Weniger, Comput. Phys. Rep. 10 (1989) 189.
  * [19] A. Pelissetto, E. Vicari, Phys. Rep. 368 (2002) 549.
  * [20] V. Abhignan, R. Sankaranarayanan, J. Stat. Phys. 183 (2021) 4.
  * [21] A.M. Shalaby, Eur. Phys. J. C 81 (2021) 87.
  * [22] V.I. Yukalov, S. Gluzman, D. Sornette, Physica A 328 (2003) 409.
  * [23] S. Gluzman, V.I. Yukalov, D. Sornette, Phys. Rev. E 67 (2003) 026109.
  * [24] H. Kleinert, V. Schulte-Frohlinde, Critical Properties of $\varphi^{4}$-Theories, World Scientific, Singapore, 2006.
  * [25] T.A. Ryttov, J. High Energy Phys. 04 (2020) 072.
  * [26] O. Schnetz, Phys. Rev. D 97 (2018) 085018.
  * [27] A.M. Shalaby, Eur. Phys. J. C 81 (2021) 87.
  * [28] M. Hasenbusch, J. Phys. A 34 (2001) 8221.
  * [29] M. Hasenbusch, Phys. Rev. B 82 (2010) 174433.
  * [30] M. Hasenbusch, E. Vicari, Phys. Rev. B 84 (2011) 125136.
  * [31] N. Clisby, B. Dunweg, Phys. Rev. E 94 (2016) 052102.
  * [32] N. Clisby, J. Phys. A 50 (2017) 264003.
  * [33] M. Hasenbusch, Phys. Rev. B 100 (2019) 224517.
  * [34] S. El-Showk, M.F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin, A. Vichi, J. Stat. Phys. 157 (2014) 869.
  * [35] H. Shimada, S. Hikami, J. Stat. Phys. 165 (2016) 1006.
  * [36] F. Kos, D. Poland, D. Simmons-Duffin, A. Vichi, J. High Energy Phys. 08 (2016) 036.
  * [37] A.C. Echeverri, B. von Harling, M. Serone, J. High Energy Phys. 09 (2016) 097.
  * [38] D. Simmons-Duffin, J. High Energy Phys. 03 (2017) 086.
  * [39] M.V. Kompaniets, K.J. Wiese, Phys. Rev. E 101 (2020) 012104.
  * [40] M.V. Kompaniets, E. Panzer, Phys. Rev. D 96 (2017) 036016.
  * [41] N. Hasselmann, Phys. Rev. E 86 (2012) 041118.
  * [42] F. Rose, N. Dupuis, Phys. Rev. B 97 (2018) 174514.
  * [43] G. De Polsi, I. Balog, M. Tissier, N. Wschebor, Phys. Rev. E 101 (2020) 042113.
  * [44] H. Kleinert, V.I. Yukalov, Phys. Rev. E 71 (2005) 026131.
  * [45] A.M. Shalaby, Phys. Rev. D 102 (2020) 105017.
  * [46] J.A. Lipa, J.A. Nissen, D.A. Stricker, D.R. Swanson, T.C. Chui, Phys. Rev. B 68 (2003) 174518.
  * [47] A.N. Vasilev, The Field Theoretic Renormalization Group in Critical Behavior Theory and Stochastic Dynamics, Chapman and Hall, Boca Raton, 2004.
  * [48] A. Olega, A. Salazar, D. Prabhakaran, J.G. Cheng, J.S. Zhou, Phys. Rev. B 85 (2012) 184425.
  * [49] A. Olega, A. Salazar, Y.M. Bunkov, J. Phys. Condens. Matter 26 (2014) 096001.
  * [50] A. Pelissetto, E. Vicari, Nucl. Phys. B 519 (1998) 626.
  * [51] B. Nienhuis, Phys. Rev. Lett. 49 (1982) 1062.
  * [52] P. Calabrese, M. Caselle, A. Celi, A. Pelissetto, E. Vicari, J. Phys. A 33 (2000) 8155.

Table 1: Self-similar factor approximants of order $k=2,3,4,5,6,7$ for the exponent $\nu$ and different number of components $N$. $k$ | $N=0$ | $N=1$ | $N=2$ | $N=3$ | $N=4$
---|---|---|---|---|---
2 | 0.64688 | 0.73938 | 0.83405 | 0.91870 | 0.98471
3 | 0.57949 | 0.61666 | 0.65176 | 0.68457 | 0.71492
4 | 0.59026 | 0.63394 | 0.67562 | 0.71452 | 0.75006
5 | 0.58665 | 0.62845 | 0.66875 | 0.70688 | 0.74221
6 | 0.58789 | 0.63030 | 0.67127 | 0.71000 | 0.74573
7 | 0.58744 | 0.62968 | 0.67073 | 0.70985 | 0.74614
Table 2: Critical exponents for $N=0$, found by different methods: Self-similar factor approximants (FA), Monte Carlo simulations (MC), Conformal bootstrap (CB), Hypergeometric Meijer summation (HGM), Borel summation with additional conjectures on the behaviour of coefficients (BAC), Borel summation with conformal mapping (BCM), and Nonperturbative renormalization group (NPRG). $Method$ | $\nu$ | $\eta$ | $\omega$
---|---|---|---
FA | 0.5877 (2) | 0.0301 (2) | 0.821 (15)
MC | 0.5875970 (4) | 0.031043 (3) | 0.899 (12)
CB | 0.5877 (12) | 0.0282 (4) | $-$
HGM | 0.5877 (2) | 0.0312 (7) | 0.8484 (17)
BAC | 0.5874 (2) | 0.0304 (2) | 0.846 (15)
BCM | 0.5874 (3) | 0.0310 (7) | 0.841 (13)
NPRG | 0.5876 (2) | 0.0312 (9) | 0.901 (24)
Table 3: Critical exponents for $N=1$, found by different methods listed in Table 2. $Method$ | $\nu$ | $\eta$ | $\omega$
---|---|---|---
FA | 0.6300 (3) | 0.0353 (3) | 0.808 (9)
MC | 0.63002 (10) | 0.03627 (10) | 0.832 (6)
CB | 0.62999 (5) | 0.03631 (3) | 0.830 (2)
HGM | 0.6298 (2) | 0.0365 (7) | 0.8231 (5)
BAC | 0.6296 (3) | 0.0355 (3) | 0.827 (13)
BCM | 0.6292 (5) | 0.0362 (6) | 0.820 (7)
NPRG | 0.63012 (16) | 0.0361 (11) | 0.832 (14)
Table 4: Critical exponents for $N=2$, found by different methods listed in Table 2. $Method$ | $\nu$ | $\eta$ | $\omega$
---|---|---|---
FA | 0.6710 (3) | 0.0372 (4) | 0.809 (11)
MC | 0.67169 (7) | 0.03810 (8) | 0.789 (4)
CB | 0.67175 (10) | 0.0385 (6) | 0.811 (10)
HGM | 0.6708 (4) | 0.0381 (6) | 0.789 (13)
BAC | 0.6706 (2) | 0.0374 (3) | 0.808 (7)
BCM | 0.6690 (10) | 0.0380 (6) | 0.804 (3)
NPRG | 0.6716 (6) | 0.0380 (13) | 0.791 (8)
Table 5: Critical exponents for $N=3$, found by different methods listed in Table 2. $Method$ | $\nu$ | $\eta$ | $\omega$
---|---|---|---
FA | 0.7099 (1) | 0.0372 (4) | 0.7919 (3)
MC | 0.7116 (10) | 0.0378 (3) | 0.773
CB | 0.7121 (28) | 0.0386 (12) | 0.791 (22)
HGM | 0.7091 (2) | 0.0381 (6) | 0.764 (18)
BAC | 0.70944 (2) | 0.0373 (3) | 0.794 (4)
BCM | 0.7059 (20) | 0.0378 (5) | 0.795 (7)
NPRG | 0.7114 (9) | 0.0376 (13) | 0.796 (11)
Table 6: Critical exponents for $N=4$, found by different methods listed in Table 2. $Method$ | $\nu$ | $\eta$ | $\omega$
---|---|---|---
FA | 0.7459 (2) | 0.0361 (4) | 0.7913 (8)
MC | 0.750 (2) | 0.0360 (3) | 0.765 (30)
CB | 0.751 (3) | 0.0378 (32) | 0.817 (30)
HGM | 0.7443 (3) | 0.0367 (4) | 0.7519 (13)
BAC | 0.7449 (4) | 0.0363 (2) | 0.7863 (9)
BCM | 0.7397 (35) | 0.0366 (4) | 0.794 (9)
NPRG | 0.7478 (9) | 0.0360 (12) | 0.761 (12)
Table 7: Critical exponents for the three-dimensional $O(N)$-symmetric $\varphi^{4}$ field theory, calculated by means of self-similar factor approximants $N$ | $\alpha$ | $\beta$ | $\gamma$ | $\delta$ | $\nu$ | $\eta$ | $\omega$
---|---|---|---|---|---|---|---
0 | 0.2369 | 0.3027 | 1.15771 | 4.8247 | 0.5877 | 0.0301 | 0.821
1 | 0.1100 | 0.3261 | 1.23776 | 4.7954 | 0.6300 | 0.0353 | 0.808
2 |  $-$0.0130 | 0.3480 | 1.31704 | 4.7848 | 0.6710 | 0.0372 | 0.809
3 |  $-$0.1297 | 0.3682 | 1.39339 | 4.7848 | 0.7099 | 0.0372 | 0.792
4 |  $-$0.2377 | 0.3864 | 1.46487 | 4.7910 | 0.7459 | 0.0361 | 0.791
Table 8: Critical exponents $\alpha$ for $N=2$, found by different methods listed in Table 2, and experimental data. $Method$ | $\alpha$
---|---
FA |  $-$0.0130 (8)
MC |  $-$0.0151 (2)
CB |  $-$0.0152 (3)
HGM |  $-$0.0124 (12)
BAC |  $-$0.0118 (6)
BCM |  $-$0.0070 (30)
NPRG |  $-$0.0148 (18)
experiment [46] |  $-$0.0127 (3)
experiment [47] |  $-$0.0130 (30)
experiment [48] |  $-$0.0130 (10)
experiment [49] |  $-$0.0130 (20)


11institutetext: Georgia State University, Atlanta, GA, USA
11email: {cpandey1, rangryk<EMAIL_ADDRESS>
22institutetext: Research Center for Astronomy and Applied Mathematics,
Academy of Athens, Athens, Greece
22email<EMAIL_ADDRESS>
# Explainable Deep Learning-based Solar Flare Prediction with post hoc
Attention for Operational Forecasting

Chetraj Pandey✉ 11 0000-0002-4699-4050 Rafal A. Angryk 11 0000-0001-9598-8207
Manolis K. Georgoulis 22 0000-0001-6913-1330 Berkay Aydin 11
0000-0002-9799-9265

###### Abstract

This paper presents a post hoc analysis of a deep learning-based full-disk
solar flare prediction model. We used hourly full-disk line-of-sight
magnetogram images and selected binary prediction mode to predict the
occurrence of $\geq$M1.0-class flares within 24 hours. We leveraged custom
data augmentation and sample weighting to counter the inherent class-imbalance
problem and used true skill statistic and Heidke skill score as evaluation
metrics. Recent advancements in gradient-based attention methods allow us to
interpret models by sending gradient signals to assign the burden of the
decision on the input features. We interpret our model using three post hoc
attention methods: (i) Guided Gradient-weighted Class Activation Mapping, (ii)
Deep Shapley Additive Explanations, and (iii) Integrated Gradients. Our
analysis shows that full-disk predictions of solar flares align with
characteristics related to the active regions. The key findings of this study
are: (1) We demonstrate that our full disk model can tangibly locate and
predict near-limb solar flares, which is a critical feature for operational
flare forecasting, (2) Our candidate model achieves an average
TSS=0.51$\pm$0.05 and HSS=0.38$\pm$0.08, and (3) Our evaluation suggests that
these models can learn conspicuous features corresponding to active regions
from full-disk magnetograms.

###### Keywords:

Solar flares Deep learning xAI Interpretability

##  1 Introduction

Solar flares are transient solar events of central importance to space weather
forecasting, manifested as the sudden large eruption of electromagnetic
radiation on the outermost atmosphere of the Sun. They are classified
according to their peak X-ray flux level into the following five categories by
National Oceanic and Atmospheric Administration (NOAA): X $(\geq
10^{-4}Wm^{-2})$, M $(\geq 10^{-5}Wm^{-2})$, C $(\geq 10^{-6}Wm^{-2})$, B
$(\geq 10^{-7}Wm^{-2})$, and A $(\geq 10^{-8}Wm^{-2})$, where,
X$>$M$>$C$>$B$>$A [6]. These flare classes are on a logarithmic scale, meaning
that each class represents a tenfold increase in X-ray flux compared to the
previous class. Large flares (M- and X-class) are scarce events that are more
likely to incur a terrestrial impact and, therefore, the classes of interest
that gather the attention of researchers. These flares may potentially disrupt
the electricity supply chain, airline industry, and satellite communications,
and pose radiation hazards to astronauts in space. To mitigate these risks,
the necessity of a precise and reliable flare prediction model becomes
imperative.

Active regions (ARs) on the Sun are places characterized by the largest
accumulations of dipolar magnetic flux in the solar atmosphere. Most
operational flare forecasts target these regions of interest and issue
predictions for individual ARs, which are the main initiators of space weather
events. To issue a full-disk forecast with an AR-based model, the output flare
probabilities for each active region are usually aggregated using a heuristic
function as mentioned in [20]. The heuristic function used to aggregate the
final forecast operates under the assumption of conditional independence among
ARs and that all ARs contribute equally to the aggregate forecast. This
uniform weighting scheme may not accurately reflect the true influence of each
AR on full-disk flare prediction probability. It is important to highlight
that the weights of these ARs are generally unknown; there are no established
methods to accurately determine them, nor are there any prior assumptions that
guide the assignment of these weights.

Furthermore, the magnetic field measurements, employed by the AR-based
forecasting techniques, are susceptible to severe projection effects as ARs
get closer to limbs (to the degree that after $\pm$60∘ the magnetic field
readings are distorted [5]); therefore, the aggregated full-disk flare
probability is in fact, restrictive (i.e., from ARs in central locations) as
the data in itself is limited to ARs located within $\pm$45∘ [11] to $\pm$70∘
[9] and in some cases, even $\pm$30∘ [8] due to severe projection effects [7].
As AR-based models include data up to $\pm$70∘, in the context of this paper,
this upper limit ($\pm$70∘) is used as a boundary for central (within
$\pm$70∘) and near-limb regions (beyond $\pm$70∘).

In contrast to AR-based models, which use individual AR data from central
locations, full-disk models use complete magnetogram images corresponding to
the entire disk. These images are typically compressed JP2 (JPEG 2000) 8-bit
representations (i.e., pixel values ranging from 0 to 255) derived from
original magnetogram rasters which contain magnetic field strength values
ranging from $\sim$$\pm$4500G. The compressed magnetogram images are used for
shape-based parameters, e.g., size, directionality, borders, and inversion
lines. Although projection effects still prevail in these images, full-disk
models can learn from the near-limb areas. Thus, incorporating a full-disk
model is essential to supplement AR-based models, enabling the prediction of
flares in the Sun’s near-limb areas and enhancing operational flare
forecasting systems.

With recent advancements in machine learning and deep learning methods, their
application in predicting solar flares has demonstrated great experimental
success and accelerated the efforts of many interdisciplinary researchers [8,
11, 15, 16, 17, 18, 19, 32]. Although deep learning methods have significantly
enhanced solutions to image classification and computer vision problems, these
models learn highly complex data representations, rendering them as black-box
models. Consequently, the decision-making process within these models remains
obscured, presenting a critical challenge for operational forecasting
communities that rely on transparency to make informed decisions. Recently
several empirical methods have been developed to explain and interpret the
decisions made by deep neural networks. These are post hoc analysis methods
(attribution methods) [12], meaning they focus on the analysis of trained
models and do not contribute to the models’ parameters while training. In this
work, we primarily focus on developing a CNN-based full-disk model for solar
flare prediction of $\geq$M1.0-class flares and evaluate and explain our
model’s performance by using three of the attribution methods: (i) Guided
Gradient-weighted Class Activation Mapping (Guided Grad-CAM) [25], (ii) Deep
Shapley Additive Explanations (Deep SHAP) [13], and (iii) Integrated Gradients
(IG) [31]. More specifically, we show that our model’s decisions are based on
the characteristics corresponding to ARs, and our models can tackle the flares
appearing on near-limb regions.

The rest of this paper is organized as follows. In Sec. 2, we present the
related work on flare forecasting. In Sec. 3, we present our methodology with
data preparation and model architecture. In Sec. 4 we provide the description
of all three post hoc explanation methods utilized in this work. In Sec. 5, we
present our experimental settings, and model evaluation, and discuss the
interpretation of our models, and in Sec. 6, we present our conclusions and
future work.

##  2 Related Work

There have been several attempts to predict solar flares using machine
learning and deep learning models. A multi-layer perceptron-based model was
applied for solar flare prediction of $\geq$C1.0- and $\geq$M1.0-class flares
in [15] by utilizing 79 manually selected physical precursors extracted from
multi-modal solar observations. A CNN-based flare forecasting model trained
with solar AR patches extracted from line-of-sight (LoS) magnetograms within
$\pm$30∘ of the central meridian to predict $\geq$C1.0-, $\geq$M1.0-, and
$\geq$X1.0-class flares was presented in [8]. Similarly, [11] also used a CNN-
based model to issue binary class predictions for both $\geq$C1.0- and
$\geq$M1.0-class flares within 24 hours using AR patches located within $\pm
45^{\circ}$ of the central meridian. It is important to note that both of
these models [8], [11] are limited to a small portion of the observable disk
in central locations (within $\pm 30^{\circ}$ to $\pm 45^{\circ}$) and thus
possess the limited operational capability.

More recently, we presented a deep learning-based binary full-disk flare
prediction model to predict $\geq$M1.0-class flares in [17] and to predict
$\geq$C4.0- and $\geq$M1.0-class flares in [18] using bi-daily observations
(i.e., two magnetograms per day) of full-disk LoS magnetograms. It is
important to note that in [18] all the instances that fall between the
$\geq$C4.0- and $\geq$M1.0-class flares were excluded in both training and
validation sets. These particular sets of instances lie on the border of two
binary outcomes and can be considered the harder-to-predict instances. These
models are still black-box and do not provide explanations on any global and
local variable importance. These explanations are important to understand the
capabilities of full-disk models in near-limb regions and improve their
trustworthiness in operational settings. In solar flare prediction, [2] used
an occlusion-based method to interpret a CNN-based solar flare prediction
model trained with AR patches. Similarly, [33] presented a deep learning-based
flare prediction model for predicting C-, M-, and X-class flares and provided
visual explanations using Grad-CAM [25], and Guided Backpropagation [28]. They
used daily observations of solar full-disk LoS magnetograms at 00:00 UT, and
their models show limitations for the near-limb flares. Moreover, in [30],
DeepLIFT [27] and IG [31] were evaluated for explaining CNN-based flare
prediction model trained using tracked AR patches within $\pm 70^{\circ}$.

This paper presents a CNN-based model to predict $\geq$M1.0-class flares,
trained with full-disk LoS magnetograms images. The novel contributions of
this paper are as follows: (i) We show an improved overall performance of a
full-disk solar flare prediction model, (ii) We utilized contemporary
attribution methods to explain and interpret the decisions of our deep
learning model, and (iii) More importantly, we show that our models can
predict solar flares appearing on difficult-to-predict near-limb regions of
the Sun.

##  3 Data and Model

Figure 1: Data distribution used in this study with four tri-monthly
partitions for training $\geq$M1.0-class flare prediction models. Note: The
length of the bars is in logarithmic scale.

We used full-disk LoS solar magnetograms obtained from the Helioseismic and
Magnetic Imager (HMI) [24] instrument onboard Solar Dynamics Observatory (SDO)
[21] available as compressed JP2 images in near real-time publicly via
Helioviewer111 Helioviewer: https://api.helioviewer.org. To enhance
computational efficiency for training the deep learning model, these
compressed images are resized to a smaller resolution of 512x512 pixels. We
sampled hourly instances of magnetogram images at [00:00, 01:00, …, 23:00]
each day from Dec 2010 to Dec 2018. We labeled our data with a prediction
window of 24 hours. The images are labeled based on the maximum peak X-ray
flux (converted to NOAA flare classes) within the next 24 hours. We collect a
total of 63,649 images and label them such that if the maximum X-ray intensity
of flare is weaker than M1.0, the observations are labeled as ”No Flare” (NF:
$<$M1.0) and $\geq$M1.0 ones are labeled as ”Flare” (FL: $\geq$M1.0). This
results in 54,649 instances for the NF-class and 9,000 instances (8,120
instances of M-class and 880 instances of X-class flares) for the FL-class.

Figure 2: The architecture of our full-disk flare prediction model.

We finally split our data into four temporally non-overlapping tri-monthly
partitions for the cross-validation experiments. This partitioning of the
dataset is created by dividing the data timeline from Dec 2010 to Dec 2018
into four partitions, where Partition-1 contains data from Jan to Mar,
Partition-2 contains data from Apr to Jun, Partition-3 contains data from Jul
to Sep, and finally, Partition-4 contains data from Oct to Dec as shown in
Fig. 1. As a result of the infrequent occurrence of $\geq$M1.0-class flares,
the dataset exhibits a significant imbalance, with the ratio of FL to NF class
being approximately 1:6.

In this work, we extend the AlexNet [10] model by concatenating a
convolutional layer at the beginning of the network to make use of the pre-
trained weights for our 1-channel input magnetogram images as the pre-trained
model requires a 3-channel image as input to the network. Our added
convolutional layer uses a 3$\times$3 kernel, size-1 stride, and outputs a
3-channel feature map which is then integrated into the standard AlexNet
architecture as shown in Fig. 2. Furthermore, to efficiently utilize the pre-
trained weights regardless of the architecture of the AlexNet model, which
expects 224$\times$224, 3-channel image as input, we use the adaptive average
pooling after feature extraction before the fully-connected layer to match the
dimension on our 1-channel, 512$\times$512 magnetogram image. Overall, our
model has six convolutional layers, three max-pool layers, one average-pool
layer, and two fully-connected layers.

##  4 Interpretation Methods

Deep learning models are often deemed black-box due to their complex
representations, resulting in interpretability, transparency, and consistency
challenges concerning the patterns they learn [12]. To address this, various
methods [34] have been proposed to interpret CNNs. One common approach is
using attribution methods, which visualize how specific parts of the input
influence the model’s decisions. Attribution methods generate attribution
vectors (heat maps) representing the contribution of each input element to the
model’s decision. These methods can be perturbation-based (e.g., Local
Interpretable Model-Agnostic Explanations (LIME) [23]), involving altering the
input and measuring the difference in output, or gradient-based, calculating
gradients via backpropagation to estimate attribution scores. While
perturbation-based methods suffer from inconsistency issues due to creating
Out-of-Distribution data [22], gradient-based methods are more robust to input
perturbations and computationally efficient [14]. Therefore, in this work, we
employed three recent gradient-based methods to assess the interpretability of
our models. By leveraging gradient-based techniques, known for their
computational efficiency and robustness compared to perturbation-based
methods, we aimed to visualize the decisions made by our model and gain
insights into the specific characteristics in a magnetogram image that trigger
the models’ decisions. These methods allowed us to cross-validate and ensure
the consistency of the explanations provided by our models, contributing to a
more reliable and robust interpretation.

Guided Grad-CAM:  The Guided Gradient-weighted Class Activation Mapping
(Guided Grad-CAM) method [25] leverages the benefits of the Grad-CAM and
guided backpropagation [28]. Grad-CAM is a model-agnostic method that uses the
class-specific gradient information flowing into the final convolutional layer
of a CNN to produce a coarse localization map of the important regions in the
image. Guided Backpropagation is based on the premise that the neurons act as
detectors of certain image features, so it computes the gradient of the output
with respect to the input, except that when propagating through ReLU
functions, it only backpropagates the non-negative gradients and highlights
the pixels that are important in the image. Attributions from Grad-CAM are
class-discriminative and localize relevant image regions; however, do not
highlight the fine-grained pixel importance as guided backpropagation [3].
Guided Grad-CAM combines the fine-grained details of guided backpropagation
with the course localization advantages of Grad-CAM and is computed as the
element-wise product of guided backpropagation with the upsampled Grad-CAM
attributions.

Deep SHAP:  SHAP values (SHapley Additive exPlanations) [13] is a method based
on cooperative game theory[26] and used to increase the transparency and
interpretability of machine learning models. SHAP shows the contribution of
each feature to the prediction of the model, it does not evaluate the quality
of the prediction itself. The contribution of each feature is calculated using
cooperative game theory and Shapley values to assess how much each feature
adds to the difference between the actual prediction and the average
prediction. For deep-learning models, Deep SHAP [13] is considered an enhanced
version of the DeepLIFT algorithm [27], where we approximate the conditional
expectations of SHAP values using a selection of baseline samples from the
dataset. The baselines typically contain a set of representative samples from
the same distribution as the input data. For each input sample, it computes
DeepLIFT attribution with respect to each baseline and averages resulting
attributions. This method assumes that input features are independent of one
another, and the explanations are modeled through the additive composition of
feature effects.

Integrated Gradients:  The last method we will analyze in this study is
Integrated Gradients (IG) [31], which quantifies feature attributions by
integrating the gradients of the model’s output along a straight-line path
from a baseline reference to the input feature under consideration. This
method requires an extra input as the baseline, representing the non-
appearance of the feature in the original image which is typically an all-zero
vector. IG is favored for its completeness property, where the sum of
integrated gradients for all features precisely equals the difference between
the model’s output for the given input and the baseline input values. This
property ensures that the feature attributions accurately represent each
feature’s individual contribution to the model output, allowing us to reliably
recover the model’s output value by summing these contributions [29].

##  5 Experimental Evaluation

###  5.1 Experimental Settings

We trained a full-disk flare prediction model with stochastic gradient descent
(SGD) as an optimizer and negative log-likelihood (NLL) as the objective
function. Our model is initialized with pre-trained weights of AlexNet Model
[10], and then we make use of a dynamic learning rate (initialized at 0.0099
and reduced 5$\%$) to further train the model to 40 epochs with a batch size
of 64. We address the class-imbalance issue using data augmentation and class
weights to the loss function. We use three augmentation techniques: vertical
flipping, horizontal flipping, and +5∘ to -5∘ rotations. We augment the data
for both classes (where the entire FL-class data are augmented three times
with three augmentation techniques and NF-class is augmented once randomly).
We then adjust class weights inversely proportional to the class frequencies
after augmentations. The use of class weights penalizes the misclassification
made in the minority class. Our models are trained as 4-fold cross-validation
experiments with each fold representing a different partition serving as the
test set. Specifically, Fold-1 corresponds to Partition-1, Fold-2 corresponds
to Partition-2, and so on.

We evaluate the performance of our models using two widely-used forecast
skills scores: True Skill Statistics (TSS, in Eq. 1) and Heidke Skill Score
(HSS, in Eq. 2), derived from the elements of confusion matrix: True Positives
(TP), True Negatives (TN), False Positives (FP), and False Negatives (FN). In
the context of our paper, the FL class is the positive outcome and NF is the
negative.

$TSS=\frac{TP}{TP+FN}-\frac{FP}{FP+TN}$ (1)

$HSS=2\times\frac{TP\times TN-FN\times FP}{((P\times(FN+TN)+(TP+FP)\times
N))},$ (2)

where $N=TN+FP$ and $P=TP+FN$.

$Recall=\frac{TP}{TP+FN}$ (3)

TSS and HSS values range from -1 to 1, where 1 indicates all correct
predictions, -1 represents all incorrect predictions, and 0 represents no
skill. In contrast to TSS, HSS is an imbalance-aware metric, and it is common
practice to use HSS in combination with TSS for the solar flare prediction
models due to the high class-imbalance ratio present in the datasets. For a
balanced dataset, these metrics are equivalent [1]. In solar flare prediction,
TSS and HSS are the preferred choices of evaluation metrics compared to
commonly used metrics in image classification (e.g., accuracy) as they ensure
a comprehensive and reliable evaluation of predictive capabilities, especially
in scenarios with imbalanced class distributions. Lastly, we report the
subclass and overall recall (shown in Eq. 3) for flaring instances (M- and
X-class) to assess the prediction sensitivity of our models in central and
near-limb regions. To reproduce this work, the source code and experimental
results can be accessed from our open-source repository [4].

###  5.2 Model Evaluation

Our models have on average TSS$\sim$0.51$\pm$0.05 and HSS$\sim$0.38$\pm$0.08,
which improves over the performance of [17] by $\sim$4% in terms of TSS
(reported 0.47$\pm$0.06) and by $\sim$3% in terms of HSS (reported
0.35$\pm$0.05) 222While there are several other works (mentioned in Sec. 2) in
solar flare prediction, the results of these models are not directly
comparable since they employ different datasets, data timelines, and data
partitioning strategies.. The detailed experimental results for each fold are
shown in Table. 1.

Table 1: A comprehensive overview of 4-fold cross-validation experiments, showing all the four outcomes of confusion matrices (TP, FP, TN, FN) evaluated on the test sets, and performance of our models in terms of two skill scores (TSS and HSS). Folds | TP | FP | TN | FN | TSS | HSS
---|---|---|---|---|---|---
Fold-1 | 1,720 | 1,943 | 10,511 | 614 | 0.58 | 0.47
Fold-2 | 1,155 | 3,083 | 10,772 | 457 | 0.49 | 0.29
Fold-3 | 1,585 | 2,668 | 11,640 | 779 | 0.48 | 0.36
Fold-4 | 1,706 | 2,241 | 11,791 | 984 | 0.47 | 0.40
Aggregated | 6,166 | 9,935 | 44,714 | 2,834 | 0.51$\pm$0.05 | 0.38$\pm$0.08
Table 2: Counts of correctly (TP) and incorrectly (FN) classified X- and M-class flares in central ($|longitude|$$\leq\pm 70^{\circ}$) and near-limb locations. The recall across different location groups is also presented. Counts are aggregated across folds. | Within $\pm$70∘ | Beyond $\pm$70∘
---|---|---
Flare-Class | TP | FN | Recall | TP | FN | Recall
X-Class | 637 | 31 | 0.95 | 157 | 55 | 0.74
M-Class | 4,229 | 1,601 | 0.73 | 1,143 | 1,147 | 0.50
Total (X&M) | 4,866 | 1,632 | 0.75 | 1,300 | 1,202 | 0.52
Figure 3: A scatterplot to quantify the performance of our models in terms of
True Positives (TP) and False Negatives (FN) for X-class flares grouped by
flare locations. The flare events beyond $\pm$70∘ longitude are represented as
near-limb events. Note: (i) Red marker is for locations with zero TP. (ii) For
some locations, TP$+$FN$<$24, given that we used hourly instances, is due to
the unavailable instances from the source.

In addition, we evaluate our results for correctly predicted and missed flare
counts for class-specific flares (X-class and M-class) in central locations
(within $\pm$70∘) and near-limb locations (beyond $\pm$70∘) of the Sun as
shown in Table 2. We observe that our models made correct predictions for
$\sim$95% of the X-class flares and $\sim$73% of the M-class flares in central
locations. Similarly, our models show a compelling performance for flares
appearing on near-limb locations of the Sun, where $\sim$74% of the X-class
and $\sim$50% of the M-class flares are predicted correctly. This is important
because, to our knowledge, the prediction of near-limb flares is often
overlooked. More false positives in M-class are expected because of the
model’s inability to distinguish bordering class flares (C4.0 to C9.9) from
$\geq$M1.0-class flares, which we have observed empirically in our prior work
[18] as well. Overall, we observed that $\sim$90% and $\sim$66% of the X-class
and M-class flares, respectively, are predicted correctly by our models.

Furthermore, given that we sample our data with a 1-hour cadence resulting in
24 instances per day unless there are gaps due to unavailable data instances,
any given flare instance is expected to be in the prediction window of 24
instances. X-class flares are relatively large flares that often dominate the
prediction window. Therefore, we analyzed the predictions on X-class flares
and observed that from a total of 45 X-class flare locations, our models
correctly predict the occurrence of a flare at least once for 44 of them, as
shown in Fig. 3. In particular, we show that the full-disk model presented in
this paper can predict flares appearing on near-limb locations of the Sun at
great accuracy, which provides a crucial addition to operational flare
forecasting systems.

###  5.3 Model Interpretation

Figure 4: A visual explanation for a correctly predicted near-limb FL-class
instance. (a) Actual magnetogram from the dataset used as the input image. (b)
Annotated full-disk magnetogram at flare start time, showing flare location
(green flag) and NOAA ARs (red flags). Overlays (GGCAM, SHAP, IG) depict the
input image overlayed with attributions, and Maps (GGCAM, SHAP, IG) showcase
the attribution maps obtained from Guided Grad-CAM, Deep SHAP, and Integrated
Gradients respectively. Figure 5: A visual explanation for an incorrectly
predicted NF-class instance. (a) Actual magnetogram from the dataset used as
the input image. (b) Annotated full-disk magnetogram at flare start time,
showing flare location (green flag) and NOAA ARs (red flags). Overlays (GGCAM,
SHAP, IG) depict the input image overlayed with attributions, and Maps (GGCAM,
SHAP, IG) showcase the attribution maps obtained from Guided Grad-CAM, Deep
SHAP, and Integrated Gradients respectively.

In this section, we present a case study, interpreting the visual explanations
generated by our model, and also discuss the implications of these
explanations in the operational forecasting scenario. For this, we use the
visualizations generated using all three post hoc explanation methods
mentioned earlier in Sec. 4 for two instances: (i) a correctly predicted (TP)
near-limb flare instance and (ii) an incorrectly predicted (FP) instance.

Firstly, we interpret the predictions of our model for a correctly predicted
X1.4-class flare observed on 2011-09-22 at 10:29:00 UTC on the East limb (note
that East and West are reversed in solar coordinates). We generate a visual
explanation using all three attribution methods. We utilized an input image
from 2011-09-22 05:00:00 UTC (approximately 5.5 hours prior to the flare
event) where the sunspot corresponding to the flare becomes visible in the
magnetogram image. Interestingly, we observed that the pixels covering the AR
on the East limb, which is responsible for the eventual X1.4 flare, are
activated, as shown in Fig. 4. Note that the location of the flare is
indicated by a green flag and all visible NOAA ARs are indicated by red flags
in Fig. 4 (b). The model focuses on specific ARs, including the relatively
smaller AR on the East limb, even though other ARs are present in the
magnetogram image. The visualization of attribution maps suggests that, for
this particular prediction, the region responsible for the flare event is
attributed as important, contributing to the consequent decision. This finding
is consistent across all three methods, corroborating the explanation’s
reliability. However, Guided Grad-CAM and Deep SHAP provide finer details by
suppressing noise compared to IG.

Similarly, to analyze a false positive case, we present an example of a C7.1
flare observed on 2014-01-06 at 00:08:00 UTC. To explain the result, we used
an input magnetogram instance from 2014-01-05 06:00:00 UTC ($\sim$18 hours
prior to the event). The model’s prediction probability for this instance
being an FL-class is $\sim$0.97. Therefore, we seek a visual explanation of
this prediction using all three interpretation methods. Upon analysis, we
observed that the prediction mainly relies on only one AR, which indeed
corresponds to the location of the eventual C7.1 flare (indicated by the green
flag) when visualized with all three attribution methods, as shown in Fig. 5.
This incorrect prediction can be attributed to the interference of the
bordering class flares mentioned in [18]. Such interference poses a problem
for binary flare prediction models. We noticed that out of 25,150 C-class
flares, 9,240 flares led to incorrect predictions, accounting for
approximately 37% of the total C-class flares in our dataset.

These two examples, although not exhaustive, carry significant implications
for operational forecasting systems. By incorporating visual explanations into
the forecasting process, in addition to providing a full-disk flare prediction
probability, we have the capability to identify potential flare event
locations among all visible ARs precisely. This is invaluable for improving
the accuracy and reliability of solar flare forecasts, aiding in effective
risk assessment and mitigation strategies. Furthermore, it provides a deeper
understanding of the underlying factors contributing to flare occurrences,
empowering researchers and space weather experts to make more informed
decisions and take timely actions to safeguard critical infrastructure and
space assets.

##  6 Conclusion and Future Work

In this work, we used three recent gradient-based methods to interpret the
predictions of our AlexNet-based binary flare prediction model trained for the
prediction of $\geq$M1.0-class flares. We addressed the highly overlooked
problem of flares appearing in near-limb locations of the Sun, and our model
shows a compelling performance for such events. Furthermore, we evaluated our
model’s predictions with visual explanations, showing that the decisions are
primarily capturing characteristics corresponding to the active regions in the
magnetogram instance. Although our model shows improved capability, still
suffers from high false positives attributed to high C-class flares. As an
extension, we plan to study the individual class characteristics to obtain a
better way of segregating these flare classes considering the background flux
and generate a new set of labels that can better address the issue with border
class flares. Furthermore, at this point, the models are only looking at the
spatial patterns in our data, and we intend to widen this work toward
spatiotemporal models to improve the performance.

Acknowledgements:  This project is supported in part under two NSF awards
#2104004 and #1931555, jointly by the Office of Advanced Cyberinfrastructure
within the Directorate for Computer and Information Science and Engineering,
the Division of Astronomical Sciences within the Directorate for Mathematical
and Physical Sciences, and the Solar Terrestrial Physics Program and the
Division of Integrative and Collaborative Education and Research within the
Directorate for Geosciences. This work is also partially supported by the
National Aeronautics and Space Administration (NASA) grant award
#80NSSC22K0272. Data used in this study is a courtesy of NASA/SDO and the AIA,
EVE, and HMI science teams and NOAA National Geophysical Data Center (NGDC).

## References

  * [1] Ahmadzadeh, A., Aydin, B., Georgoulis, M., Kempton, D., Mahajan, S., Angryk, R.: How to train your flare prediction model: Revisiting robust sampling of rare events. The ApJ Supplement Series 254(2), 23 (May 2021)
  * [2] Bhattacharjee, S., Alshehhi, R., Dhuri, D.B., Hanasoge, S.M.: Supervised convolutional neural networks for classification of flaring and nonflaring active regions using line-of-sight magnetograms. The ApJ 898(2), 98 (Jul 2020)
  * [3] Chattopadhay, A., Sarkar, A., Howlader, P., Balasubramanian, V.N.: Grad-cam++: Generalized gradient-based visual explanations for deep convolutional networks. In: 2018 IEEE Winter Conference on Applications of Computer Vision (WACV). IEEE (Mar 2018), https://doi.org/10.1109/wacv.2018.00097
  * [4] DMLab: Source Code. https://bitbucket.org/gsudmlab/explainingfulldisk/src/main/
  * [5] Falconer, D.A., Tiwari, S.K., Moore, R.L., Khazanov, I.: A new method to quantify and reduce the net projection error in whole-solar-active-region parameters measured from vector magnetograms. The ApJ 833(2), L31 (Dec 2016)
  * [6] Fletcher, L., Dennis, B.R., Hudson, H.S., Krucker, S., Phillips, K., Veronig, A., Battaglia, M., Bone, L., Caspi, A., Chen, Q., Gallagher, P., Grigis, P.T., Ji, H., Liu, W., Milligan, R.O., Temmer, M.: An observational overview of solar flares. Space Science Reviews 159(1-4), 19–106 (Aug 2011)
  * [7] Hoeksema, J.T., Liu, Y., Hayashi, K., Sun, X., Schou, J., Couvidat, S., Norton, A., Bobra, M., Centeno, R., Leka, K.D., Barnes, G., Turmon, M.: The helioseismic and magnetic imager (HMI) vector magnetic field pipeline: Overview and performance. Solar Physics 289(9), 3483–3530 (Mar 2014)
  * [8] Huang, X., Wang, H., Xu, L., Liu, J., Li, R., Dai, X.: Deep learning based solar flare forecasting model. i. results for line-of-sight magnetograms. The ApJ 856(1), 7 (Mar 2018), https://doi.org/10.3847/1538-4357/aaae00
  * [9] Ji, A., Aydin, B., Georgoulis, M.K., Angryk, R.: All-clear flare prediction using interval-based time series classifiers. In: 2020 IEEE International Conference on Big Data (Big Data). pp. 4218–4225. IEEE (Dec 2020)
  * [10] Krizhevsky, A.: One weird trick for parallelizing convolutional neural networks (2014)
  * [11] Li, X., Zheng, Y., Wang, X., Wang, L.: Predicting solar flares using a novel deep convolutional neural network. The ApJ 891(1), 10 (Feb 2020)
  * [12] Linardatos, P., Papastefanopoulos, V., Kotsiantis, S.: Explainable AI: A review of machine learning interpretability methods. Entropy 23(1), 18 (Dec 2020)
  * [13] Lundberg, S.M., Lee, S.I.: A unified approach to interpreting model predictions. In: Proceedings of the 31st International Conference on Neural Information Processing Systems. p. 4768–4777. NIPS’17, Curran Associates Inc., Red Hook, NY, USA (2017)
  * [14] Nielsen, I.E., Dera, D., Rasool, G., Ramachandran, R.P., Bouaynaya, N.C.: Robust explainability: A tutorial on gradient-based attribution methods for deep neural networks. IEEE Signal Processing Magazine 39(4), 73–84 (Jul 2022)
  * [15] Nishizuka, N., Sugiura, K., Kubo, Y., Den, M., Ishii, M.: Deep flare net (DeFN) model for solar flare prediction. The ApJ 858(2), 113 (May 2018)
  * [16] Nishizuka, N., Sugiura, K., Kubo, Y., Den, M., Watari, S., Ishii, M.: Solar flare prediction model with three machine-learning algorithms using ultraviolet brightening and vector magnetograms. The ApJ 835(2), 156 (jan 2017)
  * [17] Pandey, C., Angryk, R.A., Aydin, B.: Solar flare forecasting with deep neural networks using compressed full-disk HMI magnetograms. In: 2021 IEEE International Conference on Big Data (Big Data). pp. 1725–1730. IEEE (Dec 2021), https://doi.org/10.1109/bigdata52589.2021.9671322
  * [18] Pandey, C., Angryk, R.A., Aydin, B.: Deep neural networks based solar flare prediction using compressed full-disk line-of-sight magnetograms. In: Information Management and Big Data, pp. 380–396. Springer International Publishing (2022), https://doi.org/10.1007/978-3-031-04447-2_26
  * [19] Pandey, C., Angryk, R.A., Aydin, B.: Explaining full-disk deep learning model for solar flare prediction using attribution methods (2023), https://arxiv.org/abs/2307.15878
  * [20] Pandey, C., Ji, A., Angryk, R.A., Georgoulis, M.K., Aydin, B.: Towards coupling full-disk and active region-based flare prediction for operational space weather forecasting. Frontiers in Astronomy and Space Sciences 9 (Aug 2022), https://doi.org/10.3389/fspas.2022.897301
  * [21] Pesnell, W., Thompson, B.J., Chamberlin, P.C.: The solar dynamics observatory (SDO). Solar Physics 275(1-2), 3–15 (Oct 2011)
  * [22] Qiu, L., Yang, Y., Cao, C.C., Zheng, Y., Ngai, H., Hsiao, J., Chen, L.: Generating perturbation-based explanations with robustness to out-of-distribution data. In: Proceedings of the ACM Web Conference 2022. ACM (Apr 2022), https://doi.org/10.1145/3485447.3512254
  * [23] Ribeiro, M.T., Singh, S., Guestrin, C.: ”why should i trust you?”. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM (Aug 2016), https://doi.org/10.1145/2939672.2939778
  * [24] Schou, J., Scherrer, P.H., Bush, R.I., Wachter, R., Couvidat, S., Rabello-Soares, M.C., Bogart, R.S., Hoeksema, J.T., Liu, Y., Duvall, T.L., Akin, D.J., Allard, B.A., Miles, J.W., Rairden, R., Shine, R.A., Tarbell, T.D., Title, A.M., Wolfson, C.J., Elmore, D.F., Norton, A.A., Tomczyk, S.: Design and ground calibration of the helioseismic and magnetic imager (HMI) instrument on the solar dynamics observatory (SDO). Solar Physics 275(1-2), 229–259 (Oct 2011)
  * [25] Selvaraju, R.R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., Batra, D.: Grad-CAM: Visual explanations from deep networks via gradient-based localization. In: 2017 IEEE International Conference on Computer Vision (ICCV). IEEE (Oct 2017), https://doi.org/10.1109/iccv.2017.74
  * [26] Shapley, L.: A Value for N-Person Games. RAND Corporation (1952)
  * [27] Shrikumar, A., Greenside, P., Kundaje, A.: Learning important features through propagating activation differences (2019)
  * [28] Springenberg, J.T., Dosovitskiy, A., Brox, T., Riedmiller, M.: Striving for simplicity: The all convolutional net (2014), https://arxiv.org/abs/1412.6806
  * [29] Sturmfels, P., Lundberg, S., Lee, S.I.: Visualizing the impact of feature attribution baselines. Distill 5(1) (Jan 2020), https://doi.org/10.23915/distill.00022
  * [30] Sun, Z., Bobra, M.G., Wang, X., Wang, Y., Sun, H., Gombosi, T., Chen, Y., Hero, A.: Predicting solar flares using CNN and LSTM on two solar cycles of active region data. The ApJ 931(2), 163 (Jun 2022), https://doi.org/10.3847/1538-4357/ac64a6
  * [31] Sundararajan, M., Taly, A., Yan, Q.: Axiomatic attribution for deep networks (2017), https://arxiv.org/abs/1703.01365
  * [32] Whitman, K., Egeland, R., Richardson, I.G., Allison, C., Quinn, P., Barzilla, J., Kitiashvili, I., Sadykov, V., Bain, H.M., Dierckxsens, M., Mays, M.L., Tadesse, T., Lee, K.T., Semones, E., Luhmann, J.G., Núñez, M., White, S.M., Kahler, S.W., Ling, A.G., Smart, D.F., Shea, M.A., Tenishev, V., Boubrahimi, S.F., Aydin, B., Martens, P., Angryk, R., Marsh, M.S., Dalla, S., Crosby, N., Schwadron, N.A., Kozarev, K., Gorby, M., Young, M.A., Laurenza, M., Cliver, E.W., Alberti, T., Stumpo, M., Benella, S., Papaioannou, A., Anastasiadis, A., Sandberg, I., Georgoulis, M.K., Ji, A., Kempton, D., Pandey, C., Li, G., Hu, J., Zank, G.P., Lavasa, E., Giannopoulos, G., Falconer, D., Kadadi, Y., Fernandes, I., Dayeh, M.A., Muñoz-Jaramillo, A., Chatterjee, S., Moreland, K.D., Sokolov, I.V., Roussev, I.I., Taktakishvili, A., Effenberger, F., Gombosi, T., Huang, Z., Zhao, L., Wijsen, N., Aran, A., Poedts, S., Kouloumvakos, A., Paassilta, M., Vainio, R., Belov, A., Eroshenko, E.A., Abunina, M.A., Abunin, A.A., Balch, C.C., Malandraki, O., Karavolos, M., Heber, B., Labrenz, J., Kühl, P., Kosovichev, A.G., Oria, V., Nita, G.M., Illarionov, E., O’Keefe, P.M., Jiang, Y., Fereira, S.H., Ali, A., Paouris, E., Aminalragia-Giamini, S., Jiggens, P., Jin, M., Lee, C.O., Palmerio, E., Bruno, A., Kasapis, S., Wang, X., Chen, Y., Sanahuja, B., Lario, D., Jacobs, C., Strauss, D.T., Steyn, R., van den Berg, J., Swalwell, B., Waterfall, C., Nedal, M., Miteva, R., Dechev, M., Zucca, P., Engell, A., Maze, B., Farmer, H., Kerber, T., Barnett, B., Loomis, J., Grey, N., Thompson, B.J., Linker, J.A., Caplan, R.M., Downs, C., Török, T., Lionello, R., Titov, V., Zhang, M., Hosseinzadeh, P.: Review of solar energetic particle models. Advances in Space Research (Aug 2022), https://doi.org/10.1016/j.asr.2022.08.006
  * [33] Yi, K., Moon, Y.J., Lim, D., Park, E., Lee, H.: Visual explanation of a deep learning solar flare forecast model and its relationship to physical parameters. The ApJ 910(1), 8 (Mar 2021), https://doi.org/10.3847/1538-4357/abdebe
  * [34] Zhou, J., Gandomi, A.H., Chen, F., Holzinger, A.: Evaluating the quality of machine learning explanations: A survey on methods and metrics. Electronics 10(5), 593 (Mar 2021), https://doi.org/10.3390/electronics10050593