# Greater Transferability and Accuracy of Norm-conserving Pseudopotentials
using Nonlinear Core Corrections

Wan-Lu Li1-3 Kaixuan Chen1-3 Elliot Rossomme1,2 Martin Head-Gordon1-3 Teresa
Head-Gordon1-5 1Kenneth S. Pitzer Center for Theoretical Chemistry,
2Department of Chemistry, 3Chemical Sciences Division, Lawrence Berkeley
National Laboratory, 4Department of Chemical and Biomolecular Engineering,
5Department of Bioengineering, University of California, Berkeley,
Berkeley, California 94720, USA
<EMAIL_ADDRESS>
###### Abstract

We present an investigation into the transferability of pseudopotentials (PPs)
with a nonlinear core correction (NLCC) using the Goedecker, Teter, and Hutter
(GTH) protocol across a range of pure GGA, meta-GGA and hybrid functionals,
and their impact on thermochemical and non-thermochemical properties. The GTH-
NLCC PP for the PBE density functional demonstrates remarkable transferability
to the PBE0 and $\omega$B97X-V exchange-correlation functionals, and relative
to no NLCC, improves agreement significantly for thermochemical benchmarks
compared to all-electron calculations. On the other hand, the B97M-rV meta-GGA
functional performs poorly with the PBE-derived GTH-NLCC PP, which is
mitigated by reoptimizing the NLCC parameters for this specific functional.
The findings reveal that atomization energies exhibit the greatest
improvements from use of the NLCC, which thus provides an important correction
needed for covalent interactions relevant to applications involving chemical
reactivity. Finally we test the NLCC-GTH PPs when combined with medium-size
TZV2P molecularly optimized (MOLOPT) basis sets which are typically utilized
in condensed phase simulations, and show that they lead to consistently good
results when compared to all-electron calculations for atomization energies,
ionization potentials, barrier heights, and non-covalent interactions, but
lead to somewhat larger errors for electron affinities.

##  1 Introduction

Pseudopotentials are mathematical representations of the interactions between
explicitly treated valence electrons and effective cores that replace atomic
nuclei and inner shell electrons, allowing for a substantial reduction in
computational resources required to model large chemical and materials
systems, while also accounting for relativistic effects in some cases.1 The
Goedecker, Teter, and Hutter (GTH) formalism and optimization protocol has
produced separable, norm-conserving, Gaussian-type dual-space pseudopotentials
(PPs)2, 3 that are widely used in ab initio molecular dynamics (AIMD)
simulations with a mixed Gaussian-planewave (GPW) strategy4. The GTH PP
construction and optimization is generally performed utilizing the neutral
atomic state as a reference with spherical symmetry. To make an effective PP,
it is essential to divide the space into two regions: muffin-tin spheres
centered on the atom in a molecule or solid, and the interstitial region that
encompasses the rest of the charge density5. Knowing the scattering properties
on the surface of the muffin-tin spheres then allows for exact solutions of
the Schrödinger equation in the interstitial region.6 By satisfying the norm
conservation condition in GTH PPs7, the logarithmic derivative function
associated with the energy accurately represents the scattering properties of
a muffin-tin sphere that contains the charge distribution of the reference
configuration. Screening effects give rise to an approximately invariant
muffin-tin sphere, wherein the total electronic charge distribution remains
largely unaffected by changes in the outer chemical environment.8 The
construction of the traditional GTH2 and Hartwigsen-Goedecker- Hutter (HGH)3
PPs followed this approach, resulting in PPs that exhibit excellent
transferability for non-spin-polarized systems.

However, in a self-consistent field calculation, the charge distribution
undergoes changes when a free atom is inserted into a molecule or solid, and
the potential within the muffin-tin region will generally differ from the
potential within a muffin-tin sphere of the same radius around the reference
atom. As a result, the scattering properties undergo modifications, and the PP
constructed based on the isolated atom might not accurately reproduce the
altered scattering properties due to the new chemical environment. In spin-
polarized calculations, the concept of an invariant muffin-tin sphere is also
not applicable since the spin polarization deviates across different chemical
environments, despite the similarity in the total charge within the invariant
muffin-tin sphere. In other words, the typical assumption for PPs that assume
a linear relationship between the charge distribution and potential leads to
inaccuracies because the electron-electron interactions in the core region are
highly localized such that they exhibit a nonlinear relationship between the
electron density and the associated potential.

One approach to address this is by incorporating nonlinear core corrections
(NLCC)9 into the PP optimization. For the NLCC schemes, the spin and charge
densities within the muffin-tin sphere are not solely determined by the
valence electrons, but instead encompass the combined contributions of the
valence charge and the core charge. The NLCC method allows for a more
realistic treatment of the core electron behavior by taking into account the
strong electron-electron interactions and the asymmetric charge density
distribution near the atomic nucleus. In accordance with the principles
underlying the GTH PPs, Willand et al. in 2013 put forward a NLCC method in
which the core charge density is represented by a single Gaussian function10,
which captures the essential characteristics of the core charge. The amplitude
and width of this Gaussian core charge distribution were subsequently
optimized through a rigorous fitting procedure similar to the regular GTH PP
optimization11, 12, 13 but including not just ground state but also excited
states and ionized electronic configurations. By employing this methodology,
the NLCC-PP achieves a high level of precision in describing the atomization
energies in molecular systems compared to all-electron calculations and
reliability for high-pressure phases of crystalline solids. Together with an
adequate treatment of dispersion, such as the empirical Grimme D214 and D315
schemes, the inclusion of NLCC corrections allows for a comprehensive
description of weakly bound intermolecular interactions, achieving an average
error of approximately 0.5 kcal/mol.

Recently, we presented a study on the optimization and transferability of
small and large core GTH PPs for Density Functional Theory (DFT) functionals
applied to molecular systems and condensed phase simulations for
electrocatalysis13. We also have recently conducted a systematic exploration
of the extent of PP inconsistency errors (PPIEs)16 that arise from the
transfer of a different PP relative to the chosen level of DFT when evaluating
energy differences. To address these errors, we have developed and implemented
empirical atom- and density functional approximation specific PP
corrections.16 But in neither case have we considered the role of NLCC PPs and
their influence on the accuracy possible with different DFT functionals
outside the original parameterization of the PBE functional10 and when
combined with less complete basis sets such as the molecularly optimized
(MOLOPT) basis sets used in AIMD application studies in CP2K17.

In this work, we directly use the 2013 NLCC HGH-PP developed for the PBE
functional10, and investigate its transferability and accuracy across various
DFT functionals including the hybrid GGA PBE018, the semi-empirical hybrid
$\omega$B97M-rV19, and the meta-GGA B97M-rV20, 21, as measured on performance
for thermochemical properties, barrier heights, isomerization energies, and
non-covalent interactions. In addition, we also combine the GTH-NLCC method
with the medium-sized TZVP MOLOPT basis sets that are used to enhance the
computational efficiency in CP2K22 and compare them to the more complete
def2-TZVPPD basis set results. Although the NLCC GTH-PP results are more
transferable and yield excellent results for the hybrid functionals regardless
of basis set, the B97M-rV functional exhibits a large PPIE. Hence we undertook
a reoptimization process to enhance the accuracy of the NLCC parameters for
this meta-GGA functional which brings it into line with the other DFT/basis
set combinations. As we show below, the use of the NLCC correction improves
overall accuracy over the full range of DFT functionals and basis sets
considered here, minimizing errors across all data sets but especially for
thermochemical properties that are most relevant for applications involving
chemical reactivity.

##  2 Computational Details

### GTH PPs and Optimization Algorithm of NLCC Parameters.

The norm-conserving, dual-space, and explicitly separable GTH pseudopotential
$V^{pp}$ contains local and nonlocal parts,2, 23 where the local term is
expressed as

$\text{V}^{pp}_{loc}(r)=-\dfrac{\text{Z}_{\text{ion}}}{r}\operatorname{erf}\left(\dfrac{r}{\sqrt{2}\text{r}^{\text{pp}}_{\text{loc}}}\right)+\sum_{i}^{4}\text{C}^{\text{pp}}_{i}\left(\dfrac{r}{\text{r}^{\text{pp}}_{\text{loc}}}\right)^{2i-2}\exp\left[-\dfrac{1}{2}\left(\dfrac{r}{\text{r}^{\text{pp}}_{\text{loc}}}\right)^{2}\right]$
(1)

and erf represents the error function, $\text{Z}_{ion}$ is the ionic charge of
the atom, $\text{r}^{\text{pp}}_{\text{loc}}$ controls the range of the
Gaussian ionic charge distribution, and the $\text{C}^{\text{pp}}_{i}$ are the
coefficients to be optimized.

The nonlocal part is expressed as

$\text{V}^{\text{pp}}_{nl}(\boldsymbol{r},\boldsymbol{r}^{\prime})=\sum_{lm}\sum_{ij}\left\langle\boldsymbol{r}|\text{p}^{lm}_{i}\right\rangle\text{h}^{l}_{ij}\left\langle\text{p}^{lm}_{j}|\boldsymbol{r^{\prime}}\right\rangle$
(2)

where $\text{h}^{l}_{ij}$ denotes a scattering matrix element to be fitted and
the Gaussian-type projectors are given as

$\left\langle\boldsymbol{r}|\text{p}^{lm}_{i}\right\rangle=\text{N}^{l}_{i}\text{Y}^{lm}(r)\text{r}^{l+2i-2}\exp\left[-\dfrac{1}{2}\left(\dfrac{r}{r_{l}}\right)^{2}\right]$
(3)

where l denotes the angular momentum quantum number of subshells,
$\text{N}^{l}_{i}$ denotes a normalization constant and $\text{Y}^{lm}$ is a
spherical harmonic function. Therefore, parameters rloc and Ci in the local
term, and hij and rl in the nonlocal part need to be fitted in the GTH PP
optimization.

The nonlinear core correction is given by

$\rho_{c}(r)=c_{core}\frac{Z-Z_{ion}}{(\sqrt{2\pi}r_{core})^{3}}e^{-\frac{r^{2}}{2r_{core}^{2}}}$
(4)

in which the core charge is described by a single Gaussian function, which is
efficient for numerical integration, and is then incorporated into the
generalized Kohn-Sham total energy equation:

$E_{KS}=\sum_{i}\langle{\psi}_{i}|\\{-\frac{1}{2}[\nabla^{2}+\nabla\cdot(V_{\tau}[\rho](\mathbf{r})\nabla)]+V_{H}[\rho]+V_{xc}[\rho+\rho_{c}]+V^{pp}\\}|\psi_{i}\rangle\\\
-E_{H}[\rho]+E_{xc}[\rho+\rho_{c}]-\int
d\mathbf{r}\rho(\mathbf{r})V_{xc}[\rho+\rho_{c}](\mathbf{r})$ (5)

where $E_{xc}$, $V_{xc}$ and $V_{\tau}$ are the exchange-correlation energy,
potential, and kinetic energy term, respectively, and $V_{H}$ is the Hartree
potential. In the CP2K code17, it is assumed that the kinetic energy density
originating from the core orbitals is not included. Therefore, when performing
meta-GGA calculations, the variables used are the core charge augmented
density and density gradient, along with the kinetic energy density solely
from the valence orbitals. As a result, in our NLCC optimization for the meta-
GGA functional, we develop a model for the core density and its gradient,
adding it to the valence density. However, it is important to note that the
correction for the kinetic energy density is still missing due to this
assumption.

The amplitude ($c_{core}$) and width ($r_{core}$) of the function are
optimized based on the usual procedure of modulating initial step size, weight
for configurations in different regions such as core, valence and virtual
states, and target accuracy for Kohn-Sham eigenvalues, with reference to the
all-electron results with the consideration of scalar relativistic effects.3
However, contrary to the optimization procedure of a traditional GTH PP, which
is fitted to a single atomic configuration considered as a symmetric sphere,
the optimized NLCC PPs are parameterized not only with the ground state, but
also with excited states and ionized electronic configurations. Several atomic
properties were chosen in the fitting procedure with respect to the all-
electron (AE) calculations: 1) Atomic eigenvalues of the occupied and first
few unoccupied valence orbitals, 2) Charge within the inert region ($r_{loc}$)
of the pseudo atom matching the charge in the same region for all the
orbitals; this criterion illustrates that the PP is norm conserving for all
the configurations used in the optimization, 3) Total energy difference, and
4) Spin polarization energy of all-electron calculation is reproduced. This is
incorporated into the objective penalty function expressed in Eq. 6

