Extended pan-spectrum fitting of energetic particle spectra

1.   Introduction

Energetic particles are commonly observed in space, including solar energetic particles (SEPs) (Ellison and Ramaty, 1985; Lin RP, 1985; Lin RP and Schwartz, 1987; Mason et al., 2002; Mewaldt et al., 2005; Wang LH et al., 2006, 2022; Krucker et al., 2009; Wang XD et al., 2022; Wei WW et al., 2022), solar wind suprathermal particles (Lin RP, 1998; Gloeckler, 2003; Mason and Gloeckler, 2012; Yoon et al., 2012; Wang LH et al., 2015, 2022; Tao JW et al., 2016, 2021), galactic cosmic rays (Cummings et al., 2016), anomalous cosmic rays (ACRs) (Senanayake et al., 2015), energetic particles in magnetosphere (Zong Q-G et al., 2017; Ye YG et al., 2021; Chen JL et al., 2024), etc. The space-borne and ground observations suggest that, in a wide energy range, these energy particle populations show a non-thermal spectrum of particle intensity versus energy, often with one or two spectral transitions/breaks (Lin RP and Schwartz, 1987; Krucker et al., 2009; Wang LH, 2022). For example, Lin RP (1985) reported that the SEP electrons exhibited a power-law flux energy spectrum bending upward at energies above ~10 keV, then turning downward at energies above ~100 keV. The observed energy spectral shape can reflect the physical nature of these energetic particles’ origins, acceleration, and transportation processes.

It is customary to characterize the spectral features by fitting the observed particle energy spectrum to a parameterized formula. However, many previous studies obtained the spectral features by fitting the energetic particle observation to a function with an assumed spectral shape involving no more than one spectral transition. For instance, the SEP electrons (Wang LH et al., 2006; Krucker et al., 2009) were often fitted to a single-power-law (SPL) energy spectrum expressed as $J\left(E\right)=A\times {E}^{-\beta }$, or a classical double-power-law (DPL) spectrum expressed as

where E is the particle energy, J is the particle intensity (flux or fluence), A is the amplitude, ${E}_{0}$ is the break energy, and $\beta$, ${\beta }_{1}$ or ${\beta }_{2}$ is the power-law spectral index. The SEP protons (Ellison and Ramaty, 1985; Mewaldt et al., 2005) were typically fitted to an Ellison$-$Ramaty (ER) spectrum as

where ${E}_{c}$ is the cutoff energy. In addition, the SEP ³He-rich ions often show a power-law or curved energy spectrum (Mason et al., 2002). On the other hand, the solar wind suprathermal electrons (Tao JW et al., 2016, 2021; Wang LH, 2022) consist of halo/strahl populations at ~0.1$-$2 keV that were generally fitted to a Kappa distribution, expressed as

with the Kappa index $\kappa$, the number density ${n}_{0}$ and the efficient temperature ${W}_{0}$, as well as a superhalo population at energies above 2 keV that is fitted well to an SPL energy spectrum.

Recently, Liu ZX et al. (2020) constructed a pan-spectrum fitting formula with continuous energy derivatives of J at any orders, which incorporates a generalized DPL function, SPL function, ER function, Kappa function, Maxwellian function, logarithmic-parabola (LP) function, etc. Therefore, the pan-spectrum formula can be utilized for self-consistently determining the spectral shape of energetic/suprathermal particles with no more than one spectral transition/break, as well as characterizing the sharpness/width of a spectral transition. Furthermore, they improved the spectrum fitting method by considering the uncertainties in both particle intensity and energy.

In this paper, we further propose an extended pan-spectrum fitting formula with eight parameters, in order to self-consistently determine the spectral shape of energetic particles over a wide energy range with two spectral transitions/breaks. Moreover, we develop a fitting method that can help to comparatively investigate the spectral properties between various energetic particle phenomena.

2.   Extended Pan-Spectrum (EPS) Formula

2.1   Construction of EPS Formula

To fit the single-break energy spectra, Liu ZX et al. (2020) proposed a pan-spectrum formula,

integrated from

with continuous derivatives at any orders. It’s noted that the pan-spectrum formula was so named because it can incorporate many previously used spectrum formulae, typically with a single break energy, including classical DPL, ER, Kappa, Maxwellian function, etc., as special cases.

In order to achieve a full-range energy spectrum fitting for the common situations of one or two break energies, we extend the formalism of pan-spectrum to a formula named the “Extended pan-spectrum (EPS) Formula”, that can fit a spectral shape with 0, 1, or 2 spectral transitions/breaks. The EPS formula simultaneously meets the following requirements: (1) it can be divided into three segments clearly; (2) it can merge the three segments into one or two parts. We utilize two smooth tanh functions at break energies:

Based on the property that when $x\to -\infty \left(+\infty \right),\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(x\right)\to -1(+1)$, it satisfies the piecewise requirement:

After integrating Equation (3) over the energy E, we obtain the formula of $J\left(E\right)$ as:

where A is the amplitude coefficient, ${\alpha }_{1}$ and ${\alpha }_{2}$ ($> 0$) determines the sharpness and width of energy transitions centering at ${E}_{1}$ and ${E}_{2}$, ${\beta }_{1}$, ${\beta }_{2}$, and ${\beta }_{3}$, respectively indicate the spectral indexes of three segments. As defined by Liu ZX et al. (2020), the spectral transition/break becomes smoother with increased width as the parameter $\alpha$ (>0) decreases. Figure 1 compares the sharpness/width of spectral transition when ${\alpha }_{1}=2$ (top panels) and ${\alpha }_{2}=10$ (bottom panels). In addition, the all-order derivatives of the EPS formula are continuous, with the 1^st-order derivative being the combination of tanh functions and the 2^nd-order derivative containing 2 dips.

Figure  1.  Curves of EPS formula (a, d), and its first-order (b, e) and second-order (c, f) derivative. The dashed lines indicate the central transition, and the double-ended arrows denote the transition widths. The top (bottom) panels show a smooth (sharp) transition with ${\alpha }_{1}=2\;({\alpha }_{1}=10)$, when ${\beta }_{1}=1.5,{\beta }_{2}=6,{\beta }_{3}=1.5,{E}_{1}=1.8\;\mathrm{k}\mathrm{e}\mathrm{V}$, and ${E}_{2}=60\;\mathrm{k}\mathrm{e}\mathrm{V}$.

DownLoad: Full-Size Img PowerPoint

2.2   Different Limiting Approximations

2.2.1   Three segments

(1) Generalized triple-power-law (TPL)

The formula can be rewritten as a continuous function (Wang W et al., 2023):

where the three segments are connected by two smooth transitions (when ${\alpha }_{1}$ and ${\alpha }_{2}$ are finite). The transition energies here are defined as

When both ${\alpha }_{1}$ and ${\alpha }_{2}$ approach infinity, the formula is equivalent to a strict triple-power-law (TPL), which is simply described as:

(2) Kappa+SPL

When ${\beta }_{1}=-1$ and ${\alpha }_{1}=1$, the formula is written as a continuous function:

At low energies, $J\left(E\right)$ is a kappa-like distribution and the dominated part is $E{\left(1+\dfrac{E}{{E}_{1}}\right)}^{-\left({\beta }_{2}+1\right)}$, where the kappa index $\kappa$ is equivalent to ${\beta }_{2}+1$ and the efficient temperature is described as ${W}_{0}=\dfrac{{E}_{1}}{{\beta }_{2}}$. The middle energies are in transition between Kappa-like distribution and power-law distribution, shown by the second part of function (8). At high energies, $J\left(E\right)$ is a power-law distribution, the spectral index of which is $-{\beta }_{3}$ times a distribution ${E}^{\left(1+{\beta }_{2}\right)}{\left(1+\dfrac{E}{{E}_{1}}\right)}^{-\left({\beta }_{2}+1\right)}$, that tends to approach a constant when $E\gg {E}_{1}$. Therefore, when ${\alpha }_{2}$ approaches infinity and ${E}_{2}\gg {E}_{1}$, the limitation of the formula is:

(3) SPL+ER

When${\beta }_{3}\to \infty$ and ${E}_{2}\to \infty$, the formula is written as a continuous function:

where ${E}_{c}$, defined by ${E}_{2},{\alpha }_{2},{\beta }_{2}$, and ${\beta }_{3}$, represents the cut-off energy that divides the ER function into a power-law below and an exponential cut-off form above. This formula can be written in three parts, containing a power-law dominated part, a transition part, as well as an ER-dominated part:

when $E < {E}_{1}\ll {E}_{c}$, we can obtain ${\mathrm{e}}^{-{\left(\frac{E}{{E}_{c}}\right)}^{{\alpha }_{2}}}\to 1$, and when ${\alpha }_{1}\to \infty$, we can obtain the strict SPL + ER form:

(4) LP+LP

When ${\alpha }_{2}\to 0$, via Taylor expression of $\mathrm{l}\mathrm{n}J$, the formula is written as a continuous function:

For simplicity, we use

to refer to the quadratic term of Taylor expression. In Equation (13),

is rewritten as ${\mathrm{e}}^{\varphi \left({\alpha }_{2},E\right)}$, equivalent to $\mathrm{l}\mathrm{n}J=\varphi \left({\alpha }_{2},E\right)$. Given that the logarithmic parabola term acts on the entire energy range $\left[0, +{\infty }\right]$, it should not be neglected when rewriting Equation (13) to three parts:

When ${\alpha }_{1}\to \infty$, we can obtain a simplified form like:

that is a combination of two logarithmic parabola at energies below or above ${E}_{1}$, accompanied with different energy axes of symmetry.

2.2.2   One or two segments

On one hand, when ${\beta }_{2}={\beta }_{3}$, the formula degenerates to pan-spectrum form (Liu ZX et al., 2020),

which can fit well to a one- or two-segment spectrum. On the other hand, when ${\beta }_{2}\ne {\beta }_{3}$, the formula can also degenerate to fewer segments. For instance, when ${\alpha }_{1}\to 0$ and ${\alpha }_{2}\to 0$, the EPS formula is simplified to a logarithmic parabola:

or the form:

with an infinite transition energy interval of [0, $+{\infty }$]. Note that the logarithmic parabola formula can't be characterized by fitting parameters themselves, and we usually use the slope at ${E}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and ${E}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ or the energy axis of symmetry to depict the spectrum shape.

To summarize, the EPS formula can cover several formulae mentioned above (Equations (16)−(17)), under different limiting approximations, representing a powerful and versatile, generalized formula for interplanetary energetic spectra.

3.   Methods

3.1   Fitting Method

Suggested by Liu ZX et al. (2020), we utilize the non-linear least-square algorithms to fit a series of n data points $\left({x}_{i},{y}_{i}\right)$ to a parameterized function, $y=f\left(x,\boldsymbol {a}\right)$, where $\boldsymbol {a}=\left({a}_{1},{a}_{2},\dots ,{a}_{m}\right)$ is a vector of m undetermined parameters. By minimizing the reduced chi-square statistic (e.g., Bevington and Robinson, 2003; Liu ZX et al., 2020), the best fit is obtained.

In this study, for a serious of n data points (${E}_{i},{J}_{i}$), where E refers to the central energy of every energy channel and J refers to the flux at every energy channel, we fit the data points to our proposed EPS formula, $J=f\left(E,\boldsymbol {a}\right)$, where $\boldsymbol {a}=\left({A}_{0},{\beta }_{1},{\beta }_{2},{\beta }_{3},{\mathrm{l}\mathrm{n}E}_{1},{\mathrm{l}\mathrm{n}E}_{2},{\alpha }_{1},{\alpha }_{2}\right)$ is an 8-dimensional vector. The reduced chi-square statistic is used to measure the error between the fitted and observed values. It is the sum of normalized residuals. Considering the uncertainties in both energies and fluxes suggested by Liu ZX et al. (2020), the reduced chi-square statistic is defined as,

where ${{\sigma }}_{{E}_{i}}^{2}$ (${{\sigma }}_{{J}_{i}}^{2}$) is the standard deviation of observed ${E}_{i}$ (${J}_{i}$), m is the number of parameters, n is the number of data points, n−m is the statistical degree of freedom, and ${f}_{{E}_{i}}^{\;'2}=\dfrac{\partial f}{\partial E}{|}_{E={E}_{i}}$.

We utilize the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Broyden, 1970; Fletcher, 1987) to find the minimum of ${\chi }_{\nu }^{2}\left(\boldsymbol {a}\right)$ in the multi-dimensional parameter space, which is an iterative method for solving unconstrained nonlinear optimization problems. The optimization problem is to minimize ${\chi }_{{\nu }}^{2}\left(\boldsymbol {a}\right)$. It proceeds according to these steps:

Step (1) starts with a given initial estimate ${\boldsymbol {a}}^{0}$, a threshold $\epsilon$ where the iteration stops, and an iteration maximum number.

Step (2) calculates the gradient $\nabla f\left({\boldsymbol {a}}_{k}\right)$. If $\|\nabla f\left({\boldsymbol {a}}_{k}\right)\| < \epsilon$, the iteration ends. Otherwise, the iteration moves on. (Note that $f\left({\boldsymbol {a}}_{k}\right)$ here is the ${\chi }_{\nu }^{2}$($\boldsymbol {a}$)).