$S=\sum_{n,l}\text{w}^{2}_{n,l}(\text{p}^{\text{AE}}_{n,l}-\text{p}^{\text{PP}}_{n,l})$
(6)

where $p$ denotes different atomic properties described above, $w$ is the
weight of each property, and $n$ and $l$ are principal and angular momentum
quantum numbers, respectively.

Computational and Theoretical Details. The ATOM module implemented in CP2K
package17 was used to optimize the new NLCC GTH PPs at the level of B97M-V.
The NLCC terms optimized on top of Hartwigsen-Goedecker-Hutter (HGH)
pseudopotential served as an initial guess 10. In the calculations for
different applications of the new PPs, the energy cutoff was set as 1200 Ry
with a box size of 20 Å$\times$ 20 Å$\times$ 20 Å. During the Self-Consistent-
Field (SCF) calculation, a strategy of orbital transform (OT)24, 25 is
utilized to accelerate the convergence, of which the threshold is set to be
$10^{-6}$ hartree. The reference results are obtained from all-electron
calculations with Gaussian-type def2-TZVPPD basis sets26, 27, 28.

In order to examine the NLCC correction for different chemical properties, we
chose a number of benchmark data sets. We follow the convention outlined by
Mardirossian and Head-Gordon29 to organize datasets into categories: (1)
thermochemistry (TC), (2) barrier heights (BH), (3) isomerization energies
(I), and (4) non-covalent interactions (NC), ranging from most difficult to
easier in regards demands on the PP accuracy. Note however we are not using
their MGDB84 dataset here. For the TC category we used the whole series of G2
datasets30, 31, 32 that evaluates the covalent bond formation energy, as well
as the TAE140nonMR and TAE140MR33 for easy and difficult atomization energies,
respectively. For the BH category we used HTBH3834, NHTBH3835, and WCPT2736
datasets representing hydrogen transfer barrier heights, non-hydrogen transfer
barrier heights, and barrier heights of water-catalysed proton-transfer
reactions, respectively. For the I category we used ID (isomerization energies
”difficult”) and ISOMERIZATION2033 for isomerization energies; G21EA37, 38 and
G21IP37, 38 for adiabatic electron affinities and ionization potentials of
atoms and small molecules, respectively. For the NC category we used the NCED
(non-covalent ”easy” dimers) and S6639, 40 describing binding energies between
organic molecules and biomolecules, and a NCD (non-covalent ”difficult”)
CT2041 dataset for binding energies of charge-transfer complexes.

##  3 Results

The current standard in condensed phase simulations for codes such as CP2K is
to use GTH PPs that are combined with compact MOLOPT basis sets, such as Table
1 and Figure 1(a-c) shows the resulting errors made on atomization energies
with GTH-PP/MOLOPT relative to def2-TZVPPD all-electron calculations, which is
the most representative combination of PPs and basis sets used in AIMD
applications where chemical transformations are critical. The standard PBE PP
and its transferability to the two hybrid functionals for the GTH/MOLOPT
combination is overall superior

Table 1: Errors in atomization energies predicted by GTH PPs compared to def2-TZVPPD all-electron calculations with or without NLCC corrections and using MOLOPT or def2-TZVPPD basis sets. Mean Absolute Deviations (MADS, in kcal/mol) of the different DFT functionals and basis set combinations analyzed over the G2 series30, 31, 32 and the TAE140 datasets33 divided into 124 ”easy” molecules without multi-configurational character and 16 multi-reference atomization energies. PP |  | GTH | GTH | GTH-NLCC | GTH-NLCC
---|---|---|---|---|---
Basis Set | Dataset | MOLOPT | def2-TZVPPD | MOLOPT | def2-TZVPPD
PBE | G2 | 10.13 | 15.55 | 6.06 | 1.90
TAE140nonMR | 10.20 | 12.25 | 3.36 | 1.32
TAE140MR | 12.61 | 9.92 | 2.61 | 1.70
PBE0 | G2 | 6.68 | 11.05 | 1.80 | 5.68
TAE140nonMR | 7.69 | 8.81 | 2.34 | 3.53
TAE140MR | 10.65 | 6.96 | 3.11 | 2.76
$\omega$B97X-V | G2 | 4.98 | 9.66 | 4.03 | 3.99
TAE140nonMR | 4.05 | 4.16 | 4.48 | 2.60
TAE140MR | 8.25 | 4.35 | 4.26 | 4.82
Figure 1: Nonlinear core corrections and basis set effects on atomization
energies for the G2 and TAE140 datasets. We decomposed the total TAE140
dataset into 124 TAE140nonMR and 16 TAE140MR subsets. All the GTH
pseudopotentials applied here are at PBE level. (a-c): violin diagrams of the
energy error compared with all-electron def2-TZVPPD results for density
functionals of PBE, PBE0 and $\omega$B97X-V applied to G2, TAE140nonMR and
TAE140MR data set, respectively. (d-f): distribution of the all-electron
results predicted for different functionals for G2, TAE140nonMR and TAE140MR
datasets, respectively.

to the larger basis set for molecules without multi-reference character, but
mostly mean absolute deviations (MADs) are surprisingly large, i.e. 5-10% in
the mean absolute percentage error (MAPE) normalized by the corresponding
average value of each dataset, assessed by the energy scales shown Figure
1(d-f).

To assess the errors incurred for the GTH-NLCC parameters optimized for PBE10,
Table 1 and Figure 1(a-c) shows a significant reduction in the MADs and MAPE
regardless of basis set relative to the all-electron calculations. This
correction effectively addresses the issue of underestimating the covalent
interactions arising from the original GTH/MOLOPT parameterization. For the
PBE functional, the GTH-NLCC combined with more complete basis set def2-TZVPPD
yields the best results, with the MAD generally less than 2 kcal/mol and MAPE
less than 1.2%. Such consistency with the all-electron calculation has been
corroborated in the study conducted by Willand et al.10 where they use the
standard G2-1 dataset containing only 54 molecules. However the results using
GTH-NLCC/MOLOPT, the basis set relied upon in condensed phase AIMD simulation
work, is an impressive reduction in overall error as well.

Often PPs developed for one DFT functional, typically PBE, are transferred for
use with other DFT functionals. Table 1 and Figure 1(a-c) show that the
underestimation error incurred by PBE can be alleviated by utilizing higher
hierarchical density functionals, such as hybrid PBE0 and $\omega$B97X-V
functionals that incorporate a more accurate exchange-correlation term. And
yet the hybrid functionals also benefit from the NLCC PPs transferred from
PBE. Notably, the combination of GTH-NLCC and MOLOPT basis sets exhibits a
remarkable capacity to faithfully reproduce all-electron outcomes, when
compared to the larger basis sets used in electronic structure codes. This
underscores the extensive utility of MOLOPT basis in practical computations,
as it adeptly strikes a fine balance between precision and computational
efficiency.

Table 2 and Figure 2 depict the errors in atomization energies predicted by
various combinations of PPs and basis sets when evaluated across the G2,
TAE140nonMR, and TAE140MR datasets for the semi-empirical meta-GGA B97M-rV
functional. Although the GTH-NLCC PP for PBE was highly transferable to the
two hybrid functionals, one of which is also a semi-empirical

Table 2: Errors in atomization energies predicted by GTH PPs compared to def2-TZVPPD all-electron calculations with or without NLCC corrections and using MOLOPT or def2-TZVPPD basis sets for the meta-GGa B97M-rV functional. Mean Absolute Deviations (MADS, in kcal/mol) of the different DFT functionals and basis set combinations analyzed over the G2 series30, 31, 32 and the TAE140 datasets33 divided into the 124 ”easy” molecules without multi-configurational character and 16 multi-reference atomization energies. PP | GTH | GTH | GTH-NLCC | GTH-NLCC | GTH-NLCC-OPT | GTH-NLCC-OPT
---|---|---|---|---|---|---
Basis Set | MOLOPT | def2-TZVPPD | MOLOPT | def2-TZVPPD | MOLOPT | def2-TZVPPD
Dataset |  |  |  |  |  |
G2 | 25.42 | 29.96 | 68.76 | 61.39 | 6.63 | 5.95
TAE140nonMR | 19.38 | 20.36 | 55.23 | 51.93 | 5.62 | 4.16
TAE140MR | 11.87 | 9.07 | 38.71 | 42.43 | 6.25 | 5.24
Figure 2: Energy error distribution of atomization energies for the G2 and
TAE140 datasets predicted by B97M-rV functional. The reference is all-electron
calculation with def2-TZVPPD basis sets.

functional developed by the same data-driven procedure42, we find the
transferability is quite poor for B97M-rV meta-GGA functional20, 21, yielding
large MAD errors ($>$ 50 kcal/mol) regardless of the completeness of the basis
sets when compared to all-electron calculations (Table 2). Notably two clear
trends emerge, demonstrating first that the original GTH/MOLOPT approach
significantly underestimates the atomically covalent interaction, whereas the
existing NLCC correction (GTH-NLCC-2013) exhibits an adverse and now
substantial overestimation of the interaction. By comparing the MADs predicted
by the original GTH/MOLOPT approach (Table S2), it is also apparent that under
the existing NLCC correction the lighter elements exhibit significantly larger
MADs compared to their heavier counterparts. This disparity is likely
attributed to the greater sensitivity to the core region for lighter elements,
while conversely the stronger screening effect8 exhibited by heavier elements
helps maintain the invariance of the charge density in the core region,
resulting in lower MAD values. Nonetheless, we found it imperative in this
case to reparameterize the NLCC PPs specifically for the B97M-rV functional.
After reoptimization of the NLCC parameters at the level of B97M-rV for only
the lighter main-group elements (B, C, N, O, F), leaving the original GTH
parameters in place for the heavier atoms, we get more consistent results with
a MAD less than 6 kcal/mol using either the MOLOPT or more complete basis sets
(Table 2 and Figure 2).

Figure 3: Sensitivity of electron affinity (EA) and ionization potential (IP)
to basis size and NLCC effect. In Figure 3c from B97M-rV method, the
additional bars are corresponding to the results predicted by optimized NLCC
parameters for lighter elements (B, C, N, O, F) and original GTH for heavier
elements (Al, Si, P, S, Cl).

Additionally, we investigated the transferability and errors for the standard
and NLCC PPs and basis set combinations for other thermochemical properties
including electron affinities (EA) and ionization potentials (IP) derived from
the G21EA and G21IP datasets38, 43 as shown in Figure 3. As detailed in
Figures S1 and S2, ionization involves a larger energy scale than that found
for electron affinity energies, as it is harder to remove an electron from the
highest occupied molecular orbital (HOMO) than to attach an electron into the
lowest unoccupied molecular orbital (LUMO) because the former is much deeper
than the latter relative to the vacuum level. Overall we find that IPs are
fairly insensitive to both PPs and basis set completeness for the pure GGA and
hybrid DFT functionals, whereas the B97M-rV meta GGA benefits from NLCC
optimization, after which it yields comparable results to the other DFT
functionals and regardless of basis set. These figures also highlight a
significant trend in which the EA is far more sensitive to the completeness of
the basis set, and that the NLCC has a more minor corrective effect on the
errors due to overbinding.

Table 3: Errors in non-thermal chemistry properties predicted by GTH PPs compared to def2-TZVPPD all-electron calculations with or without NLCC corrections and using TZV2P MOLOPT or def2-TZVPPD basis sets. MAD results (kcal/mol) predicted by PBE, PBE0 and $\omega$B97M-rV functionals for S66, CT20, ISOMERIZATION20, HTBH38, NHTBH38, and WCPT27. Functional | PP/Basis | S66 | CT20 | ISOMER20 | HTBH38 | NHTBH38 | WCPT27
---|---|---|---|---|---|---|---
PBE | GTH/MOLOPT | 0.12 | 0.11 | 0.63 | 1.13 | 2.00 | 0.43
GTH/def2-TZVPPD | 0.06 | 0.04 | 0.99 | 0.95 | 1.00 | 0.90
GTH-NLCC/MOLOPT | 0.14 | 0.10 | 1.57 | 0.38 | 1.84 | 0.49
GTH-NLCC/def2-TZVPPD | 0.08 | 0.05 | 0.72 | 0.61 | 0.88 | 1.15
PBE0 | GTH/MOLOPT | 0.22 | 0.12 | 0.74 | 0.83 | 1.67 | 0.44
GTH/def2-TZVPPD | 0.07 | 0.07 | 0.78 | 0.71 | 1.03 | 1.04
GTH-NLCC/MOLOPT | 0.23 | 0.12 | 0.68 | 0.37 | 1.87 | 0.45
GTH-NLCC/def2-TZVPPD | 0.12 | 0.08 | 0.78 | 0.76 | 1.47 | 1.31
$\omega$B97M-rV | GTH/MOLOPT | 0.25 | 0.11 | 0.89 | 0.52 | 1.41 | 0.62
GTH/def2-TZVPPD | 0.10 | 0.09 | 0.38 | 0.51 | 0.58 | 1.03
GTH-NLCC/MOLOPT | 0.24 | 0.11 | 0.93 | 0.23 | 1.63 | 0.57
GTH-NLCC/def2-TZVPPD | 0.11 | 0.08 | 0.37 | 0.51 | 1.03 | 1.28
B97M-rV | GTH/MOLOPT | 0.10 | 0.14 | 1.56 | 1.36 | 2.19 | 0.50
GTH/def2-TZVPPD | 0.22 | 0.09 | 1.46 | 1.20 | 1.53 | 0.80
GTH-NLCC-2013/MOLOPT | 0.30 | 0.21 | 5.74 | 3.96 | 2.40 | 2.05
GTH-NLCC-2013/def2-TZVPPD | 0.24 | 0.17 | 4.67 | 4.00 | 2.06 | 2.70
B97M-rV | GTH-NLCC-OPT/MOLOPT | 0.27 | 0.16 | 0.61 | 0.79 | 2.58 | 0.62
GTH-NLCC-OPT/def2-TZVPPD | 0.23 | 0.11 | 0.42 | 0.67 | 0.75 | 0.57

Table 3 provides an assessment of PP/basis set combinations for other non-
reactive properties including binding energies of non-covalent interactions in
organic molecules and biomolecules (S66)39, 40, binding energies of charge-
transfer complexes (CT20)41, isomerization energies (ISOMERIZATION20)33,
hydrogen transfer barrier heights (HTBH38)34, non-hydrogen transfer barrier
heights (NHTBH38)35, and barrier heights for water-catalyzed proton-transfer
reactions (WCPT27)36. In Table 3 we consider the accuracy and transferability
of the PBE PP to the other functionals. Clearly, the NLCC effect is negligible
for these non-thermochemical properties as the chemical environment is not
dramatically changed before and after the reaction and the spin polarized
state is less involved compared to the atomization process. Similar to the
other properties examined previously, the GTH-NLCC pseudopotential at the PBE
level exhibits poor transferability to the B97M-rV functional, resulting in
significantly worse results compared to the GTH/MOLOPT approach without
correction. However, with utilization of optimized NLCC parameters for the
B97M-rV functional, substantial improvements can be achieved for non-thermal
properties as well.

##  4 Discussion and Conclusion

Pseudopotentials are widely employed in theoretical chemistry to provide a
smooth potential profile within a specified core region, that permits removal
of a specified set of core electrons from explicit consideration, thereby
reducing the computational overhead and often enabling efficient inclusion of
relativistic effects. However, their utilization can introduce non-negligible
deviations, particularly in spin-polarized systems, as spin polarization
exhibits distinct variations in response to different chemical environments,
despite the apparent similarity in charge density within the invariant muffin-
tin sphere. To overcome this challenge, the NLCC strategy developed for HGH
pseudopotentials reported by Willand et al. 10 introduced a single Gaussian
function as the core charge density, thereby attaining chemical accuracy when
compared to all-electron calculations for the standard G2-1 dataset for the
PBE functional as well as for semiempirical models corrected by Grimme
dispersion terms.

In this work we tested the transferability and accuracy of the original NLCC
parameters to the GTH pseudopotential developed for the PBE GGA functional44
and associated with the TZV2P MOLOPT basis set22, to hybrid DFT functionals
including PBE018 and the range-seperated hybrid, meta-GGA $\omega$B9X-V19, as
well as the B97M-rV meta-GGA functional. Our investigation of transferability
and accuracy encompassed a wide range of both thermochemistry and non-
thermochemistry datasets, and our findings indicate that the NLCC correction
has the largest impact on atomization energies, effectively rectifying the
original GTH/MOLOPT method and resulting in an error reduction of less than 1%
when compared to all-electron calculations.

For all properties investigated, the existing NLCC parameters optimized at the
level of PBE have good transferability to PBE0 and $\omega$B97M-rV while
maintaining satisfactory accuracy using the TZV2P MOLOPT basis sets, although
electron affinities were found to be more sensitive to the completeness of
basis sets than the NLCC correction. In contrast, the GTH-NLCC PP for PBE
shows poor transferability to the meta-GGA B97M-rV functional, with
significant deviations of over 50 kcal/mol in atomization energies when
compared to all-electron calculations, thus requiring a reoptimization of the
NLCC parameters specifically for this DFT functional. We note that similar
PPIEs were found for the SCAN functional as described by Rossomme and co-
workers16, indicating that meta-GGA functionals may be unusually sensitive to
PP replacement of core electrons. Fortunately, our optimization of the GTH-
NLCC PP for B97M-rV overcame the drastic errors using the standard GTH/MOLOPT,
in which optimization was only necessary for lighter second-row elements while
heavier elements did not require any changes, which may be insightful for
improving other meta-GGA functionals. We note that the errors in the PPs for
the meta-GGA may arise from the assumption made regarding the kinetic energy
density, which neglects the contribution of the core charge and focuses solely
on the valence density (which is also an assumption in the CP2K program). This
highlights the need, and opportunity, for developing an NLCC correction for
the kinetic energy density since meta-GGAs may offer better DFT accuracy at an
affordable cost.

In addition to transferability, the GTH-NLCC PPs combined with the MOLOPT
basis sets performed overall as well as the complete basis set calculations
for PBE and hybrid functionals, and comparably for B97M-rV after optimization
of its non-transferable GTH-NLCC PP. Thus our conclusion is that the GTH-NLCC
PP/MOLOPT combination can be reliably employed in large-scale systems,
effectively reducing computational costs while maintaining good accuracy. In
summary, we believe these findings contribute to a better understanding of the
sources of error in calculations using PPs and provide a protocol for
enhancing the reliability of other DFT functionals when combined with chosen
PPs and basis sets.

Supporting Information: Electronic configurations of core and valence regions
for pseudopotentials of elements studied; Comparison of mean absolute
deviations for a specific atom type within the G2 dataset; Distributions of
MADs and all-electron calculations for IP and EA datasets; NLCC performance on
non-thermochemical properties.

This work was supported by the U.S. Department of Energy, Office of Science,
Office of Advanced Scientific Computing, and Office of Basic Energy Sciences,
via the Scientific Discovery through Advanced Computing (SciDAC) program. This
work used computational resources provided by the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility operated under Contract DE-AC02-05CH11231. W.-L.L. and
T.H.-G. acknowledge discussions with Prof. Jürg Hutter (University of Zurich)
when optimizing the NLCC parameters.