Step (3) obtains the searching direction by solving the quasi-Newton equation,

where ${{\boldsymbol{B}}}_{k}$ is an approximation of the Hessian matrix at ${\boldsymbol {a}}_{k}$, which is updated iteratively at each stage, and $\nabla f\left({\boldsymbol {a}}_{k}\right)$ is the gradient of the function evaluated at ${\boldsymbol {a}}_{k}$.

Step (4) performs a one-dimensional optimization (line search) to find an acceptable step size α_k in the direction obtained in Step (3). α_k = $\mathrm{arg}\mathrm{min}\;f\left({\boldsymbol {a}}_{k}+\alpha {\boldsymbol {p}}_{k}\right)$, which means that α_k minimizes the function $f\left({\boldsymbol {a}}_{k}+\alpha {\boldsymbol {p}}_{k}\right)$.

Step (5) sets ${\boldsymbol {s}}_{k}={\boldsymbol {\alpha }}_{k}{\boldsymbol {p}}_{k}$ and updates ${\boldsymbol {a}}_{k+1}={\boldsymbol {a}}_{k}+{\boldsymbol {s}}_{k}$.

Step (6) calculates

Step (7) takes ${\boldsymbol {a}}_{k+1}={\boldsymbol {a}}_{k}$ and iteratively moves to Step (2).

Finally, the best fit ${\boldsymbol {a}}^{{\mathrm{fit}}}$ is obtained. As for estimating the uncertainties of fitted parameters, as suggested by Liu ZX et al. (2020) and Press (2007), we derive the covariance matrix C of $\boldsymbol {a}$ from the Hessian matrix of

and

The fitted ${\boldsymbol {a}}^{{\mathrm{fit}}}$ is the diagonal vector of matrix C, and the standard deviation of fitted parameters is defined as

Besides using ${\chi }_{{\nu }}^{2}$ to select the best results under a fitting model, in our fitting procedure, we also calculate

(Expected Cross-Validation Index) as a model selection index and a smaller ECVI indicates a better fit model (Liu ZX et al., 2020).

3.2   Parameter Analysis

We utilize parametric bootstrap (Efron, 1992; MacKenzie et al., 2017) to characterize the uncertainties of parameter estimation and test the correlations between parameters. Parametric bootstrap is based on the assumption that the initial data set is from a distribution of a specific parametric type. In this paper, the parametric model $J\left(E\right)$ is fitted by parameter vector $\boldsymbol {a}$, then bootstrapped random samples are drawn from this model based upon the initial data set. Repeating many times of fitting to the generated samples produces the probability distribution of fitted parameters.

In this study, given an initial parameter vector ${\boldsymbol {a}}=({\mathrm{e}}^{20},2.0, 1.0,3.5,1.1\;\mathrm{k}\mathrm{e}\mathrm{V},13.4\;\mathrm{k}\mathrm{e}\mathrm{V},5.0,3.0)$, we construct a group of 40 initial data points $\left({E}_{i},{J}_{i}\right)$. Then we generate 1000 groups of random data points by adding a relative Gaussian noise of 0.05 to each data point in both ${E}_{i}$ and ${J}_{i}$. Finally, we fit these 1000 groups of data to obtain a distribution of minimal chi-square and all fitted parameters. During the bootstrap simulation, around 180 samples are not converged and they are excluded in the statistical results that follow.

Figure 2 is the bootstrap simulation results for EPS formula with all eight parameters effective, one of which is closest to the initial parameters, shown in panel (b). The simulated ${\chi }_{{\nu }}^{2}$ probability distribution agrees with the ideal reduced chi-square distribution for 32 degrees of freedom in panel (a) (shown as dashed lines). Figure 2 (c−j) plot the histograms of the fitted parameters, indicating that all the spectral parameters are reasonably distributed, and in accordance with the initial parameters (shown as dashed lines).

Figure  2.  Histograms of fitted parameters and one fitted sample in bootstrap simulation. (a) The probability distribution function (PDF) of simulated minimal ${\chi }_{{\nu }}^{2}$ for the 1000 test samples, compared with the ideal reduced chi-square distribution (dashed line). (b) A sample in the simulation, with fitted parameters written in, central transitions shown in red dashed lines, and transition widths shown in double-ended arrows. (c$-$j) The PDF of the eight simulated parameters, with the simulation initial parameters shown in dashed lines and the simulated average values with standard deviations shown in green arrows.

DownLoad: Full-Size Img PowerPoint