## References

  * Schwerdtfeger 2011 Schwerdtfeger, P. The Pseudopotential Approximation in Electronic Structure Theory. _ChemPhysChem_ 2011, _12_ , 3143–3155
  * Goedecker et al. 1996 Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. _Phys. Rev. B_ 1996, _54_ , 1703–1710
  * Hartwigsen et al. 1998 Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. _Phys. Rev. B_ 1998, _58_ , 3641–3662
  * LIPPERT et al. 1997 LIPPERT, B. G.; HUTTER, J.; PARRINELLO, M. A hybrid Gaussian and plane wave density functional scheme. _Mol. Phys._ 1997, _92_ , 477–488
  * Martin 2004 Martin, R. Electronic Structure–Basic Theory and Practical Methods, Cambridge Univ. _Pr., West Nyack, NY_ 2004,
  * Austin et al. 1962 Austin, B. J.; Heine, V.; Sham, L. J. General Theory of Pseudopotentials. _Phys. Rev._ 1962, _127_ , 276–282
  * Bachelet et al. 1982 Bachelet, G. B.; Hamann, D. R.; Schlüter, M. Pseudopotentials that work: From H to Pu. _Phys. Rev. B_ 1982, _26_ , 4199–4228
  * Goedecker and Maschke 1992 Goedecker, S.; Maschke, K. Transferability of pseudopotentials. _Phys. Rev. A_ 1992, _45_ , 88–93
  * Louie et al. 1982 Louie, S. G.; Froyen, S.; Cohen, M. L. Nonlinear ionic pseudopotentials in spin-density-functional calculations. _Phys. Rev. B_ 1982, _26_ , 1738–1742
  * Willand et al. 2013 Willand, A.; Kvashnin, Y. O.; Genovese, L.; Vazquez-Mayagoitia, A.; Deb, A. K.; Sadeghi, A.; Deutsch, T.; Goedecker, S. Norm-conserving pseudopotentials with chemical accuracy compared to all-electron calculations. _J. Chem. Phys._ 2013, _138_ , 104109
  * Lu et al. 2019 Lu, J.-B.; Cantu, D. C.; Nguyen, M.-T.; Li, J.; Glezakou, V.-A.; Rousseau, R. Norm-Conserving Pseudopotentials and Basis Sets To Explore Lanthanide Chemistry in Complex Environments. _J. Chem. Theory Comput._ 2019, _15_ , 5987–5997
  * Lu et al. 2021 Lu, J.-B.; Cantu, D. C.; Xu, C.-Q.; Nguyen, M.-T.; Hu, H.-S.; Glezakou, V.-A.; Rousseau, R.; Li, J. Norm-Conserving Pseudopotentials and Basis Sets to Explore Actinide Chemistry in Complex Environments. _J. Chem. Theory Comput._ 2021, _17_ , 3360–3371
  * Li et al. 2021 Li, W.-L.; Chen, K.; Rossomme, E.; Head-Gordon, M.; Head-Gordon, T. Optimized Pseudopotentials and Basis Sets for Semiempirical Density Functional Theory for Electrocatalysis Applications. _J. Phys. Chem. Lett._ 2021, _12_ , 10304–10309
  * Grimme 2006 Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _J. Comput. Chem._ 2006, _27_ , 1787–1799
  * Grimme et al. 2010 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _J. Chem. Phys._ 2010, _132_ , 154104
  * Rossomme et al. 2023 Rossomme, E.; Cunha, L. A.; Li, W.; Chen, K.; McIsaac, A. R.; Head-Gordon, T.; Head-Gordon, M. The Good, the Bad, and the Ugly: Pseudopotential Inconsistency Errors in Molecular Applications of Density Functional Theory. _J. Chem. Theory Comput._ 2023, _19_ , 2827–2841
  * Hutter et al. 2014 Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: atomistic simulations of condensed matter systems. _WIREs Comput. Mol. Sci._ 2014, _4_ , 15–25
  * Adamo and Barone 1999 Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. _J. Chem. Phys._ 1999, _110_ , 6158–6170
  * Mardirossian and Head-Gordon 2016 Mardirossian, N.; Head-Gordon, M. omegaB97M-V: A combinatorially optimized, range-separated hybrid, meta-GGA density functional with VV10 nonlocal correlation. _J. Chem. Phys._ 2016, _144_ , 214110
  * Mardirossian et al. 2017 Mardirossian, N.; Ruiz Pestana, L.; Womack, J. C.; Skylaris, C.-K.; Head-Gordon, T.; Head-Gordon, M. Use of the rVV10 Nonlocal Correlation Functional in the B97M-V Density Functional: Defining B97M-rV and Related Functionals. _J. Phys. Chem. Lett._ 2017, _8_ , 35–40
  * Mardirossian and Head-Gordon 2015 Mardirossian, N.; Head-Gordon, M. Mapping the genome of meta-generalized gradient approximation density functionals: The search for B97M-V. _J. Chem. Phys._ 2015, _142_ , 074111
  * VandeVondele and Hutter 2007 VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. _J. Chem. Phys._ 2007, _127_ , 114105
  * Krack 2005 Krack, M. Pseudopotentials for H to Kr optimized for gradient-corrected exchange-correlation functionals. _Theor. Chem. Acc._ 2005, _114_ , 145–152
  * VandeVondele and Hutter 2003 VandeVondele, J.; Hutter, J. An efficient orbital transformation method for electronic structure calculations. _J. Chem. Phys._ 2003, _118_ , 4365–4369
  * Weber et al. 2008 Weber, V.; VandeVondele, J.; Hutter, J.; Niklasson, A. M. N. Direct energy functional minimization under orthogonality constraints. _J. Chem. Phys._ 2008, _128_ , 084113
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. _Phys. Chem. Chem. Phys._ 2005, _7_ , 3297–3305
  * Rappoport and Furche 2010 Rappoport, D.; Furche, F. Property-optimized Gaussian basis sets for molecular response calculations. _J. Chem. Phys._ 2010, _133_ , 134105
  * Weigend et al. 2003 Weigend, F.; Furche, F.; Ahlrichs, R. Gaussian basis sets of quadruple zeta valence quality for atoms H–Kr. _J. Chem. Phys._ 2003, _119_ , 12753–12762
  * Mardirossian and Head-Gordon 2017 Mardirossian, N.; Head-Gordon, M. Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals. _Mol. Phys._ 2017, _115_ , 2315–2372
  * Curtiss et al. 1997 Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. Assessment of Gaussian-2 and density functional theories for the computation of enthalpies of formation. _J. Chem. Phys._ 1997, _106_ , 1063–1079
  * Pople et al. 1989 Pople, J. A.; Head‐Gordon, M.; Fox, D. J.; Raghavachari, K.; Curtiss, L. A. Gaussian‐1 theory: A general procedure for prediction of molecular energies. _J. Chem. Phys._ 1989, _90_ , 5622–5629
  * Curtiss et al. 1990 Curtiss, L. A.; Jones, C.; Trucks, G. W.; Raghavachari, K.; Pople, J. A. Gaussian‐1 theory of molecular energies for second‐row compounds. _J. Chem. Phys._ 1990, _93_ , 2537–2545
  * Karton et al. 2011 Karton, A.; Daon, S.; Martin, J. M. W4-11: A high-confidence benchmark dataset for computational thermochemistry derived from first-principles W4 data. _Chem. Phys. Lett._ 2011, _510_ , 165–178
  * Zhao et al. 2005 Zhao, Y.; Lynch, B. J.; Truhlar, D. G. Multi-coefficient extrapolated density functional theory for thermochemistry and thermochemical kinetics. _Phys. Chem. Chem. Phys._ 2005, _7_ , 43–52
  * Zhao et al. 2005 Zhao, Y.; González-García, N.; Truhlar, D. G. Benchmark Database of Barrier Heights for Heavy Atom Transfer, Nucleophilic Substitution, Association, and Unimolecular Reactions and Its Use to Test Theoretical Methods. _J. Phys. Chem. A_ 2005, _109_ , 2012–2018
  * Karton et al. 2012 Karton, A.; O’Reilly, R. J.; Radom, L. Assessment of Theoretical Procedures for Calculating Barrier Heights for a Diverse Set of Water-Catalyzed Proton-Transfer Reactions. _J. Phys. Chem. A_ 2012, _116_ , 4211–4221
  * Goerigk and Grimme 2010 Goerigk, L.; Grimme, S. A General Database for Main Group Thermochemistry, Kinetics, and Noncovalent Interactions Assessment of Common and Reparameterized (meta-)GGA Density Functionals. _J. Chem. Theory Comput._ 2010, _6_ , 107–126
  * Curtiss et al. 1991 Curtiss, L. A.; Raghavachari, K.; Trucks, G. W.; Pople, J. A. Gaussian‐2 theory for molecular energies of first‐ and second‐row compounds. _J. Chem. Phys._ 1991, _94_ , 7221–7230
  * Řezáč et al. 2011 Řezáč, J.; Riley, K. E.; Hobza, P. Extensions of the S66 Data Set: More Accurate Interaction Energies and Angular-Displaced Nonequilibrium Geometries. _J. Chem. Theory Comput._ 2011, _7_ , 3466–3470
  * Lao and Herbert 2015 Lao, K. U.; Herbert, J. M. Accurate and Efficient Quantum Chemistry Calculations for Noncovalent Interactions in Many-Body Systems: The XSAPT Family of Methods. _J. Phys. Chem. A_ 2015, _119_ , 235–252
  * Steinmann et al. 2012 Steinmann, S. N.; Piemontesi, C.; Delachat, A.; Corminboeuf, C. Why are the Interaction Energies of Charge-Transfer Complexes Challenging for DFT? _J. Chem. Theory Comput._ 2012, _8_ , 1629–1640
  * Mardirossian and Head-Gordon 2014 Mardirossian, N.; Head-Gordon, M. omegaB97X-V: A 10-parameter, range-separated hybrid, generalized gradient approximation density functional with nonlocal correlation, designed by a survival-of-the-fittest strategy. _Phys. Chem. Chem. Phys._ 2014, _16_ , 9904–9924
  * Goerigk and Grimme 2011 Goerigk, L.; Grimme, S. Efficient and Accurate Double-Hybrid-Meta-GGA Density Functionals—Evaluation with the Extended GMTKN30 Database for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. _J. Chem. Theory Comput._ 2011, _7_ , 291–309
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868


# Predicting the energetic proton flux with a machine learning regression
algorithm

Mirko Stumpo INAF-Istituto di Astrofisica e Planetologia Spaziali, 00133 Roma,
Italy  Monica Laurenza INAF-Istituto di Astrofisica e Planetologia Spaziali,
00133 Roma, Italy  Simone Benella INAF-Istituto di Astrofisica e Planetologia
Spaziali, 00133 Roma, Italy  Maria Federica Marcucci INAF-Istituto di
Astrofisica e Planetologia Spaziali, 00133 Roma, Italy

###### Abstract

The need of real-time of monitoring and alerting systems for Space Weather
hazards has grown significantly in the last two decades. One of the most
important challenge for space mission operations and planning is the
prediction of solar proton events (SPEs). In this context, artificial
intelligence and machine learning techniques have opened a new frontier,
providing a new paradigm for statistical forecasting algorithms. The great
majority of these models aim to predict the occurrence of a SPE, i.e., they
are based on the classification approach. In this work we present a simple and
efficient machine learning regression algorithm which is able to forecast the
energetic proton flux up to 1 hour ahead by exploiting features derived from
the electron flux only. This approach could be helpful to improve monitoring
systems of the radiation risk in both deep space and near-Earth environments.
The model is very relevant for mission operations and planning, especially
when flare characteristics and source location are not available in real time,
as at Mars distance.

Solar Proton Events — Radiation environment — Space Weather forecasting —
Machine Learning

##  1 Introduction

Solar Proton Events (SPEs) are pronounced enhancements of the energetic proton
flux measured by instruments placed on different space probes across the
Heliosphere. Solar protons can reach high energies, say tens of GeVs, as a
consequence of different acceleration processes occurring at the Sun in
association with transient phenomena like solar flares and coronal mass
ejections (CMEs; Kahler et al., 1984; Shea & Smart, 1990; Aschwanden, 2002;
Iucci et al., 2005). Then, particles travel along interplanetary magnetic
field lines and can produce a geoeffective SPE that can be detected by
instruments placed on Earth-orbiting satellites, such as the Geostationary
Operational Environmental Satellite (GOES). Geoeffective SPEs are defined
operationally as $\geq$S1 events, i.e., when the $>$10 MeV proton flux
overcomes the threshold of 10 pfu. When SPEs reach relativistic energies,
secondary particles are produced in the Earth atmosphere, causing sharp
increases in the count rates of ground-based particle detectors, called ground
level enhancements (Forbush, 1946). SPEs have been traditionally divided in
impulsive and gradual, depending on their time profile after the onset and on
different properties such as composition and acceleration mechanisms.
Nevertheless, a more complex scenario emerged after the launch of Advanced
Composition Explorer (ACE) and Solar and Heliospheric Observatory (SOHO) in
the late 90s, when hybrid events were identified, possessing properties of the
both fundamental types, leading to an expanded classification system (e.g.,
Cohen et al., 1999; Kocharov & Torsti, 2002; Cane et al., 2003; Cliver, 2008).
The primary acceleration source producing protons with energies $>$10 MeV is
still a subject of debate (e.g., Zhang et al., 2021, and references therein),
with theories suggesting CME-driven shocks or post-CME reconnection processes
(Cane et al., 2003; Reames, 2000; Desai & Giacalone, 2016).

From a Space Weather perspective, SPEs pose significant hazards within the
interplanetary and near-Earth space (Eastwood et al., 2017). Indeed, they can
damage spacecraft instruments (Iucci et al., 2005), interfere with radio
communications in Earth’s atmosphere and spacecraft operations (Iucci et al.,
2005). Moreover, they present radiation threats to astronauts (Furukawa et
al., 2020), and even to airline crews, and passengers that traverse polar
routes when they reach relativistic energies (Jones et al., 2005). Therefore,
the development of warning systems is vital to predict the occurrence of SPEs
and ensure an effective mitigation against them. Such systems would involve
real-time monitoring of the solar activity, and/or advanced modeling of SPE
generation and transport, as well as implementing dedicated safety protocols
for the protection of technologies and human health. The need of forecasting
systems or at least monitoring of the radiation environment is increasingly
growing also in view of the new era of sustainable human space exploration
(e.g. Moon and Mars; Zhang et al., 2023; Fogtman et al., 2023). Indeed, no
natural shielding like the Earth magnetosphere is present in the deep space.
Radiation exposure may cause several diseases to humans (Furukawa et al.,
2020), including the induction of cataract (Cucinotta et al., 2001), cancer
(Sridharan et al., 2015; Cucinotta, 2015; Chang et al., 2016), damages to
nervous and heart systems (Jandial et al., 2018; Hughson et al., 2018) and
acute radiation syndrome (Wu et al., 2009; Carnell, 2021). Recent studies have
identified that space radiations also affect plants, which are an important
part of the environment that would support humans for the long-term habitation
of spacecraft (Ferl et al., 2002; Arena et al., 2014). Therefore, the planning
of this type of missions must include a radiation risk assessment, including
the warning system for SPEs (Fogtman et al., 2023).

The prediction of SPEs constitutes a challenging task due to different factors
at the origin of the particle acceleration, including solar conditions, seed
particles, shock geometry, coronal field configuration, and interplanetary
transport. Moreover, a main hurdle for SPE forecasting is the short transit
time of protons travelling from Sun to Earth, going from $\sim$15 min up to
few hours in the energy range between 10 MeV and 20 GeV. Several models
devoted to the forecasting of SPEs have been introduced in the last decades,
as reported in a comprehensive review by Whitman et al. (2023). SPE
forecasting models can be mainly divided into two classes: physics-based
models, and empirical models. In recent years, the latter class of models has
also taken advantage from the application of machine learning algorithms,
thanks to which it is possible to recognize hidden characteristics in the data
and use them to automatically adapt the prediction to changing physical
conditions.

Physics-based models usually simulate processes of particle
release/acceleration at the source and transport/diffusion in the
interplanetary space. A significant shortcoming of these models is the high
computational cost needed to achieve comprehensive and accurate numerical
simulations, thus making them operationally challenging. In addition, their
performances strongly depend on the precision of the input data and on the
validity of the underlying physical assumptions. In recent years, significant
efforts have been made with the goal of improving the performance of physics-
based models and making them suitable for operational forecasting. This is the
case of the SEPMOD for instance, which approximates the observed SPEs through
a simple shock source description and a scatter-free particle propagation
(Luhmann et al., 2007, 2010), later integrated with the Wang–Sheeley–Arge
coronal model (Arge et al., 2004) and with the 3D Enlil heliospheric model
(Odstrcil, 2003), allowing for simulating the ICME propagation (Palmerio et
al., 2022).

Empirical models primarily rely on identifying statistical relations between
SPEs and possible precursors. These models leverage historical data and
statistical analysis to make predictions, instead of performing numerical
simulations grounded on the knowledge of the physical processes involved.
Hence, models belonging to this class are computationally less demanding and
more suitable in operational settings. Early models to predict SPEs include:
the “proton prediction system” (PPS76; Smart & Shea, 1979, 1989; Kahler et
al., 2007) based on solar flare parameters (microwave/X-ray flux, flare
location, etc.); the Protons model (Balch, 2008, currently in use on the NOAA
Space Weather Prediction Center) using also the occurrence of radio type II
and type IV bursts as input parameters, being indicative of the presence of a
CME driven shock. The Empirical model for Solar Proton Events Real Time Alert
(ESPERTA; Laurenza et al., 2009; Alberti et al., 2017) performs a logistic
regression analysis on three input parameters, namely, the flare location, the
1–8 Åsoft X-ray (SXR) fluence, that is indicative of flare strength, and the 1
MHz radio fluence, which allows for the identification of Type III radio
bursts, indicative of fast electron beams escaping from the acceleration site
into the interplanetary space. ESPERTA aims to provide timely warnings within
a time interval of 10 minute following the SXR peak time of flares of $\geq$M2
class. It has also been adapted to forecast only the largest radiation storms
($\geq$S2) produced by $\geq$100 pfu SPEs (Laurenza et al., 2018), and
reinterpreted within a machine learning framework (Stumpo et al., 2021;
Laurenza et al., 2024). Other empirical models incorporate also information
about CMEs (e.g., Papaioannou et al., 2016; St. Cyr et al., 2017; Anastasiadis
et al., 2017; Papaioannou et al., 2018, 2022) or exploit the arrival times of
relativistic electrons and/or high-energy protons at 1 AU (e.g., REleASE and
UMASEP; Posner, 2007; Núñez, 2011). The REleASE model initially relied only on
SOHO data, and it has been adapted later to use also ACE observations
(Malandraki & Crosby, 2018). The UMASEP scheme, introduced by Núñez (2011,
2015), is based on the lag-correlation between the time derivatives of the SXR
flux and the near-Earth proton differential fluxes (e.g., those provided by
GOES). Machine learning models, such as decision trees, random forest and
neural networks, have been introduced in the landscape of SPE forecasting
models (Boubrahimi et al., 2017; Engell et al., 2017; Núñez & Paul-Pena, 2020;
Lavasa et al., 2021; Hosseinzadeh et al., 2024; Rotti et al., 2024).

The models introduced so far aim to predict mainly whether a SPE will occur or
not and/or some SPE properties. Hence, generally they do not provide the
forecast of the proton flux, which is a slightly different problem in terms of
model development and search for input features. Several authors privileged
the flux prediction, instead of making classification, providing a timely
estimate of the proton flux and also providing the possibility of more
detailed analysis in an operational setting. Such methodologies have recently
been also introduced in the time series reconstruction of cosmic-ray intensity
in the heliosphere (e.g., see Sabbatini & Grimani, 2022; Sabbatini et al.,
2024). Concerning the SPE forecasting, recently Nedal et al. (2023) employed
the bidirectional long short-term memory (bi-LSTM) model architecture in order
to make predictions of the daily $>10$Mev, $>30$MeV and $>60$MeV proton flux.
They take as input features the sunspot number, the F10.7 index, the X-ray
flux, the solar wind speed, and the interplanetary magnetic field magnitude.
Laurenza et al. (2024) introduced an LSTM-based algorithm for predicting the
GOES $>10$ MeV proton flux 1 hour in advance based on measurements of the
electron flux data in the energy range 0.25–0.7 MeV obtained from the Electron
Proton Helium INstrument (EPHIN) on board the Solar and Heliospheric
Observatory (SOHO). Here, authors compare the scores of the LSTM-based model
with the ESPERTA classifier, showing that the former significantly improves
the performance, lowering the false alarm rate, albeit the reduced warning
time due to the introduction of the electron flux as main input feature.

The great majority of the existing models reviewed so far are based on flare
and cme precursors/characteristics. However, warning systems for supporting
space exploration missions are even more challenging because flare precursors
for predicting SPEs cannot be always exploited. For example at Mars distance,
the nominal magnetic field line connection to the Sun is at $\sim$ W90,
implying that the full flare signature may not be visible as the source region
could be often located behind the solar limb (Posner & Strauss, 2020). In
addition, communications delays or gaps require innovative and high-
performance methods entirely based on in situ data for future space missions
designed for Space Weather activities and beyond. With this aim, we introduce
a novel machine learning model based on the random forest (RF) algorithm for
energetic proton flux forecasting. The main goal of this work is to create a
model capable of giving an estimation of the energetic proton flux in a
desired energy range based on near relativistic electron flux. The model
production pipeline can be adapted to build and train similar models involving
integral proton flux (e.g. $>$10 MeV) or differential proton flux. This makes
our approach very helpful in view of space missions carrying on board particle
detectors that allows for electron and proton flux measurements, with the
option of providing a dedicated SPEs prediction system to be directly
implemented on board. Therefore, this work contribution may be twofold.
Firstly it contributes to the need of efficient and accurate algorithms for
predicting the proton flux forward in time directly onboard. Secondly, this
work may be helpful for designing warning systems not based on flare
characteristics, thus potentially supporting space exploration missions.

##  2 Dataset

In real-time monitoring and space weather related contexts, several space
missions collected a huge amount of data during the last solar cycles which
are sufficient for training machine learning algorithms. In this work we used
data from the COSTEP/EPHIN instrument placed on board the ESA/NASA SOHO
spaceprobe, orbiting around the L1 libration point (see Müller-Mellin et al.,
1995, for a description of the COSTEP experiment). The EPHIN telescope is
meant for measuring electrons in the energy range from 250 keV to $>$8.7 MeV,
and protons and alpha particles from 4 MeV/nucleon to $>$53 MeV/nucleon (see
Table 1) with a 1 minute cadence. Note that SPEs definition is usually based
on GOES integral proton flux ($>10$ MeV), thus our results are not directly
comparable with the intensities reported by GOES. However, in order to show
the capabilities of our approach, we integrated the COSTEP/EPHIN differential
channels and merged the P8, P25 and P41 channels. We compared the resulting
time series with the GOES $>10$ MeV integral channel and they are very similar
in terms of both flux value and profile (not shown). Thus the energetic proton
flux we are going to consider in this work is in the range of energies between
7.8 MeV and 53 MeV. The COSTEP/EPHIN instrument is also provided of an
integral channel measuring together i) $>53$ MeV protons and alpha particles;
ii) $>8.7$ MeV electrons. We did not include this channel in the merged one
because of the presence of both electrons and alpha particles together with
protons.

Channel | Species | Energy range
---|---|---
E150 | Electron | (0.25 - 0.70) MeV
E300 | Electron | (0.70 - 3.00) MeV
E1300 | Electron | (2.60 - 6.20) MeV
E3000 | Electron | (4.80 - 10.40) MeV
P4 | Proton | (4.30 - 7.80) MeV
P8 | Proton | (7.80 - 25.00) MeV
P25 | Proton | (25.00 - 40.90) MeV
P41 | Proton | (40.90 - 53.00) MeV
INT | p + e + He |  $>$53 MeV
Table 1: Proton/electron differential channels and the integral channel (INT)
of the EPHIN instrument. Adapted from Müller-Mellin et al. (1995).

In order to make reliable and efficient predictions for increases in the
energetic proton flux due to solar particles we need information about the
magnetic connection between the source and the observer. Some classification
algorithms, based on flare precursors (e.g., Laurenza et al., 2009; Stumpo et
al., 2021), use the heliolongitude of the flare, inferred from observations of
the solar flare region in the $H\alpha$ spectral line. In this work the
information on magnetic connectivity is indirectly provided by the earlier
arrival of quasi-relativistic electrons. Therefore, we base our model on the
hypothesis that electron and proton events share the same source region.
Figure 1 shows an event included in the dataset used in this work. The onset
of the event in the near-relativistic electron data anticipates the one in the
proton data at several MeVs, as expected. Furthermore, the time profiles of
electrons and protons are very similar (see also Posner, 2007).

Figure 1: EPHIN proton/electron differential channels P8 (blue), E150 (orange)
and E300 (green) during the event of October 25, 2000. The black vertical
lines mark the onset in proton and electron data. In this event, the electrons
measured in E150 and E300 differential channels arrive 01:30 hours earlier
than protons in the P8 channel. Furthermore, the time profiles are very
similar.

###  2.1 Feature selection

The target of the model is the logarithm of the proton flux integrated in the
channels P8, P25 and P41 of the COSTEP/EPHIN instrument (see Table 1). From
now on, for the sake of simplicity, we will refer to logarithm of the
integrated proton flux at such energies as only proton flux. The physical
observables used in the training set are:

  * •

electron differential channel E150 (0.25-0.70) MeV

  * •

electron differential channel E300 (0.67-3.00) MeV

  * •

electron differential channel E1300 (2.64-6.18) MeV

  * •

electron differential channel E3000 (4.80-10.4) MeV

  * •

derivative of each electron differential channel

  * •

derivative of the proton flux

All these quantities have been log-transformed in order to reduce the range of
variability of data between the time profiles produced by background protons
and SPEs. The derivatives of i) each electron differential channel and ii) the
proton flux, are here considered in the training dataset in order to provide
the model of the information about the rise parameter/grow rate (and thus
about the magnetic connectivity, indirectly Posner, 2007). Note that the log-
transformation has been applied to the proton and electron fluxes before
calculating the derivatives in order to avoid negative values. Each feature in
principle may appear with different temporal lags with respect to the proton
flux (as also explained in Section 3). Finally, the training dataset is
completed by the average and the standard deviation of each physical
observable, computed in moving windows of different widths. In particular, we
set the number of windows to be the same as the number of lags
$n_{\text{lags}}$ and their temporal widths are the same as the lags.
Therefore, the total number of features is given by $N=27\cdot
n_{\text{lags}}$, where $27=9+9\cdot 2$, i.e., obtained by summing up the 9
physical observables and the 9 physical observables multiplied by the 2
statistic measures (mean and standard deviation) considered for each feature.
However, even though the procedure described above is quite general, we choose
to use only one lag (1 hour) in order to avoid overfitting. Thus we have a
total of $N=27$ features and the lookback is $t_{b}=1$ hour. Regularization
techniques are necessary for longer lookbacks, which are not currently
implemented in our pipeline.

In principle any lookforward can be chosen. Our choice ($t_{f}=1$ hour) is
motivated by the onset delay analysis reported by Posner (2007), according to
which the forecasting time of the energetic proton flux based only on the
relativistic electron flux is up to 1 hour. In general, there is a trade-off
between the lookback and the lookforward of the model. The prediction becomes
smoother and more accurate as the number of lags of the features increases,
but the model tends to overfit the training dataset. For example a time offset
between the target and the prediction may appear.

$\mathbf{n_{tree}}$ | 50
---|---
max. depth | 50
criterion | Squared error
max features | sqrt
min. samples split | 10
min. samples leaf | 10
Table 2: Hyperparameters of the RF regressor used in this work, referred to
the Python implementation given by scikit-learn (version-1.2.2; Pedregosa et
al., 2011). All the non-mentioned hyperparameters are kept to the default
values.

##  3 Random Forest Regressor

The algorithm used to build the model is the Random Forest (RF). The RF is an
ensemble learning algorithm that operates by creating an ensemble of decision
trees through bootstrapped sampling and random feature selection. A decision
tree is a hierarchical structure used in both classification and regression
tasks. The training process of a decision tree involves recursively splitting
the data according with feature values to make predictions. These splits are
determined by optimization criteria such as the Gini impurity for
classification or the mean squared error (MSE) for regression. Each internal
node of a decision tree represents a feature and a decision rule associated
with the feature, while each leaf node represents a prediction. For example,
given a feature vector $\mathbf{X}$, the root node may be asking the value of
the feature $i$-th feature, $x_{i}$. The algorithm reads the value $x_{i}$ and
compare it to the threshold value $\epsilon$. If $x_{i}>\epsilon$, then the
algorithm decides to go through the branch 1 composed by additional nodes
checking the other features otherwise it goes through the branch 2 composed by
some other nodes. The hierarchical structure of each decision tree (i.e. the
value of the thresholds $\epsilon_{i}$ and the nodes) is automatically created
during the learning phase. However, a decision tree alone is susceptible to
overfitting when the tree becomes too complex and closely fits the training
data. The RF algorithm has been introduced with the aim of overcoming this
problem. The idea of the RF algorithm is to create an ensemble of independent
decision trees, each one making its own prediction (after the learning phase).
Then the final prediction is the arithmetic average of decision trees outcomes
(see Figure 2 for a sketch).