In addition, Figure 3 plots the scatter diagrams between fitted parameters, also indicating that the averaged parameters (shown as green horizontal and vertical solid lines) agree with the initial parameters (shown as purple stars). We note: (1) the strong positive correlation between ${A}_{0}$ and ${\beta }_{1}$ with a correlation coefficient CC > 0.7; (2) the moderate positive correlation between ${\beta }_{2}$ and $\mathrm{l}\mathrm{n}{E}_{2}$ with 0.5 < CC < 0.7; (3) the moderate negative correlation between $\mathrm{l}\mathrm{n}{E}_{1}$ and $\mathrm{l}\mathrm{n}{E}_{2}$ with 0.5 < −CC < 0.7. These correlations are reasonable. For example, if ${A}_{0}$ increases, ${\beta }_{1}$ will also increase, leading to the formation of a steeper spectrum, thus reducing the residuals over the full energy range.

Figure  3.  Scatter plots between parameters in bootstrap simulation. The purple stars denote the initial parameters, and the green horizontal and vertical solid lines denote the standard deviations of the simulated parameters.

DownLoad: Full-Size Img PowerPoint

4.   Applications of EPS Formula

We apply the fitting method to the energy spectra of various representative particle phenomena that achieve the best and most reliable fit, with reasonable parameters, acceptable uncertainty estimations, as well as randomly distributed residuals. As for the fitting uncertainties, since the measurements given by different instruments behave differently, here we use both ${{\sigma }}_{{E}_{i}}=0.05{E}_{i}$ and ${{\sigma }}_{{J}_{i}}=0.05{J}_{i}$ for consistency.

Figure 4 and Table 1 show the fitted results of several typical particle phenomena. In Figure 4a, solar energetic electrons (SEEs) (Lin RP, 1985) are fitted to a strict TPL with two sharp transitions, that can be interpreted as a power-law of ${E}^{\;-{\beta }_{1}}$ at very low energies, ${E}^{\;-{\beta }_{2}}$ at middle energies that are far enough from both ${E}_{1}$ and ${E}_{2}$, and ${E}^{\;-{\beta }_{3}}$ at very high energies. By comparison, in Figure 4c, the O⁺ in an SEP event (Mewaldt et al., 2005) shows a generalized TPL, with two smooth transitions, that are 0.28−0.54 and 1.2−15.5 MeV/nuc. In Figure 4e, the halo and superhalo electrons (Yoon et al., 2012) in quiet solar wind are fitted to the combination (red line) of a kappa distribution at low energies (shown in purple text), and an SPL tail at high energies (shown in blue text). The fitting is performed at suprathermal energies with a flux intensity at least five times higher than the thermal Maxwellian distribution estimate (Guo XN et al., 2024). The close fit observed here validates our limit hypothesis. Furthermore, previous studies chose to fit the halo electrons and superhalo tail separately by dividing the energies (Yoon et al., 2012; Wang LH, 2022), with similar results to our monolithic fitting. In Figure 4g, the ions in a stream interaction region (SIR) event (Wei WW et al., 2022) seem to show a combined spectrum of SPL and ER-like tail at energies below/above 0.56 MeV/nuc. Figures 4i and k apply the method to other usual phenomena: the pick-up ions (PUI), the suprathermal tail (Gloeckler, 2003; Mason and Gloeckler, 2012), and the ACRs (Senanayake et al., 2015). Each are generally well-fitted to the EPS formula. Note that the EPS fitting method has already taken into account the effects of data from different instruments with different energy channel resolutions, by considering the energy bandwidth (energy error) of observations. In addition, the presence of data gaps could affect the uncertainty estimates (see Figure 4).

Figure  4.  Best fits and the normalized residuals of the representative energy spectra of SEE (a$-$b), SEP O⁺ (c$-$d), halo and super halo electrons in solar wind (e$-$f), SIR ions (g$-$h), pick-up ions and the suprathermal tail (i$-$j), and ACRs (k$-$l). The colored lines are explained by the text in each panel, either for different formula parts or for different ion species.

DownLoad: Full-Size Img PowerPoint

Table  1.  Fitted parameters of Figure 4.

		${\beta }_{1}$	${\beta }_{2}$	${\beta }_{3}$	${E}_{1}$	${E}_{2}$	${\alpha }_{1}$	${\alpha }_{2}$