The RF algorithm can be effectively adapted to timeseries forecasting, by
redefining the problem as a regression task. To achieve this, we create a
modified dataset where the target variable $y$ is paired with lagged versions
of each single feature $x$, shifted forward in time by a fixed amount
$t^{\prime}=n$. In this setup, $y_{t}$ is compared to $x_{t-n}$, $y_{t+1}$ is
compared to $x_{t-n+1}$, and so forth. In essence, we establish a relationship
between the current target variable $y_{t}$ and the historical values of each
feature $x$, allowing the model to make predictions of $y$’s values $n$ time
steps into the future. By training the RF on this modified dataset, the model
learns the dependencies between the past values of each feature $x$ and future
values of $y$ and automatically builds the hierarchical structure of each
decision tree. Formally the model is expressed by a function $f$ mapping the
total feature vector $\mathbf{X}_{t,t_{b}}$ into the target $y_{t}$, i.e.,

$y_{t+t_{f}}=f(\mathbf{X}_{t,t_{b}};\mathbf{W}),$ (1)

where $t_{f}$ is the lookforward and $t_{b}$ is the lookback. $\mathbf{W}$ is
a set of parameters that define the structure of the trees (nodes, thresholds
and the associated decision rules) and are learnt from the training dataset.
Note that, for a fixed time $t$ the features are represented by the vector
$\mathbf{X}$, which incorporates the lags of each physical observable (i.e.,
past values of the observables are flattened in a $N$-dimensional row vector,
where $N$ is the total number of the features). Practically, in this way, we
are able to define the lookback $t_{b}$ for the RF. Therefore, for example,
suppose to have a single feature $x$. Then the features vector would be
$\mathbf{X}_{t,t_{b}}=(x_{t},x_{t-t_{1}},...,x_{t-t_{b}})$, where $t_{i}$
denotes the i-th lag of the observable $x$, i.e.,
$t_{i}\in\left[t_{1},t_{b}\right]$. The whole vector $\mathbf{X}_{t,t_{b}}$ is
associated to the single output $y_{t+t_{f}}$, i.e.,

$\mathbf{X}_{t,t_{b}}\mapsto y_{t+t_{f}}.$ (2)

In principle, the lags $t_{i}$ need not to be equally spaced.

Finally, note also that the function defined in Eq. (1) implicitly depends on
a set of hyperparameters which control for example the number of trees in the
forest or the loss function used for optimizing the weights $\mathbf{W}$. In
the following we are going to consider the hyperparameters reported in Table 2
and a single-output model.

for tree=l sep=2em, s sep=1em, anchor=center, inner sep=0.4em, fill=blue!50,
circle, where level=2no edge [ Training Dataset, node box [Sample and feature
bagging, node box, alias=bagging, above=3em
[,red!70,alias=a1[[,alias=a2][]][,red!70,edge label=node[above=1ex,red
arrow][[][]][,red!70,edge label=node[above=1ex,red arrow][,red!70,edge
label=node[below=1ex,red arrow]][,alias=a3]]]] [,red!70,alias=b1[,red!70,edge
label=node[below=1ex,red arrow][[,alias=b2][]][,red!70,edge
label=node[above=1ex,red arrow]]][[][[][,alias=b3]]]] [ $\dots$ ,scale=2,no
edge,fill=none,yshift=-3em] [,red!70,alias=c1[[,alias=c2][]][,red!70,edge
label=node[above=1ex,red arrow][,red!70,edge label=node[above=1ex,red
arrow][,alias=c3][,red!70,edge label=node[above=1ex,red arrow]]][,alias=c4]]]]
] [tree box, fit=(a1)(a2)(a3)] (t1) ; [tree box, fit=(b1)(b2)(b3)] (t2) ;
[tree box, fit=(c1)(c2)(c3)(c4)] (tn) ; [below right=0.5em, inner sep=0pt] at
(t1.north west) Tree 1; [below right=0.5em, inner sep=0pt] at (t2.north west)
Tree 2; [below right=0.5em, inner sep=0pt] at (tn.north west) Tree $n$;
t1.southwest)--tn.south east) node[midway,below=4em, node box] (mean) Average;
[below=3em of mean, node box] (pred) Prediction; [black arrow=5mm4mm]
(bagging) – (t1.north); [black arrow] (bagging) – (t2.north); [black
arrow=5mm4mm] (bagging) – (tn.north); [black arrow=5mm5mm] (t1.south) –
(mean); [black arrow] (t2.south) – (mean); [black arrow=5mm5mm] (tn.south) –
(mean); [black arrow] (mean) – (pred);

Figure 2: RF algorithm diagram adapted from Riebesell & Bringuier (2023). The
paths highlighted in red represent the series of decisions made by each tree
in the forest to make the final prediction. The RF prediction is the average
of the trees final decisions. Each node represents a logical operation on a
single feature. For example, a node may be asking if the variable E300 $<$ 100
(cm2 s sr MeV)-1. If yes go to the next node on the right, otherwise go to the
next node on the left. The graph and the logical operations are inferred from
data during the training phase of the model.

##  4 Results

To asses the model performance, the dataset has been split in two parts. The
first part used for the learning phase (training dataset) and the second one
used for testing only (test dataset). Historical electron and proton flux
timeseries have been considered in the analysis, i.e., without separating
transient events (SPEs) from the background, thus enabling the model to
accurately learn the differences between them. The dataset has been used with
full cadence sampling of 1 minute. The training dataset ranges from 1996
January 1 to 2012 February 12, whereas the test dataset ranges from 2012
February 13 to 2017 September 2. Data from 2017 October 4 onwards are not
included because of instrument reconfiguration (see http://www2.physik.uni-
kiel.de/SOHO/phpeph/EPHIN.htm).

Figure 3 shows a sample from the test dataset. The inset panel of the same
Figure shows a closer view of the event of 2012 March 7. All the SPEs observed
in the time intervals depicted in Figure 3 are precisely reproduced from the
model in terms of both proton flux value and timing, i.e., event onsets and
flux values. By inspecting the SPE of March 2012 (inset of Figure 3), the
overall trend of the proton flux is nicely reproduced by the model, which
provides also realistic fluctuations overlaid on a global trend of the event.
Hence, the features selected for the model provide extensive information about
energetic proton flux characteristics. However, such features are mainly
derived from electron flux measurements and the proton derivative, thus being
informative of the location/magnetic connectivity of the probe with the Sun
and of the particle acceleration processes at the source. Signatures of the
proton local energization associated with shock propagation, i.e., the so
called energetic storm particles (ESP), are then not expected in the predicted
timeseries. As a confirmation of this, the EPHIN in situ measurements of the
March 2012 event reported in the inset panel of Figure 3 exhibit the typical
peak of an ESP event, which, conversely, is not present in the predicted
proton flux.

To better asses the quality of the predictions we refer to: i) the root mean
squared error (RMSE), which is in the same units as the target; ii) the
coefficient of determination $R^{2}$; iii) the cross correlation between the
target timeseries and the predicted one. The i) and ii) provide the
performance metrics which quantify the information explained by the model,
i.e., how well the details of the target time series (e.g. trend and
fluctuations) are reproduced by the predicted one. The RMSE is defined as the
mean distance between the target $y_{t_{k}+n}$ and the predicted value
$\hat{y}_{t_{k}+n}$

$\text{RMSE}=\sqrt{\frac{\sum_{k=1}^{N}(y_{t_{k}+n}-\hat{y}_{t_{k}+n})^{2}}{N}},$
(3)

and provides a measure of the mean amplitude of the fluctuations around the
target $y_{t_{k}+n}$.

The coefficient of determination $R^{2}$ is in essence the explained variance
of the predictions and is defined as

$R^{2}=1-\frac{\sum_{k=1}^{N}(y_{t_{k}+n}-\hat{y}_{t_{k}+n})^{2}}{\sum_{k=1}^{N}(y_{t_{k}+n}-\bar{y})^{2}},$
(4)

where $\bar{y}$ is the mean of the target. The value of $R^{2}$ computed on a
dataset (not used for fitting the model) ranges from $-\infty$ to 1, where in
particular $R^{2}=1$ means perfect predictions.

Finally, we computed the iii) cross correlation for measuring the mean time
shift between the prediction and the target timeseries and is defined by

$C(\tau)\propto\langle y(t+\tau)\hat{y}(t)\rangle,$ (5)

provided that $y$ and $\hat{y}$ timeseries have been subtracted to their mean
values. We used the cross correlation function to measure the warning time
statistically. In particular, we defined the warning time $\Delta t$ as a
function of the lag $\bar{\tau}$ that maximizes the cross correlation
function, i.e., $\Delta t=t_{f}-\bar{\tau}$ is given by the lookforward of the
model (which is $t_{f}=1$ hour), minus the lag $\bar{\tau}$. Thus, by
definition, the warning time $\Delta t$ may be $\Delta t\leq 1$ hr.

To test the model we selected a set of testing samples. Each sample contains
the signature of a proton flux enhancement, not necessarily corresponding to
the overcoming of the 10 pfu threshold (defined for GOES proton flux at $>10$
MeV). Furthermore, the samples used for testing have been excluded from the
training dataset, in order to avoid any bias in the estimation of performance
metrics. The results are summarized in Table 3 where, from left to right, we
report: the start date, the end date, the coefficient of determination
$R^{2}$, the RMSE, the peak flux reached by the target time series in the
selected time interval, the peak reached by the predicted time series in the
same time interval, and the statistical warning time based on the estimation
of the cross correlation function. Generally, the model shows a quite good
performance in terms of $R^{2}$, with the exception of sample 31, for which
$R^{2}=-0.49$ is negative. However, the bottom panel of Figure 4, which
compares the model’s prediction with the target for sample 31, shows that even
though the flux is overestimated, the trend is quite well reproduced,
especially during the enhancement of the proton flux observed on 2016-07-23.
In terms of the RMSE, which quantifies the fluctuations of the prediction with
respect to the target, the model has better performance for higher values of
the proton flux (see e.g. samples 12 and 34). This may be due to the effect of
the log-transformation. Indeed, as can be also seen from Figure 3, the
fluctuations of the target are generally more pronounced for lower values.

The warning time $\Delta t$ is predominantly greater than zero, meaning that
the model forecast the enhancement of the proton flux in advance. Almost for
all the testing samples the warning time is equal to the lookforward of the
model (1 hr), i.e., predicted and target timeseries are synchronized on
average. The only exception is sample 24 which is not timely predicted.

Figure 4 shows the results in terms of the timeseries of three test samples
(e.g. from top to bottom, 23, 18 and 31). The left, top and middle panels show
two testing samples for which the model is quite good in terms of estimation
of the flux and timing of both the onset and the peak. The right top and
middle panels of each row also show the correlation plot between the
prediction and the target. As can be seen, the points are almost distributed
around the bisector, meaning that the target and the prediction are
comparable. In the bottom panel, however, we show the worse sample in terms of
model performance (sample 31). Even though the flux is here overestimated and
the points are not distributed along the bisector in the correlation plot
(right bottom panel of Figure 4), the enhancement is still predicted by the
model as can be seen also from the fact that the points are linearly related.

Figure 3: A 3-months sample dataset used for testing/visualization purposes. Green line represents the target of the model (true values), while blue line represents the output of the model (forecasting). Both are representative of the logarithm of the proton flux integrated from 7.8 MeV and 53 MeV energy channels (P8, P25, P41; see Table 1). The proton flux is given in pfu, while the logarithm is adimensional. The inset shows a closer look around the ESP event of 2012-03-08, which is not well reproduced by the model. This is due to the fact that the electrons do not exhibit the ESP signature (see the light blue line, which shows the track of the E300 electron differential channel). Sample N. | Start date | End date | $R^{2}$ | RMSE | target peak | pred. peak | Warning time [min]
---|---|---|---|---|---|---|---
1 | 2012-02-24 | 2012-03-01 | 0.92 | 0.33 | 0.96 | 1.28 | 60
2 | 2012-03-04 | 2012-03-22 | 0.97 | 0.26 | 3.89 | 3.22 | 60
3 | 2012-05-16 | 2012-05-24 | 0.93 | 0.34 | 1.55 | 1.89 | 60
4 | 2012-05-26 | 2012-06-01 | 0.94 | 0.25 | 1.21 | 0.80 | 60
5 | 2012-06-14 | 2012-06-21 | 0.71 | 0.56 | 1.77 | 0.34 | 60
6 | 2012-07-05 | 2012-08-05 | 0.81 | 0.37 | 2.25 | 1.69 | 20
7 | 2012-08-27 | 2012-09-18 | 0.89 | 0.41 | 1.97 | 1.58 | 60
8 | 2012-09-17 | 2012-10-13 | 0.88 | 0.36 | 1.25 | 1.35 | 60
9 | 2013-02-21 | 2013-03-02 | 0.69 | 0.46 | 0.13 | -0.13 | 60
10 | 2013-03-05 | 2013-04-01 | 0.93 | 0.28 | 1.54 | 1.05 | 60
11 | 2013-04-08 | 2013-05-01 | 0.91 | 0.34 | 2.00 | 1.58 | 60
12 | 2013-05-09 | 2013-06-05 | 0.94 | 0.39 | 3.19 | 3.14 | 60
13 | 2013-06-21 | 2013-07-03 | 0.93 | 0.32 | 1.27 | 0.89 | 60
14 | 2013-08-15 | 2013-09-08 | 0.87 | 0.26 | 0.35 | 0.24 | 60
15 | 2013-09-28 | 2013-10-10 | 0.96 | 0.31 | 2.55 | 2.23 | 60
16 | 2013-10-11 | 2013-11-25 | 0.83 | 0.34 | 0.72 | 1.40 | 60
17 | 2013-12-27 | 2013-12-31 | 0.48 | 0.44 | 1.43 | 0.52 | 60
18 | 2014-01-02 | 2014-01-22 | 0.96 | 0.30 | 3.22 | 3.09 | 60
19 | 2014-02-16 | 2014-02-24 | 0.80 | 0.40 | 1.43 | 1.53 | 30
20 | 2014-02-24 | 2014-03-21 | 0.97 | 0.26 | 2.18 | 2.24 | 60
21 | 2014-03-28 | 2014-04-02 | 0.84 | 0.28 | 0.45 | 1.06 | 59
22 | 2014-04-16 | 2014-04-30 | 0.90 | 0.43 | 2.17 | 1.66 | 60
23 | 2014-08-29 | 2014-09-30 | 0.94 | 0.26 | 2.78 | 1.97 | 60
24 | 2014-10-31 | 2014-11-08 | 0.90 | 0.36 | 1.28 | 0.77 | 0
25 | 2015-06-17 | 2015-07-08 | 0.83 | 0.58 | 3.45 | 2.40 | 60
26 | 2015-09-18 | 2015-09-26 | 0.89 | 0.25 | 0.32 | 0.66 | 60
27 | 2015-10-26 | 2015-11-04 | 0.91 | 0.32 | 1.00 | 1.56 | 60
28 | 2015-12-27 | 2015-12-31 | 0.68 | 0.60 | 0.76 | 0.17 | 60
29 | 2016-01-01 | 2016-01-07 | 0.78 | 0.56 | 1.80 | 0.70 | 27
30 | 2016-05-15 | 2016-06-19 | 0.56 | 0.32 | 0.41 | 0.07 | 60
31 | 2016-07-22 | 2016-07-25 | -0.49 | 0.73 | -0.32 | 1.49 | 41
32 | 2017-07-13 | 2017-07-19 | 0.91 | 0.41 | 1.42 | 1.09 | 60
33 | 2017-07-22 | 2017-08-06 | 0.88 | 0.27 | 0.15 | 0.12 | 60
34 | 2017-09-02 | 2017-09-30 | 0.90 | 0.48 | 3.30 | 3.17 | 32
35 | 2021-10-28∗ | 2021-11-05∗ | 0.84 | 0.45 | 0.87 | 1.08 | 58
36 | 2022-03-28∗ | 2022-04-07∗ | 0.82 | 0.46 | 0.86 | 1.50 | 48
Table 3: Error metrics for sample events compared with the peak flux. From
left to right: the start date, the end date, the coefficient of determination
$R^{2}$, the RMSE, the target peak, the predicted peak and the warning time.
Note that the target peak, the predicted peak and the RMSE are referred to the
logarithm of the proton flux. Thus, given for example the predicted peak $x$,
the proton flux in pfu is $10^{x}$. The warning time $\Delta t$ is computed as
the difference between the lookforward and the cross correlation lag, i.e.,
$\Delta t=t_{f}-\bar{\tau}$. We also included two events after 2017 October 4
(denoted by asterisks). The results for these two samples have to be taken
with care.

Figure 4: (Left panels) Comparison between the predicted time series and the
target one for three samples: 23, 18 and 31 (from top to bottom). (Right
panel) correlation plot between the predicted and the target points. The grey
line represents the bisector (not the linear regression line). Note that the
sample 31 is the worst possible example (see Table 3). However, even though
the model overestimates the flux, the trend is reproduced quite well during
the enhancement of 2016-07-23.

##  5 Discussion and conclusion

In this work we have presented a pipeline for predicting the energetic proton
flux by using a machine learning model based on the RF algorithm. We
considered data from the SOHO/EPHIN telescope for training and testing the
model. The reason why we used the SOHO/EPHIN data for demonstrating the
effectiveness of our approach is that differential channels for both protons
and electrons are embedded in the same dataset, so there are no possible
effects due to spacecraft placed in different environments. The pipeline for
generating the most informative dataset resulted in a set of 27 features, all
extracted from the electron flux channels and derivative of the integrated
proton flux channel, both with different time lags. The model average
performance is reported in Table 4. All the metrics reported suggest that the
model has a remarkable performance in predicting the energetic proton flux up
to one hour ahead by exploiting information derived from the electron
differential channels. Nevertheless, we remark that these are the average
characteristics of the model. There may be some events anticipating or
delaying the characteristics of the time profile, making a non-timely
prediction of the onset or of the peak (see Table 3). This may be due to the
fact that the delay between protons and electrons (onset or peak) changes
dynamically event per event, depending on several factors, such as: the
magnetic connectivity with the source of particles, the different release time
and/or sources of electrons and protons, and their propagation in the
interplanetary space (Kollhoff et al., 2021; Dresing et al., 2023). According
to Posner (2007) the distribution of the delay between protons and electrons
onsets range from 20-30 minutes up to 1 hour. However, in this work, we took
into account the dependence of the delay on the magnetic connectivity by
employing the time derivative of both the electron and proton differential
channels, which in turn contain the information about the rise parameter. The
resulting mean warning time (averaged over the samples reported in Table 3) is
54 minutes.

At its current state, the model produces a single output $y_{t+t_{f}}$. In
principle, a multi-output regressor can represent a natural evolution of the
RF-based model presented in this work, providing a vector
$(y_{t+1},y_{t+2},...,y_{t+t_{f}})$ which contains the prediction of the time
profile in the interval $[t+1,t+t_{f}]$, where $t_{f}$ indicates the size of
the forecasting window (i.e. the lookforward). This perspective could be
particularly helpful for having an idea of SPEs duration and generally more
accurate information on the timings.

From a physical perspective, the predictability of the proton flux from the
electron flux is plausible, since protons and electrons propagation is along
the same interplanetary magnetic field configuration (a spiral during non
disturbed periods), given that they are accelerated during the same solar
eruption, with electrons arriving earlier than protons at the observer
location due to the mass difference. Hence, electrons may carry information
about the value of proton flux and the magnetic connectivity of the observer
with the source (Posner, 2007). On the contrary, characteristics of the proton
flux associated with local proton/ion energization processes, such as ESPs
usually observed during interplanetary shock passages, are not present in the
electron flux channels and then are not reproduced by the model.

$R^{2}$ | RMSE | Warning time [min]
---|---|---
0.82 | 0.38 | 54
Table 4: Error metrics for the model averaged over the test samples reported
in Table 4. From left to right: coefficient of determination $R^{2}$, root
mean squared error and the warning time. Note that the RMSE is referred to the
logarithm of the proton flux.

In this work we used specific differential channels for both electrons and
protons, therefore our model is very specialized to predict the protons in the
EPHIN integral channel. However, these results highlight that a very simple
algorithm such as the RF is able to capture the variability of the solar
proton flux by extracting basic features (mean, standard deviation and
derivative) from the electron channels. Thus we may speculate that this
approach can possibly be valid also for other instruments offering similar
observations. In those cases, possible differences in the performance would
depend mainly on the quality of the available data. This makes our approach
very promising in view of space missions carrying on board particle detectors
that allows for electron and proton flux measurements, with the option of
providing a dedicated SPEs prediction system to be directly implemented on
board. In principle, the model can be coupled with a decision rule to issue an
alert. For example the most naive approach would be to issue an alert if the
predicted value overcome a threshold $\epsilon$. Another approach could be to
issue an alert when the predicted time series presents an anomalous value.
More sophisticated approach are also feasible in terms of issuing alerts, like
anomaly detection algorithm running directly on the predicted timeseries.
However, the advantage of predicting the future values of the proton flux
instead of giving the classification as output is that the decision rule can
be adapted to any situation and also allows more detailed analysis of the
trend.