SEE		2.8 ± 0.2	1.24 ± 0.04	2.72 ± 0.08	6.1 ± 0.6 keV	100 ± 8 keV	30a	25a
SEP O+		1.9 ± 0.3	0.3 ± 0.3	3.4 ± 0.3	0.39 ± 0.06 MeV/nuc	4.3 ± 0.6 MeV/nuc	6.2 ± 5.3	1.2 ± 0.5
Halo & superhalo solar wind		−1a	6.8 ± 0.4	2.74 ± 0.04	0.25 ± 0.03 keV	2.6±0.2 keV	1a	10 ± 5.5
SIR ions		4.2 ± 0.2	1.43 ± 0.09	infinity	0.56 ± 0.04 MeV/nuc	infinity	3.8 ± 1.1	2.3 ± 0.1
PUI & suprathermal tail	H+	0.1 ± 0.1	7.7 ± 0.4	2.37 ± 0.01	2.86 ± 0.04 keV	4.7 ± 0.1 keV/nuc	102a	57a
	He+	0.59 ± 0.08	6.8 ± 0.4	3.27 ± 0.07	2.66 ± 0.03 keV	4.5 ± 0.3 keV/nuc	56a	70a
	He2+	4.20 ± 0.07	2.24 ± 0.06	2.90 ± 0.07	4.7 ± 0.2 keV	12.3 ± 0.7 keV/nuc	89a	77a
ACR	He	2.08 ± 0.04	0.28 ± 0.05	2.8 ± 0.4	12.8 ± 0.4 keV/nuc	53 ± 3 keV/nuc	37a	27a
	O	1.9 ± 0.5	−0.9 ± 1.3	4.1 ± 0.7	4.2 ± 0.1 keV/nuc	11 ± 2 keV/nuc	20a	2.0 ± 1.0
Note: a Parameters with no uncertainties indicate the data points are fitted under limiting approximations.

 | Show Table

DownLoad: CSV

5.   Conclusion

In this study, we construct an EPS formula (Equation (5)) to self-consistently determine the spectral shape of energetic particles over a wide energy range that often show two spectral transitions/breaks. The EPS formula has eight parameters: the amplitude coefficient A; central energy E₁ and E₂ of low-energy and high-energy spectral transition; sharpness parameter ${\alpha }_{1}$and ${\alpha }_{2}$ of low-energy and high-energy transition; and power-law spectral index ${\beta }_{1}$, ${\beta }_{3}$, and ${\beta }_{2}$, respectively, at energies below the low-energy transition, above the high-energy transition, and between the two transitions.

This EPS formula can be rewritten as a piecewise form of generalized TPL. When ${\beta }_{2}={\beta }_{3}$, it degenerates to the PS formula proposed by Liu ZX et al. (2020), which can self-consistently determine the spectral shape with one or no spectral transitions/breaks (e.g., the SPL, DPL, EP, Kappa, and LP functions). Furthermore, the EPS formula can be transformed under certain specific parameter conditions to spectral forms with two spectral transitions/breaks. For instance, a classical TPL with sharp spectral breaks when ${\alpha }_{1}\left({\alpha }_{2}\right)\to {\infty }$, a Kappa distribution (at low energies) plus an SPL (at high energies) when ${\beta }_{1}=-1,{\alpha }_{1}=1,{\alpha }_{2}\to {\infty },\mathrm{a}\mathrm{n}\mathrm{d}\;E_{1}\ll {E}_{2}$, an SPL (at low energies) plus an ER-like distribution (at high energies) when ${E}_{2}\to {\infty },{\beta }_{3}\to {\infty },\mathrm{a}\mathrm{n}\mathrm{d}\;{E}_{1}\ll {E}_{c}$, and an LP plus LP distribution when ${\alpha }_{1}\to {\infty }\;\mathrm{a}\mathrm{n}\mathrm{d}\;{\alpha }_{2}\to 0$.

Considering the uncertainties in both E and J, we utilize the BFGS algorithm (Broyden, 1970; Fletcher, 1987) in order to minimize ${\chi }_{{\nu }}^{2}\left(\boldsymbol {a}\right)$ and achieve the best fit. The EPS formula fits well to the energy spectra of various energetic particle phenomena (Table 1 and Figure 4), including SEPs (electrons, protons, ³He and heavier ions), solar wind suprathermal particles (halo and superhalo electrons; pick-up ions and the suprathermal tail), ACRs, etc. The EPS fitting method can self-consistently determine the spectral features of energetic particle phenomena without first assuming a given spectral shape. Therefore, it can help us comprehensively understand the origin, acceleration, and transport processes of energetic particles over a wide energy spectrum, as well as comparatively examine the physical properties of different energetic particle phenomena.