Therefore, this work is particularly relevant for the implementation of
Machine Learning algorithms for the next generation applications on board
space missions, which is a critical need for the future of space exploration.
For instance, the development of AI-based algorithms for the automation of on
board operations, including SPE forecasting, is a main goal of the ASAP
(Automatics in SpAce exPloration; https://asap-space.eu) project, funded by
the Horizon Europe Framework Programme. Our system can be also easily adapted
aboard the HENON (HEliospheric pioNeer for sOlar and interplanetary threats
defeNce) pathfinder mission dedidated to Space Weather operations and science.
HENON is currently in the detailed design Phase C1, foreseeing a cubesat on
the Distant Retrograde Orbit in deep space well upstream of L1 around 0.1 AU
toward the Sun. Being equipped with a radiation monitor, a Faraday cup and a
magnetometer, HENON allows, in principle, the early computation and near real
time transmission to Earth of alerts for the Space Weather events, including
SPEs (Laurenza et al., 2023; Provinciali et al., 2024). In addition, an
onboard SPE forecasting model is essential to address the issue of astronaut
radiation safety and to limit human exposure to high‐dose rates for space
exploration missions to the Moon and Mars. Possible forecasting systems
relying on remote sensing observations could face the problem at Mars distance
that for up to half of the relevant SPE events, the nominal magnetic ﬁeld line
connection to the Sun is at W90, implying that the full flare signature may
not be visible and obstructed from view. Forecasting systems based only on
near-relativistic observations, such as a two-element REleASE warning system
(Posner & Strauss, 2020) and the model presented here, can provide a valid
solution as they are not affected by this issue.

## Acknowledgements

This study has been performed in the framework of the HEliospheric pioNeer for
sOlar and interplanetary threats defeNce (HENON) mission Phase A/B and C1.
HENON is part of the Italian Space Agency (ASI) program Alcor and is being
developed under the European Space Agency General Support Technology Programme
(ESA-GSTP) through the support of the national delegations of Italy (ASI), UK,
Finland and Czech Republic. MS acknowledges the support from the ASI-SERENA
contract no. 2018-8-HH.O “Partecipazione scientifica alla missione BEPICOLOMBO
SERENA Fase E” and the ESA contract (RFP/NC/ IPL- PSS/JD/258.2016) Expert
Support to SERENA Science Operations. MS, ML, SB and MFM acknowledge the
SOHO/EPHIN team for providing the data used in this work.

## References

  * Alberti et al. (2017) Alberti, T., Laurenza, M., Cliver, E., et al. 2017, The Astrophysical Journal, 838, 59, doi: 10.3847/1538-4357/aa5cb8
  * Anastasiadis et al. (2017) Anastasiadis, A., Papaioannou, A., Sandberg, I., et al. 2017, Solar Physics, 292, 1, doi: 10.1007/s11207-017-1163-7
  * Arena et al. (2014) Arena, C., De Micco, V., Macaeva, E., & Quintens, R. 2014, Acta Astronautica, 104, 419
  * Arge et al. (2004) Arge, C. N., Luhmann, J. G., Odstrcil, D., Schrijver, C. J., & Li, Y. 2004, Journal of Atmospheric and Solar-Terrestrial Physics, 66, 1295, doi: 10.1016/j.jastp.2004.03.018
  * Aschwanden (2002) Aschwanden, M. J. 2002, Particle Acceleration and Kinematics in Solar Flares (Dordrecht: Springer Netherlands), 1–227, doi: 10.1007/978-94-017-2541-5_1
  * Balch (2008) Balch, C. C. 2008, Space Weather, 6, S01001, doi: 10.1029/2007SW000337
  * Boubrahimi et al. (2017) Boubrahimi, S. F., Aydin, B., Martens, P., & Angryk, R. 2017, in 2017 IEEE International Conference on Big Data (Big Data), IEEE, 2533–2542
  * Cane et al. (2003) Cane, H. V., von Rosenvinge, T. T., Cohen, C. M. S., & Mewaldt, R. A. 2003, Geophys. Res. Lett., 30, 8017, doi: 10.1029/2002GL016580
  * Carnell (2021) Carnell, L. S. 2021, International journal of radiation biology, 97, S125
  * Chang et al. (2016) Chang, P. Y., Cucinotta, F. A., Bjornstad, K. A., et al. 2016, Radiation research, 185, 449
  * Cliver (2008) Cliver, E. W. 2008, Proceedings of the International Astronomical Union, 4, 401–412, doi: 10.1017/S1743921309029639
  * Cohen et al. (1999) Cohen, C. M. S., Mewaldt, R. A., Leske, R. A., et al. 1999, Geophys. Res. Lett., 26, 2697, doi: 10.1029/1999GL900560
  * Cucinotta et al. (2001) Cucinotta, F., Manuel, F., Jones, J., et al. 2001, Radiation research, 156, 460
  * Cucinotta (2015) Cucinotta, F. A. 2015, Health physics, 108, 131
  * Desai & Giacalone (2016) Desai, M., & Giacalone, J. 2016, Living Reviews in Solar Physics, 13, 3, doi: 10.1007/s41116-016-0002-5
  * Dresing et al. (2023) Dresing, N., Rodríguez-García, L., Jebaraj, I. C., et al. 2023, A&A, 674, A105, doi: 10.1051/0004-6361/202345938
  * Eastwood et al. (2017) Eastwood, J., Biffis, E., Hapgood, M., et al. 2017, Risk analysis, 37, 206
  * Engell et al. (2017) Engell, A. J., Falconer, D. A., Schuh, M., Loomis, J., & Bissett, D. 2017, Space Weather, 15, 1321, doi: 10.1002/2017SW001660
  * Ferl et al. (2002) Ferl, R., Wheeler, R., Levine, H. G., & Paul, A.-L. 2002, Current opinion in plant biology, 5, 258
  * Fogtman et al. (2023) Fogtman, A., Baatout, S., Baselet, B., et al. 2023, npj Microgravity, 9, 8
  * Forbush (1946) Forbush, S. E. 1946, Physical Review, 70, 771
  * Furukawa et al. (2020) Furukawa, S., Nagamatsu, A., Nenoi, M., et al. 2020, BioMed Research International, 2020, 4703286
  * Hosseinzadeh et al. (2024) Hosseinzadeh, P., Boubrahimi, S. F., & Hamdi, S. M. 2024, The Astrophysical Journal Supplement Series, 270, 31
  * Hughson et al. (2018) Hughson, R. L., Helm, A., & Durante, M. 2018, Nature Reviews Cardiology, 15, 167
  * Iucci et al. (2005) Iucci, N., Levitin, A. E., Belov, A. V., et al. 2005, Space Weather, 3, 01001, doi: 10.1029/2003SW000056
  * Iucci et al. (2005) Iucci, N., Levitin, A., Belov, A., et al. 2005, Space weather, 3
  * Jandial et al. (2018) Jandial, R., Hoshide, R., Waters, J. D., & Limoli, C. L. 2018, Surgical neurology international, 9
  * Jones et al. (2005) Jones, J. B. L., Bentley, R., Hunter, R., et al. 2005, Advances in Space Research, 36, 2258
  * Kahler et al. (2007) Kahler, S., Cliver, E., & Ling, A. 2007, Journal of atmospheric and solar-terrestrial physics, 69, 43, doi: 10.1016/j.jastp.2006.06.009
  * Kahler et al. (1984) Kahler, S. W., Sheeley, N. R., J., Howard, R. A., et al. 1984, Journal of Geophysical Research, 89, 9683, doi: 10.1029/JA089iA11p09683
  * Kocharov & Torsti (2002) Kocharov, L., & Torsti, J. 2002, Sol. Phys., 207, 149, doi: 10.1023/A:1015540311183
  * Kollhoff et al. (2021) Kollhoff, A., Kouloumvakos, A., Lario, D., et al. 2021, A&A, 656, A20, doi: 10.1051/0004-6361/202140937
  * Laurenza et al. (2018) Laurenza, M., Alberti, T., & Cliver, E. 2018, The Astrophysical Journal, 857, 107, doi: 10.3847/1538-4357/aab712
  * Laurenza et al. (2009) Laurenza, M., Cliver, E., Hewitt, J., et al. 2009, Space Weather, 7, doi: 10.1029/2007SW000379
  * Laurenza et al. (2024) Laurenza, M., Stumpo, M., Zucca, P., et al. 2024, Journal of Space Weather and Space Climate, 14, 8, doi: 10.1051/swsc/2024007
  * Laurenza et al. (2023) Laurenza, M., Marcucci, M. F., Zimbardo, G., et al. 2023, in EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, EGU–14274, doi: 10.5194/egusphere-egu23-14274
  * Lavasa et al. (2021) Lavasa, E., Giannopoulos, G., Papaioannou, A., et al. 2021, Solar Physics, 296, 107, doi: 10.1007/s11207-021-01837-x
  * Luhmann et al. (2007) Luhmann, J. G., Ledvina, S. A., Krauss-Varban, D., Odstrcil, D., & Riley, P. 2007, Advances in Space Research, 40, 295, doi: 10.1016/j.asr.2007.03.089
  * Luhmann et al. (2010) Luhmann, J. G., Ledvina, S. A., Odstrcil, D., et al. 2010, Advances in Space Research, 46, 1, doi: 10.1016/j.asr.2010.03.011
  * Malandraki & Crosby (2018) Malandraki, O. E., & Crosby, N. B. 2018, Space Weather, 16, 591, doi: https://doi.org/10.1029/2018SW001950
  * Müller-Mellin et al. (1995) Müller-Mellin, R., Kunow, H., Fleißner, V., et al. 1995, Solar Physics, 162, 483, doi: 10.1007/BF00733437
  * Müller-Mellin et al. (1995) Müller-Mellin, R., Kunow, H., Fleißner, V., et al. 1995, The SOHO mission, 483
  * Nedal et al. (2023) Nedal, M., Kozarev, K., Arsenov, N., & Zhang, P. 2023, Journal of Space Weather and Space Climate, 13, 26, doi: 10.1051/swsc/2023026
  * Núñez (2011) Núñez, M. 2011, Space Weather, 9, S07003, doi: 10.1029/2010SW000640
  * Núñez (2015) —. 2015, Space Weather, 13, 807, doi: 10.1002/2015SW001256
  * Núñez & Paul-Pena (2020) Núñez, M., & Paul-Pena, D. 2020, Universe, 6, 161, doi: 10.3390/universe6100161
  * Odstrcil (2003) Odstrcil, D. 2003, Advances in Space Research, 32, 497, doi: 10.1016/S0273-1177(03)00332-6
  * Palmerio et al. (2022) Palmerio, E., Lee, C. O., Mays, M. L., et al. 2022, Space Weather, 20, e2021SW002993, doi: 10.1029/2021SW002993
  * Papaioannou et al. (2018) Papaioannou, A., Anastasiadis, A., Sandberg, I., & Jiggens, P. 2018, Journal of Space Weather and Space Climate, 8, A37, doi: 10.1051/swsc/2018024
  * Papaioannou et al. (2016) Papaioannou, A., Sandberg, I., Anastasiadis, A., et al. 2016, Journal of Space Weather and Space Climate, 6, A42, doi: 10.1051/swsc/2016035
  * Papaioannou et al. (2022) Papaioannou, A., Vainio, R., Raukunen, O., et al. 2022, Journal of Space Weather and Space Climate, 12, 24, doi: 10.1051/swsc/2022019
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Posner (2007) Posner, A. 2007, Space Weather, 5, doi: 10.1029/2006SW000268
  * Posner & Strauss (2020) Posner, A., & Strauss, R. 2020, Space Weather, 18, e2019SW002354, doi: 10.1029/2019SW002354
  * Provinciali et al. (2024) Provinciali, L., Calcagno, D., Amabili, P., et al. 2024, in 2024 IEEE Aerospace Conference, 1–12, doi: 10.1109/AERO58975.2024.10521299
  * Reames (2000) Reames, D. V. 2000, in American Institute of Physics Conference Series, Vol. 516, 26th International Cosmic Ray Conference, ICRC XXVI, ed. B. L. Dingus, D. B. Kieda, & M. H. Salamon (AIP), 289–300, doi: 10.1063/1.1291483
  * Riebesell & Bringuier (2023) Riebesell, J., & Bringuier, S. 2023, tikz: A Collection of TikZ Diagrams, GitHub.  https://github.com/janosh/tikz
  * Rotti et al. (2024) Rotti, S. A., Aydin, B., & Martens, P. C. 2024, The Astrophysical Journal, 966, 165
  * Sabbatini & Grimani (2022) Sabbatini, F., & Grimani, C. 2022, Aeronautics and Aerospace Open Access Journal, 6, 90, doi: 10.15406/aaoaj.2022.06.00145
  * Sabbatini et al. (2024) Sabbatini, F., Grimani, C., & Calegari, R. 2024, Advances in Space Research, 74, 505, doi: 10.1016/j.asr.2024.04.041
  * Shea & Smart (1990) Shea, M. A., & Smart, D. F. 1990, Solar Physics, 127, 297, doi: 10.1007/BF00152170
  * Smart & Shea (1989) Smart, D., & Shea, M. 1989, Advances in Space Research, 9, 281, doi: https://doi.org/10.1016/0273-1177(89)90450-X
  * Smart & Shea (1979) Smart, D. F., & Shea, M. A. 1979, in NOAA Solar-Terrestrial Predictions Proceedings. Volume 1., ed. R. F. Donnelly, Vol. 1, 406–427
  * Sridharan et al. (2015) Sridharan, D., Asaithamby, A., Bailey, S., et al. 2015, Radiation research, 183, 1
  * St. Cyr et al. (2017) St. Cyr, O., Posner, A., & Burkepile, J. 2017, Space Weather, 15, 240, doi: 10.1002/2016SW001545
  * Stumpo et al. (2021) Stumpo, M., Benella, S., Laurenza, M., et al. 2021, Space Weather, 19, e2021SW002794, doi: 10.1029/2021SW002794
  * Whitman et al. (2023) Whitman, K., Egeland, R., Richardson, I. G., et al. 2023, Advances in Space Research, 72, 5161, doi: 10.1016/j.asr.2022.08.006
  * Wu et al. (2009) Wu, H., Huff, J. L., Casey, R., Kim, M.-H., & Cucinotta, F. A. 2009, The Human Health and Performance Risks for Space Explorations. Houston, Texas: NASA Human Research Program, 171
  * Zhang et al. (2023) Zhang, J., Guo, J., & Dobynde, M. I. 2023, Space Weather, 21, e2023SW003490
  * Zhang et al. (2021) Zhang, J., Temmer, M., Gopalswamy, N., et al. 2021, Progress in Earth and Planetary Science, 8, 56, doi: 10.1186/s40645-021-00426-7