Maximum Energy of Particles in Plasmas

Mitsuo Oka; Kazuo Makishima; Toshio Terasawa

doi:10.3847/1538-4357/ad9916

1. Introduction

Particles—both ions and electrons—are accelerated to very high nonthermal energies in space, solar, and astrophysical plasma environments (e.g., R. Blandford & D. Eichler 1987; T. Terasawa 2001; S. Krucker et al. 2008; M. Oka et al. 2018, 2023; F. Guo et al. 2024, and references therein). While the mechanism of particle acceleration has been a topic of major debate, a significant question remains regarding the maximum particle energy that can be reached in each plasma environment.

In an effort to constrain potential origins of ultra-high-energy (>10¹⁸ eV) cosmic rays (UHECRs), A. M. Hillas (1984) argued that the size L₀ of the essential part of the acceleration region must be greater than 2r_g, where r_g is the gyroradius of UHECRs. If particles are accelerated stochastically, the diffusive motion needs a very large space, much larger than the gyroradius, so that L₀ > 2r_g/β, where β = V/c, and V is the characteristic velocity of the scattering centers. Based on these considerations, A. M. Hillas (1984) presented what is now known as the Hillas diagram, which illustrates possible candidates of UHECR sources in the B − L₀ parameter space, where B is the magnetic field strength. It was proposed as a way to identify which sources fail to satisfy the requirements for UHECR sources.

Evidently, the arguments can be adapted in more generic discussions of maximum energy of particles in various plasma environments, well beyond the topic of search for UHECR production sites. The gyroradius of relativistic particles is expressed as r_g ∼ ε/cqB (SI unit) where ε is the particle energy, c is the light speed, and q is the particle charge. Then, the condition of L₀ > 2r_g gives the hard limit for the maximum particle energy achievable in a plasma environment:

$\begin{eqnarray}&&{\varepsilon }_{\,\rm{limit}\,}=qcBL\end{eqnarray} \tag{ 1 }$

where we omitted the factor of 2 by introducing the half-width of the acceleration region L = L₀/2. A caveat here is that the statement L > r_g is a necessary condition, rather than a sufficient condition, for particles to be trapped in the acceleration region. This is because particles might escape from the system even when L > r_g.

Hillas's second condition of L₀ > 2r_g/β can be rewritten to provide a more stringent limit of the maximum energy of particles

$\begin{eqnarray}&&{\varepsilon }_{{\rm{H}}}=qVBL,\end{eqnarray} \tag{ 2 }$

which we hereafter refer to as the Hillas limit. While A. M. Hillas (1984) envisioned the stochastic acceleration scenario, the same expression can be obtained by assuming a one-shot, direct acceleration by the system-wide motional electric field VB. In fact, similar forms of Equation (2) have been widely used in the context of magnetic reconnection where particles can be accelerated by the coherent reconnection electric field (e.g., W. H. Matthaeus et al. 1984; M. L. Goldstein et al. 1986; P. C. H. Martens 1988; B. Kliem 1994). K. Makishima (1999) also used such an idea of direct acceleration to arrive independently at the same scaling as Equation (2), and argued that the predicted ε_H is often close to the highest particle energies ε_obs found in various plasma environments where not only magnetic reconnection but also shocks could play a major role.

Interestingly, T. Terasawa (2001) obtained essentially the same limit as Equation (2) by assuming more explicitly that particles are accelerated stochastically in a diffusive environment. Using the spatial diffusion coefficient D, the size of the acceleration region can be roughly expressed as λ = D/V (e.g., R. Blandford & D. Eichler 1987). If we normalize D by the Bohm diffusion coefficient D_B = (1/3)r_gv (where v denotes the particle speed) and introduce η ≡ D/D_B, we obtain λ = ηD_B/V = ηε/3qVB where we used the relativistic gyroradius. Then, from the condition λ < L, we obtain the maximum attainable energy as

$\begin{eqnarray}&&{\varepsilon }_{\,\rm{diffusive}\,}=\frac{3}{\eta }qVBL.\end{eqnarray} \tag{ 3 }$

In the strong scattering limit (η = 1), this maximum energy is a factor of 3 higher than ε_H, although η is likely much larger than unity in solar and space plasmas. In the nonrelativistic regime, we can readily obtain ε_diffusive = (3/2η)qVBL, which is still comparable to ε_H in the strong scattering limit. Therefore, we should keep in mind that Equation (2) implies not only one-shot direct acceleration scenario but also the stochastic acceleration scenario.

A crucial point to note is that Equation (2) remains a necessary, but not a sufficient, condition (e.g., K. V. Ptitsyna & S. V. Troitsky 2010; K. Kotera & A. V. Olinto 2011; R. Alves Batista et al. 2019; J. H. Matthews et al. 2020). More stringent conditions may arise from additional processes. In fact, A. M. Hillas (1984) already discussed, for example, possible compositions of UHECRs, associated transport processes, and energy losses due to synchrotron radiation. More up-to-date descriptions of the limitation factors for the maximum particle energy are well summarized by J. H. Matthews et al. (2020), particularly their Figure 5. In other words, with the understanding that the Hillas limit is a ballpark estimate, the cosmic-ray physics community has sought more restrictive conditions in each plasma environment.

In contrast, K. Makishima (1999) and T. Terasawa (2001, 2011) took a different approach, aiming to validate the Hillas limit itself by comparing it with observations of high-energy particles and photons. They suggested that particle energies do reach the Hillas limit in various acceleration regions, even in solar and space plasma environments. On the other hand, a recent study suggests that the Hillas limit may overestimate observed maximum energies in many plasma environments (A. Chien et al. 2023). Regardless, these studies did not fully account for additional loss processes considered in the context of UHECR studies.

This Paper presents a novel study that integrates both approaches. Leveraging more recent observations of space, solar, astrophysical, and laboratory plasmas, we first attempt to validate the Hillas limit itself. However, when interpreting the results, we consider radiative energy losses as well as the scattering intensity, parameterized by η. A large value of η represents weak scattering, less effective confinement of particles in the acceleration region, and therefore a condition more restrictive than the Hillas limit. We also treat ions and electrons separately. This is because, although the Hillas limit does not depend on the particle mass, ions and electrons could experience different acceleration and energy-loss processes. Additionally, we try to use consistent definitions for V, B, and L for all plasma environments. With these new approaches, we have achieved a much better understanding of the maximum energy of particles in plasmas.

We first describe the methodology including the definitions of the V, B, L parameters (Section 2), followed by detailed descriptions of how we collected and compiled the parameter values, as well as the observed maximum energies of particles (Section 3). Then, we compare the observations and predictions in graphs (Section 4), discuss the results (Section 5), and conclude at the end (Section 6).

2. Methodology

We conducted a systematic review of the literature to deduce the key parameters (V, B, and L), predict the maximum particle energies ε_H via Equation (2), and compare these predictions with the highest, observed energies of protons (ε_obs,p) and electrons (ε_obs,e). We selected plasma environments from which we can confidently derive these key parameters as well as ε_obs. While heavy ions are also observed in some cases, we focus on protons in this study for simplicity.

There are three caveats when conducting this study. The first is that, in many cases, the observed maximum energies (ε_obs) represent only a lower limit of the possible maximum particle energy. This is because the higher-energy cutoff in power-law spectra is often unclear, particularly in solar and space plasmas. Power laws with exponential cutoffs (D. C. Ellison & R. Ramaty 1985; D. Band et al. 1993; Z. Liu et al. 2020) can be found in, for example, solar energetic particle events (e.g., D. C. Ellison & R. Ramaty 1985; R. A. Mewaldt et al. 2012) and crossings of Earth's bow shock (e.g., M. Oka et al. 2006; T. Amano et al. 2020). However, observed energy spectra do not always exhibit such a cutoff feature in many plasma environments. The ε_obs values are typically derived from the highest-energy channel of instruments that registered significant particle or photon counts. These values may be updated in the future with improved measurements that feature a larger energy range and higher sensitivity (or a larger dynamic range). Additionally, in the case of astrophysical objects, radiative energy losses could reduce the actual maximum energy of particles (K. V. Ptitsyna & S. V. Troitsky 2010). Therefore, ε_obs values that are lower than the Hillas limit do not automatically invalidate Equation (2). Similarly, ε_obs values that are comparable to the Hillas limit ε_H do not necessarily confirm its validity, as the actual maximum energy could be exceeding the Hillas limit. However, ε_obs ∼ ε_H provides strong evidence for particle acceleration mechanisms capable of reaching the Hillas limit energies.

Second, while direct in-situ measurements are available for space plasmas, remote-sensing measurements provide the key for diagnosing astrophysical plasmas. The observed maximum energy of photons can be identified, but it needs to be converted to the corresponding particle energy after knowing the emission mechanism and relevant parameters. For example, synchrotron radiation depends on the magnetic field B, and so, we need to treat the estimation of both B and ε_obs as a set of simultaneous equations.

Finally, V, B, and L should be defined in a way that is as consistent as possible across different plasma environments. As described in previous section, the Hillas limit can imply both one-shot direct acceleration and stochastic acceleration. Also, it does not necessarily distinguish shocks and magnetic reconnection. Therefore, we prefer to have a definition that does not explicitly depend on a specific acceleration mechanism. Nevertheless, the plasma condition and geometry differ greatly from acceleration site to site. To define V, B, and L, it is therefore still convenient to categorize plasma environments based on the type of primary source of energy. Figure 1 illustrates the categorization and our definitions of V, B, and L.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Schematic illustration of different types of plasma environments. Annotated on the left of each illustration are (1) the primary source of energy, (2) our definitions of V, B, and L, and (3) example plasma environments. This categorization is not meant to restrict the acceleration mechanism. For example, in the Type 2 environment, particles may be accelerated not only by magnetic reconnection but also by the termination shock formed by the dissipation of reconnection jets at the apex of a flaring loop. For more details, see Section 2 and Table 1.
Download figure:
Standard image High-resolution image

The first type of plasma environment (Figure 1, top) has the upstream bulk flow as the primary source of energy. Examples include planetary bow shocks and the termination shocks of the heliospheric solar wind, pulsar winds, and jets in radio galaxies. This configuration does not necessarily mean particles are accelerated by shocks only. It is certainly possible that particles are accelerated by magnetic reconnection that occurs, for example, in the current sheets embedded in solar and pulsar winds (e.g., J. F. Drake et al. 2010; G. P. Zank et al. 2014, 2015; Y. Lu et al. 2021). We define V as the bulk flow speed on the upstream side and B as the strength of the magnetic field carried by the flow. In astrophysical nonthermal sources, a magnetic field strength may be obtained through observations of emissions from the postshock region, but we regard such a value as the magnetic field on the downstream side. Upon application of the Hillas limit, we reduce such a postshock value by dividing by the shock compression ratio to derive the upstream value. The spatial scale L is approximated by the curvature radius of the shock front. In the case of shell-type supernova remnants (SNRs), we use the radius of the shell instead of the thickness of the shell, although we also provide thickness-based estimates for reference. In the case of a radio galaxy, twin jets from its nucleus terminate via shocks, to form gigantic radio lobes filled with energetic particles. Although the shock front itself may not be spatially resolved via imaging in any wavelengths, we can still use the radius of the so-called hot-spots that sit at lobe centers, which is clearly visible in radio images. In terrestrial/planetary bow shocks, the curvature radius could roughly correspond to the half-width of the obstruction, i.e., the magnetotail. It should also be roughly comparable to the standoff distance between the nose of the shock and the central planet.

The second type of plasma setting (Figure 1, center) has accumulated magnetic flux in the upstream region as the primary source of energy. In such a situation, the inflowing antiparallel fields reconnect to produce bidirectional jets, one of which impinges on quasi-stationary magnetic fields anchored to some thick materials. Examples include solar flares and Earth's magnetosphere. This categorization does not necessarily mean we exclude the possibility of other particle acceleration mechanisms. In solar flares, it is certainly possible that particles are accelerated by the termination shocks located at or above the top of the magnetic loop on the downstream side (e.g., S. Masuda et al. 1994; S. Tsuneta & T. Naito 1998; B. Chen et al. 2015). We define B as the characteristic (asymptotic) magnetic field on the upstream side of reconnection, such as the lobe magnetic field in Earth's magnetotail and ambient magnetic field in the solar corona. We define V as the typical speed of the fast flows in the system to obtain an estimate of the maximum attainable energy in the system. It should be comparable to the reconnection outflow speed or the Alfvén speed V_A. Planetary radiation belts also fall into this category. Although the solar wind dynamic pressure could be contributing to a global-scale compression of the inner magnetosphere and hence particle energization in the radiation belts, observations also indicate that a substantial amount of particle and energy flux is "injected" into the inner magnetosphere from the magnetotail.

The third type of plasma configuration involves a rapidly rotating, strongly magnetized planet or star, with the electric potential induced by the corotating magnetic field. Examples include the magnetospheres of Jupiter and pulsars such as the Crab Pulsar. While the precise mechanisms of particle acceleration remain unclear in such environments, a thin layer of charge separation that we refer to as a double-layer is often considered. As in the former two categories, the acceleration could also be driven by other mechanisms, such as magnetic reconnection and pitch angle scattering by waves. We define V as the maximum corotation speed, L as the radial distance from the central object to the location of the maximum corotation speed, and B is the typical magnetic field near the radial distance of L.

3. Data Collection

Based on the methodology described above, we searched the literature and collected data, as summarized in Table 1. In this section, we describe the details of the collected data.

Table 1. The Hillas Limit and the Observed Maximum Energy of Particles in Various Plasma Environments

Plasma Environment		Flow Speed	Magnetic Field	System Size	Hillas Limit	Observations
Type	Name	V	B	L	ε_H	ε_obs,p	ε_obs,e
		(km s⁻¹)	(nT)	(m)
2	Earth's magnetotail	300–1000	15–25	(9.6−16) × 10⁷	0.4–4 MeV	1 MeV		2 MeV
2	Earth's radiation belt	300–1000	15–25	(9.6−16) × 10⁷	0.4–4 MeV	800 MeV		15 MeV
1	Earth's bow shock	400–600	5–8	(9.6−16) × 10⁷	190–760 keV	300 keV		300 keV
2	Mercury's magnetosphere	300–400	40–50	(4.9–9.8) × 10⁶	59–200 keV	550 keV		300 keV
3	Jupiter's radiation belt	400–800	5–15	(5.2–8.2) × 10⁹	10–98 MeV	80 MeV		100 MeV
3	Saturn's radiation belt	200–500	5–10	(1.2–2.9) × 10⁹	1–15 MeV	300 MeV		2 MeV
1	Saturn's bow shock	300–600	0.2–1.0	(2.9–3.5) × 10⁹	0.2–2 MeV	⋯		700 keV
1	Heliosphere (ACR)	250–350	0.02–0.1	(1.2–1.8) × 10¹³	60–630 MeV	100 MeV		⋯

2	Solar flares	1000–3000	(5–50) × 10⁶	(1.0–3.0) × 10⁷	0.05–5 TeV	30 GeV		45 MeV
1	CME (shock)	1000–4000	300–10,000	(3.5–10) × 10⁹	1–420 GeV	30 GeV		⋯

1	SN 1006	2700–4600	0.1–0.5	2.7 × 10¹⁷	73–620 TeV	100 TeV		10 TeV
1	RX J1713.7-3946	3900	0.1–1.0	1.4 × 10¹⁷	53–530 TeV	100 TeV		100 TeV
1	Crab Nebula	3 × 10⁵	1–4	4.3 × 10¹⁵	2–6 PeV	⋯		2 PeV
3	Crab Pulsar	3 × 10⁵	1.0 × 10¹¹	1.5 × 10⁶	45 PeV	⋯		15 TeV
1	Cygnus A	72,000–90,000	9–10	(3.4–4.6) × 10¹⁹	22–42 EeV	⋯		30 GeV

1	Laser plasma (shock p)	1500	2.0 × 10¹⁰	(3.0–4.0) × 10⁻³	90–120 keV	80 keV		⋯
1	Laser plasma (shock e-)	1000	1.0 × 10¹¹	5.0 × 10⁻³	500 keV	⋯		500 keV
2	Laser plasma (mrx e-)	500–1100	5.0 × 10¹⁰	1.0 × 10⁻³	25–55 keV	⋯		40–70 keV

Note. The first column shows the type of each environment (either 1, 2, or 3). See Section 2 for more details. The plasma environments are listed in the order of appearance in Section 3. The average radii of the planets are as follows: Earth R_E ~ 6.371 × 10⁶ m, Mercury R_M ~ 2.440 × 10⁶ m, Jupiter R_J ~ 6.9911 × 10⁷ m, and Saturn R_S ~ 5.8232 × 10⁷ m. The average solar radius is ~6.957 × 10⁸ m. For unit conversions, 1 nT = 10 μG, where nT is nano Tesla and μG is micro gauss. For the laser plasma environments, experiments are performed for collisionless shock (shock) and magnetic reconnection (mrx) processes for both protons (p) and electrons (e-).

Download table as: ASCII Typeset image

3.1. Space Plasmas

Let us start from Earth's magnetotail because it has been extensively investigated over the past decades. To estimate ε_H, we consider "near-Earth" reconnection that takes place in the down-tail distances of ∼15–25R_E where R_E is Earth's radius (e.g., Y. Miyashita et al. 2009). In this region, the lobe magnetic field is B = 15–25 nT (D. H. Fairfield & J. Jones 1996), and the half-width of the magnetotail is L = 15–25R_E or (9.6–16) × 10⁷ m (D. H. Fairfield 1971). For the speed V, we use the typical speed of bursty bulk flows V = 300–1000 km s⁻¹ in the central plasma sheet (e.g., V. Angelopoulos et al. 1988; L. Q. Zhang et al. 2016). While the flow channel appears to be much narrower than the full width of the magnetotail, on the order of a few R_E (e.g., V. Angelopoulos et al. 1996), we adhere to our definition described in Section 2 and use the magnetotail half-width for L. Consequently, we obtain ε_H = 0.4–4.0 MeV, while earlier observations have reported up to ∼1 MeV for both ions and electrons (e.g., T. Terasawa & A. Nishida 1976; S. M. Krimigis 1979; S. P. Christon et al. 1988; S. P. Christon et al. 1989; S. P. Christon et al. 1991). A recent study also reported a detection of ∼2 MeV electrons that were precipitating into Earth's atmosphere from the magnetotail (A. V. Artemyev et al. 2024). It is worth noting that earlier studies already argued that the direct acceleration by the reconnection electric field does not explain ∼MeV particles (e.g., T. Terasawa & A. Nishida 1976; S. M. Krimigis 1979). In fact, the polar cap potential, which is formed when the dayside magnetopause reconnection enables the solar wind to drive ionospheric convection, is known to saturate only at ∼200 kV (M. R. Hairston et al. 2005; S. G. Shepherd 2007).

Earth's radiation belt is also a distinct region of particle acceleration, but the values of B vary largely and the values of V can be fluctuating around ∼0 due to its oscillatory motion, making it difficult to choose a characteristic value of B and V. Therefore, we use the characteristic values in the magnetotail, as also described in Section 2. This is because, while electrons can be accelerated locally inside the radiation belt, a substantial fraction of their energy originates from the magnetotail in the forms of waves and moderately accelerated seed particles (e.g., A. N. Jaynes et al. 2015; K. A. Sorathia et al. 2018; D. L. Turner et al. 2021, and references therein). The highest energies of observed particles are ∼800 MeV for protons (J. E. Mazur et al. 2023) and ∼15 MeV (or possibly higher) for electrons (J. B. Blake et al. 1992; A. L. Vampola & A. Korth 1992). However, such high-energy protons are likely to be supplied by galactic cosmic rays (GCRs). It has been known that GCR particles knock neutrons out from the atmospheric atoms, and such neutrons decay to produce protons (e.g., S. F. Singer 1958a, 1958b; W. N. Hess 1959; Y. Li et al. 2023). This is called cosmic-ray albedo neutron decay (CRAND). More details of CRAND are discussed in Section 5.4.

For Earth's bow shock, we consider the typical values of the "fast" solar wind at 1 au (e.g., C. Larrodera & C. Cid 2020), i.e., B = 5–8 nT and V = 400–600 km s⁻¹. The system size L should be roughly the same as that of the magnetotail, L = 15–25R_E, or (9.6–16) × 10⁷ m. As a result, we expect ε_H ∼ 190–760 keV, while observations in the shock upstream region find energetic protons of up to ∼300 keV (e.g., M. Scholer et al. 1981; D. L. Turner et al. 2018; K. J. Trattner et al. 2023) and electrons of similar energies (e.g., L. B. Wilson et al. 2016), although these energetic protons may be partially contributed by those from the magnetosphere (M. Scholer et al. 1981). A recent observation of hot flow anomalies—transient concentrations of shock-reflected, suprathermal ions that occur when certain structures in the solar wind collide with the bow shock—reported energizations of ions (D. L. Turner et al. 2018). In their report, the energy spectrum of H⁺, He²⁺, and O⁶⁺ extended up to ∼100, ∼300, and ∼600 keV, while showing a peak at ∼30, ∼120, and ∼300 keV, respectively. While the authors successfully interpreted the peak energies by a Fermi-type acceleration "trap" between two converging magnetic mirrors, such a dependence on charge state could also be attributed to Equation (2).

In Mercury's magnetosphere with its half-width being L ∼ 2–4R_M where R_M is Mercury's radius or (4.9–9.8) × 10⁶ m (R. M. Winslow et al. 2013), protons of up to ∼550 keV (J. A. Simpson et al. 1974) and electrons of up to ∼300 keV (J. A. Simpson et al. 1974; G. C. Ho et al. 2011, 2012; D. J. Lawrence et al. 2015) have been detected. The magnetotail reconnection appears to occur around ∼3R_M down-tail, and the lobe magnetic field in such a region is ∼40–50 nT (G. Poh et al. 2017; C. F. Bowers et al. 2024). The bulk flow speed V has been neither directly nor routinely measured, but we assume 300–400 km s⁻¹ based on flow speeds associated with dipolarization (R. M. Dewey et al. 2018) and flux ropes (FTEs; J. A. Slavin et al. 2012). This is comparable to the estimated thermal speed of ions (J. M. Raines et al. 2011). As a result, we obtain ε_H = 59–200 keV.

Jupiter's plasma environment is very different from those of rocky planets such as Earth and Mercury. Its dynamics is largely driven by the rapid planetary rotation, assisted by the significant internal mass loading from the volcanic moon Io (e.g., C. M. Jackman et al. 2014; R. Guo & Z. Yao 2024, and references therein). A key picture here is that, at a certain distance from the central planet, the magnetic field is stretched outward, and the corotation breaks down due to magnetic reconnection (V. M. Vasyliunas 1983). This is consistent with Type 3 illustrated in Figure 1, and we use the typical radial distance at which the magnetic reconnection takes place, i.e., L ∼ 90–140R_J where R_J is Jupiter's radius, or (5.2–8.2) × 10⁹ m (M. F. Vogt et al. 2010, 2020). The lobe magnetic field at such radial distances is 5–15 nT (M. G. Kivelson & K. K. Khurana 2002), and the typical speed of fast plasma flows is V ∼ 400–800 km s⁻¹ (E. A. Kronberg et al. 2008). As a result, we obtain ε_H = 10–98 MeV. For observation, synchrotron emissions from ≳50 MeV electrons in the radiation belt was first reported by S. J. Bolton et al. (2002), and it was suggested that there is a maximum cutoff at ∼100 MeV (I. de Pater & D. E. Dunn 2003). We also add that nonthermal bremsstrahlung X-rays of up to 20 keV have been detected from the polar (auroral) regions of Jupiter (K. Mori et al. 2022). For protons, the maximum observed energy of ∼80 MeV was recorded near the main ring of Jupiter (H. M. Fischer et al. 1996; P. Kollmann et al. 2017). Notably, heavy ions of ≳125 MeV are also detected based on ionization signatures in the star camera on board the Juno spacecraft (H. N. Becker et al. 2021).

Saturn's plasma environment is in many ways intermediate between those of Earth and Jupiter (C. M. Jackman et al. 2014). However, the effect of planetary rotation appears to be dominant (e.g., R. Guo & Z. Yao 2024), and so, we consider Saturn's magnetosphere as one of the third type described in Section 2. On average, the X-line is located at distances of 20–30R_S where R_S is Saturn's radius, but it is highly variable and could be located as far out as ∼50R_S (H. J. McAndrews et al. 2009; A. W. Smith et al. 2016). Thus, we adopt the spatial scale of L = 20–50R_S, or (1.2–2.9) × 10⁹ m. The lobe magnetic field in such locations is 5–10 nT (C. M. Jackman et al. 2014). With the possible flow speeds V = 200–500 km s⁻¹ (e.g., T. W. Hill et al. 2008; H. J. McAndrews et al. 2009), we obtain ε_H = 1–15 MeV. In reality, protons with energies up to at least ∼300 MeV (e.g., N. Krupp et al. 2018; E. Roussos et al. 2018) and electrons up to at least ∼2 MeV (e.g., C. Paranicas et al. 2010; E. Roussos et al. 2014) have been observed in Saturn's radiation belt. Here, the high-energy protons are likely CRAND protons that also exist in Earth's radiation belt.

At Saturn's bow shock, electrons up to ∼700 keV were detected when the solar wind Alfvén Mach number M_A was unusually high ∼100 (A. Masters et al. 2013). Based on the solar wind parameters at Saturn, i.e., B = 0.2–1.0 nT and V = 300–600 km s⁻¹ (e.g., A. Masters et al. 2013; M. F. Thomsen et al. 2019; D. J. Gershman et al. 2024), as well as the half-width of the magnetosphere L = 50–60R_S, or (2.9–3.5) × 10⁹ m (e.g., H. J. McAndrews et al. 2009; C. M. Jackman et al. 2019), we estimate ε_p = 0.2–2 MeV.

For the heliosphere, we consider anomalous cosmic rays (ACRs) that appear as an excess flux in the lowest energy range of the galactic cosmic-ray spectra (e.g., J. Giacalone et al. 2022, and references therein). In today's standard scenario, neutral particles of interstellar origin are ionized deep inside the heliosphere and then "picked-up" by the solar wind to be transported back toward the outer edge of the heliosphere. They are likely accelerated by the heliospheric termination shock to become ACRs, although the precise acceleration mechanism has been a topic of debate, especially after recent observations by Voyagers 1 and 2. The energies of ACR protons can reach ∼100 MeV (e.g., A. C. Cummings & E. C. Stone 1998). J. R. Jokipii & J. Giacalone (1998) already argued that these highest energies are gained from electrostatic potential while drifting along the shock surface, which they estimated as 200–300 MeV. In our case, we use B = 0.02–0.1 nT and V = 250–350 km s⁻¹ as they were measured just upstream of the termination shock by the Voyager spacecraft (L. F. Burlaga et al. 2008). For L, we use the expected average distance between the Sun and the termination shock along the "nose" direction, i.e., 80–120 au, or (1.2–1.8) × 10¹³ m. As a result, we obtain ε_H = 60–630 MeV.

3.2. Solar Plasmas

In solar flares, the Bremsstrahlung emission from energetic electrons produces the X-ray/γ-ray continuum up to a few tens of MeV (e.g., R. P. Lin et al. 2002). In such remote-sensing measurements, the precise evaluation of the maximum energy of electrons is confounded by the photons produced by the decay of pions (e.g., N. Vilmer et al. 2003). However, in-situ measurements reported ∼45 MeV electrons from solar flares (P. Evenson et al. 1984; D. Moses et al. 1989).

To estimate ε_H for solar flares, we consider spatially large flares that fit into the standard model where magnetic reconnection in the corona plays a major role (e.g., S. Masuda et al. 1994; K. Shibata et al. 1995). This is a Type 2 plasma environment. We therefore adopt the values of V and B in the corona, even though a large part of flare emission originates from the chromospheric footpoints of flaring loop. For the flow speed V, we assume the typical Alfvén speed, 1000–3000 km s⁻¹ (e.g., M. J. Aschwanden 2005). For the magnetic field, we assume (5.0–50) × 10⁶ nT, or B = 50–500 G (e.g., M. J. Aschwanden 2005; B. Chen et al. 2020). For the system size L, we assume L = 10–30 Mm, or (1.0–3.0) × 10⁷ m, which is half the loop-size of spatially large flares (e.g., M. J. Aschwanden et al. 1996; M. J. Aschwanden 2005; S. Krucker et al. 2008; M. Oka et al. 2015). Alternatively, we could use the length of flare arcade (or flare ribbon) that extends in the direction of reconnection electric field. It can be ∼200 Mm or even larger. However, we keep the definition outlined in Section 2 and use the flare loop size. In fact, hard X-ray sources are usually isolated or fragmented and do not extend along the flare arcade (e.g., A. Asai et al. 2002; G. Shi et al. 2024), suggesting that the acceleration region does not extend as long as the length of the arcade. With these parameters, we obtain ε_H = 0.05–5 TeV.

For high-energy protons of solar origin, we consider observations by ground-based, cosmic-ray monitors. During solar eruptive events including solar flares and coronal mass ejections (CMEs), significant enhancements of GeV proton flux are observed, which are referred to as ground level enhancements (GLEs; e.g., H. Carmichael 1962). To date, protons of up to ∼30 GV in rigidity (approximately 30 GeV in energy) have been detected (J. L. Lovell et al. 1998). One intriguing aspect of these observations is that the spectra often exhibit a roll-off, indicating that the power-law spectra are unlikely to extend further, even with improved measurements of GLEs. Therefore, we do not consider the observed value of 30 GeV as a lower limit.

Unfortunately, the origin of GLE protons remains ambiguous. They can be produced by solar flares (in which magnetic reconnection plays an important role), CME-associated shocks, or both (S. Kahler 1994; E. V. Vashenyuk et al. 2011; M. J. Aschwanden 2012; J. Guo et al. 2023; K. G. McCracken et al. 2023). Some observations and modeling indicate that GLEs often exhibit a two-step evolution: a prompt enhancement possibly due to a solar flare, followed by a gradual decay possibly due to a CME shock (e.g., E. V. Vashenyuk et al. 2011). However, because of the lack of consensus on the origin, we assume that both solar flares and CMEs produce protons of ∼30 GeV. It should be noted that recent observations by Fermi's Large Area Telescope reported γ-ray spectra that are consistent with the decay of pions produced by >300 MeV protons (M. Ajello et al. 2021). They found at least two distinct types of γ-ray emission: prompt-impulsive and delayed-gradual, which could be related to the two-step evolution of GLE events.

To estimate ε_H for CMEs (and associated shocks), we take into account the statistical study by S. Kahler (1994). They reported that GLEs exhibited a peak when the height of associated CME reached 5–15R_⊙ (S. Kahler 1994). Then, we use L = 5–15R_⊙, or (3.5–10) × 10⁹ m, combined with B = 300–10,000 nT, or 0.003–0.1 G (M. K. Bird & P. Edenhofer 1990), and V = 1000–4000 km s⁻¹ (D. F. Smart & M. A. Shea 1985), to obtain ε_H = 1–420 GeV for CMEs in Table 1.

3.3. Astrophysical Plasmas

3.3.1. SN 1006

SNRs are one of the most important astrophysical plasma environments in the context of the Hillas limit. This is because SNRs are promising sites for cosmic-ray production up to the "knee" energy (∼TeV) or even beyond, and whether they can actually achieve this energy has been a subject of intense debate, both observationally and theoretically (e.g., P. O. Lagage & C. J. Cesarsky 1983; A. R. Bell 2004; A. M. Hillas 2005; L. O. Drury 2012; H. Suzuki et al. 2022; R. Diesing 2023; M. Tao et al. 2024).

The Hillas limit ε_H for SN 1006 is obtained as follows, using our definitions for Type 1 environment. The remnant of SN 1006 exhibits a shell structure with the average radius of $\sim 15^{\prime}$ (e.g., F. F. Gardner & D. K. Milne 1965; K. Koyama et al. 1995; A. Bamba et al. 2003). Using the distance of ∼2 kpc (e.g., P. F. Winkler et al. 2003; D. A. Green 2019) where 1 pc = 3.1 × 10¹⁶ m, we have L ∼ 2.7 × 10¹⁷m. The remnant is interpreted as a shock front expanding at the speed of 0 $\mathop{.}\limits^{\unicode{x02033}}$ 28–0 $\mathop{.}\limits^{\unicode{x02033}}$ 48 yr^–1 (P. F. Winkler et al. 2003; S. Katsuda et al. 2009) or 2700–4600 km s⁻¹. For the upstream magnetic field, we use the typical interstellar magnetic field of 0.1–0.5 nT, or 1–5 μG. As a result, we obtain ε_H ∼ 73–620 TeV.

It should be noted that a higher value of ε_H can be obtained by focusing on the local conditions at the shock front. In the diffusive shock acceleration scenario (e.g., R. Blandford & D. Eichler 1987), particles are accelerated within a finite region around the shock front. Consequently, the thickness of the shell measured in X-rays and γ-rays, approximately 1 $^{\prime}$ (A. Bamba et al. 2003; J.-T. Li et al. 2018), should correspond to the spatial extent of the acceleration region (e.g., R. Yamazaki et al. 2004). For the magnetic field B, recent observations have found evidence for a highly amplified value at the shell, B ≳ 20 nT (M. Tao et al. 2024), consistent with earlier expectations (e.g., E. G. Berezhko et al. 2002; E. Parizot et al. 2006). Therefore, by adopting B ∼ 20 T and $L\sim 1^{\prime} \rm{}\,\rm{}\,$ , or 1.8 × 10¹⁶ m, we obtain a predicted maximum energy ε_H ∼ 2 PeV. This is an order of magnitude higher than our first estimate.

Let us turn to the observational estimates ε_obs, using X-ray and γ-ray remote sensing. The observed X-ray spectra exhibit a clear cutoff at a photon energy of ∼0.2 keV, which can be interpreted as synchrotron emission from electrons with the maximum energy of 10 TeV in a magnetic field of 1 nT (A. Bamba et al. 2008). On the other hand, observed γ-rays can be attributed to inverse Compton scattering of cosmic microwave background (CMB) photons by electrons, or proton-induced pion-decay, or a combination of both. Based on a hybrid model that takes into account both mechanisms of γ-ray production, it has been reported that the maximum energies of electrons and protons are ε_obs,e ∼ 10 TeV and ε_obs,p ∼ 100 TeV, respectively, with a magnetic field strength of 4.5 nT (F. Acero et al. 2010). These energies match with our estimates within an order of magnitude.

3.3.2. RX J1713.7-3946

RX J1713.7-3946 is another shell-forming SNR, roughly 30 $^{\prime} \rm{}\,\rm{}\,$ in diameter at the distance of ∼1 kpc (e.g., N. Tsuji & Y. Uchiyama 2016), giving us L ∼ 1.4 × 10¹⁷m. The expansion speed is 3900 km s⁻¹ (N. Tsuji & Y. Uchiyama 2016), which we use for V. A key difference from SN 1006 is that it is located right on the Galactic plane where the interstellar medium is much denser and possibly nonuniform. In fact, a part of the shell is moving substantially slower, indicating an ongoing interaction with local clouds (T. Tanaka et al. 2020). Furthermore, the presence of dense interstellar material appears to enable the so-called hadronic channel, in which accelerated protons collide with ambient protons to produce pions. These pions subsequently decay into a pair of 70 MeV (in rest-frame) γ-ray photons, thereby contributing significantly to the observed flux of γ-ray emission. For the upstream magnetic field, we assume 0.1–1 nT, which includes values slightly higher than those we used for SN 1006, considering the presence of denser materials. Then, we predict the maximum energy of particles as ε_H ∼ 53–530 TeV, very similar to that of SN 1006. Alternatively, we can again consider the local condition with amplified magnetic field of B ∼ 10 nT (Y. Uchiyama et al. 2007) and shell thickness $L\sim 0.5^{\prime} \rm{}\,\rm{}\,$ (T. Tanaka et al. 2020) to obtain ε_H ∼ 170 TeV, a value consistent with the above estimation.

For observations, the spectra have been measured with both X-rays (e.g., T. Takahashi et al. 2008; Q. Yuan et al. 2011; N. Tsuji et al. 2019) and γ-rays (e.g., A. A. Abdo et al. 2011; Q. Yuan et al. 2011; H. Abdalla et al. 2018). The observed maximum energy of particles is obtained as ε_obs,p ∼ ε_obs,e ∼ 100 TeV, although it can vary within an order of magnitude depending on the spectral model employed, including how much of a fraction of ions contribute to the photon emissions (hadronic versus leptonic). In this regard, a recent study that takes into account the column density of the local interstellar medium concludes that ions and electrons constitute about 70% and 30% of the total γ-ray emissions, respectively (Y. Fukui et al. 2021). This result provides valuable evidence for proton acceleration in astrophysical settings.

3.3.3. Crab Nebula

The Crab Nebula is also an SNR, but the presence of a spinning pulsar at its center makes it very different from more ordinary shell-type SNRs (e.g., M. S. Longair 2011; R. Buhler & R. Blandford 2014, and references therein). A distinctive feature in its X-ray image is the torus and the associated inner ring (M. C. Weisskopf et al. 2000), which is interpreted as the termination shock formed as the pulsar wind—a magnetized, relativistic flow of electrons and positrons from the pulsar—collides with the surrounding nebula (C. F. Kennel & F. V. Coroniti 1984). To estimate ε_H, we use the radius of the inner ring ∼0.14 pc, or 4.3 × 10¹⁵ m (M. C. Weisskopf et al. 2000) for L and the light speed c ∼ 3 × 10⁵ km s⁻¹ for V. The magnetic field B in the postshock region can be estimated from the spectral break that appears in the synchrotron spectrum at ν ∼ 10¹³ Hz. Assuming that the spectral break represents the condition of the synchrotron cooling time being comparable to the age of the nebula, we can obtain B ∼ 30 nT (P. L. Marsden et al. 1984). Similar estimates in the range of 10–30 nT have been obtained by other studies (e.g., R. Buhler & R. Blandford 2014, and references therein). To further estimate the magnetic field in the upstream region of the shock, we adopt the compression ratio 7 for a relativistic, strictly perpendicular shock and obtain B ∼ 1–4 nT. Combining the parameters, we obtain the Hillas limit of ε_H ∼ 2–6 PeV.

Regarding observations, the electromagnetic spectrum from Crab Nebula has been observed over a wide range of energies from radio to γ-ray ranges, revealing both components of synchrotron and inverse Compton emissions (e.g., X. Zhang et al. 2020; F. Aharonian et al. 2024). According to these observations, the synchrotron component exhibits a clear cutoff at the photon energy of ∼50 MeV. Using the magnetic field 30 nT obtained above and the formula for the synchrotron critical frequency, we find that this photon energy corresponds to an electron energy of ε_obs,e ∼ 2 PeV. The inverse Compton component extends beyond 0.1 PeV (F. Aharonian et al. 2024), which corresponds to the electron energy of ε_obs,e ∼ 0.5 PeV. However, the cutoff is not visible in the spectra, and so, we regard this value as a lower limit of ε_obs,e.

It should be noted that there may be a nonnegligible contribution to the spectrum from energetic protons, although the leptonic model has been widely accepted (e.g., A. M. Atoyan & F. A. Aharonian 1996). In this regard, a recent broadband spectral analysis considered hadronic photons from the neutral pion-decay process and showed that the fraction of energy converted into energetic protons is likely <0.5% although it could be as high as 7% if only the γ-ray data are used (X. Zhang et al. 2020). Therefore, we consider that the the plasma in the Crab Nebula consists of predominantly electrons and positron. Consequently, we did not include protons in Table 1.

3.3.4. Crab Pulsar

The Crab Pulsar is not spatially resolved by remote-sensing measurements, and much of our evaluation of the Hillas limit is based on the arguments and parameters described in classical literature (e.g., P. Goldreich & W. H. Julian 1969; M. S. Longair 2011; R. Buhler & R. Blandford 2014). The Crab Pulsar produces intense beams of radiation through acceleration of particles within its magnetosphere. The particle acceleration mechanism remains largely unknown, and there exist different models, such as polar cap, outer gap, and slot gap models, depending on the assumed location of the acceleration site. An important spatial scale in pulsar's magnetosphere is the radius of the light cylinder R_LC. This is where the speed of the magnetic field corotating with the pulsar reaches the light speed, so that R_LC = c/Ω ∼ 1.5 × 10⁶m where Ω = 200 rad s⁻¹ is the pulsar's angular velocity. Therefore, to estimate the predicted maximum energy of particles, we adopt the definitions for Type 3 environment and use L = R_LC and V = c ∼ 3 × 10⁵ km s⁻¹. For the magnetic field strength B, it scales like ∝r⁻³ where r is the distance from the pulsar and B ∼ 4 × 10¹⁷ nT, or 4 × 10¹² G, at the equatorial surface of the neutron star. Thus, we obtain 10¹¹ nT at the surface of the cylinder. Then, the predicted maximum energy becomes ∼45 PeV.

The Crab Pulsar produces pulsed emissions, and photon spectra that extend up to ∼TeV have been measured, despite the flux being 2 orders of magnitude smaller than that of the nebula (S. Ansoldi et al. 2016; P. K. H. Yeung 2020). Assuming inverse Compton scattering of CMB photons by electrons, this maximum photon energy corresponds to the maximum particle energy of ε_obs ∼ 15 TeV. We consider this as a lower limit of ε_obs.

3.3.5. Cygnus A

As touched on in Section 2, jet-terminal lobes of radio galaxies (particularly of the so-called Fanaroff–Riley Type II, or FR-II) are also an interesting astroplasma environment. It is generally considered that particle acceleration occurs predominantly at the termination (reverse) shock, creating hot-spots in a central region of their radio lobes (e.g., C. L. Carilli & P. D. Barthel 1996; F. Casse & A. Marcowith 2005). In this study, we specifically consider Cygnus A (3C405) because it is one of the brightest radio galaxies and has been studied intensively from various approaches. Based on observations by, for example, K. Meisenheimer et al. (1997) and K. C. Steenbrugge & K. M. Blundell (2008), combined with our guideline for the Hillas limit (Section 2), we use the estimated radius of the hot-spots, L = 1.1–1.5 kpc, or (3.4–4.6) × 10¹⁹ m, and the jet speed of V = (0.24–0.30)c, or 72,000–90,000 km s⁻¹. For the magnetic field B, we assume its value in the hot-spot 40 nT (K. Meisenheimer et al. 1997), and divide it by the shock compression ratio 3.9–4.3 to obtain the shock upstream value 9–10 nT. These numbers bring the predicted maximum energy ε_H very high, 22–42 exa-electronvolt (EeV) , up to the range of UHECRs, ≳EeV. In fact, radio galaxy lobes were some of the candidates for the site of UHECR production based on earlier evaluations of the Hillas limit (e.g., A. M. Hillas 1984; J. P. Rachen & P. L. Biermann 1993; C. A. Norman et al. 1995).

In contrast to the above predictions, the observed maximum energy of electrons is much lower, ε_obs ∼ 30 GeV in the hot-spots of Cygnus A (K. Meisenheimer et al. 1997). Similar values of ε_obs have been obtained for the lobe regions (K. C. Steenbrugge & K. M. Blundell 2008; Y. Yaji et al. 2010). This issue is not specific to Cygnus A, because FR-II radio galaxies generally show synchrotron emissions predominantly in the radio range and do not extend to the higher-energy X-ray range. The observations in the X-ray range are interpreted as the inverse Compton scattering of CMB photons by relativistic electrons (e.g., E. D. Feigelson et al. 1995; H. Kaneda et al. 1995). Therefore, the observed maximum energies are robust, and it is clear that electrons do not reach the Hillas limit in the hot-spots of FR-II radio galaxies. We discuss this gap in Section 5.

3.4. Laboratory Plasmas

Laboratory experiments provide unique opportunities to study plasma phenomena in a controlled manner. In particular, recent developments in laser facilities have enabled experimental studies of particle acceleration at shocks and reconnection in collisionless conditions.

W. Yao et al. (2021) observed 80 keV protons at a shock front generated by a laser-driven, fast flow of plasma expanding into an ambient cloud of plasma. Using an earlier phase of the evolution when the shock front was still moving at a high speed, they evaluated the Hillas limit with the shock speed V ∼ 1500 km, the shock size L ∼ (3–4) × 10⁻³ m, and the ambient magnetic field 2 × 10¹⁰ nT, to obtain ε_H ∼ 90–120 keV, which is comparable to the observed energy.

F. Fiuza et al. (2020) detected 500 keV electrons at a shock front generated by an interaction of laser-driven, counterstreaming flows of plasma. They reported the Alfvén Mach number M_A ∼ 400. They also evaluated the Hillas limit with the following parameters: the flow speed V ∼ 1000 km s⁻¹, the transverse size of the shock L ∼ 5 × 10⁻³ m, and the magnetic field ∼10¹¹ nT, or 1 MG, leading to the expected maximum energy of ε_H ∼ 500 keV, well consistent with the observation.

A. Chien et al. (2023) reported a detection of 40–70 keV electrons during magnetic reconnection in a laser-driven plasma experiment. The reconnection geometry was formed by currents induced in two U-shaped coils by the same laser-driven plasma. During the initial "push" phase, when the currents were still increasing, the reconnection was more strongly driven than compared to a later phase. Their simulations of this experiment indicated that the peak rate of reconnection was α_R ∼ 0.6 where α_R is defined as the inflow speed V_in normalized by the Alfvén speed V_A, i.e., α_R = V_in/V_A. Using data from this earlier phase, along with the inferred peak value of α_R ∼ 0.6, they calculated the reconnection electric field to be 0.6V_AB₀, where V_A ∼ 500–1100 km s⁻¹ and B₀ ∼ 5 × 10¹⁰ nT. With a system size of ∼10⁻³ m (or 1.4 times the ion inertia length), they derived a predicted maximum energy of 13–30 keV, which is roughly half the observed energy.

A caveat here is that this prediction was based on the inflow speed of 0.6V_A, whereas our Hillas limit evaluation is based on the outflow speed, or equivalently V_A (Section 2). Thus, we estimate ε_H without the factor 0.6 and obtain ε_H ∼ 25–55 keV, which is comparable to the observed values. It should also be noted that, according to A. Chien et al. (2023), the system size was so small that the observed reconnection was deeply in the electron-only regime, where ions are decoupled from the process and do not exhibit ion-scale jets (T. D. Phan et al. 2018). In the presented systematic review, this is the smallest physical system used for testing the Hillas limit.

4. Results

All the values compiled in Table 1 are visualized in Figures 2 and 3 for protons and electrons, respectively. In both panels, the horizontal axis shows the maximum energy ε_H predicted by Equation (2), derived from the parameters V, B, and L, whereas the vertical axis shows the observed maximum energy ε_obs. For those plasma environments in which ε_obs is the lower limit of the observations, we added an upward arrow. The black solid line is the identity line of ε_obs = ε_H. It is labeled "direct" because the Hillas limit can be interpreted by "direct" (one-shot) acceleration.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** The maximum energy of protons in various plasma environments, as compared between the Hillas limit prediction (horizontal axis) and observations (vertical axis). See Table 1 for the values used and Section 3 for their justification. The black line is the identity line and represents the scenario of direct, one-shot acceleration. The gray shaded region represents the scenario of stochastic and diffusive acceleration with a strong scattering, η = 1–100. The dashed curve represents the trapping condition r_g < L for the case of ψ = c/V = 1000 or V = 300 km s⁻¹.
Download figure:
Standard image High-resolution image

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** The maximum energy of electrons in various plasma environments, with the same format as Figure 2.
Download figure:
Standard image High-resolution image

The gray shaded area indicates "diffusive" acceleration, as approximated by Equation (3) and more accurately expressed as follows. The spatial diffusion coefficient is D = Vλ = ηpv/3qB, where p is the particle momentum. The particle kinetic energy is ${\varepsilon }_{{\rm{obs}}}=\sqrt{{p}^{2}{c}^{2}+{m}^{2}{c}^{4}}-m{c}^{2}$ . Then, from the condition of λ < L where λ is the spatial scale of the acceleration region, combined with ε_H = qVBL, we eventually obtain the maximum energy from the case of λ = L as

$\begin{eqnarray}&&{\varepsilon }_{{\rm{obs}}}={\varepsilon }_{{\rm{H}}}^{{\prime} }+\sqrt{{\varepsilon }_{{\rm{H}}}^{{\prime} 2}+{m}^{2}{c}^{4}}-m{c}^{2}\end{eqnarray} \tag{ 4 }$

where ${\varepsilon }_{{\rm{H}}}^{{\prime} }=(3/2\eta ){\varepsilon }_{{\rm{H}}}$ . The shaded region in Figures 2 and 3 shows this condition with the range of η = 1–100. In the limit of ε_H ≫ mc², this condition reduces to Equation (3) or ε_obs = (3/η)ε_H.

The dashed line indicates the "trapping" condition or the hard limit, as approximated by Equation (1) and more accurately expressed as follows. From the condition of r_g < L where r_g = p/qB is the gyroradius, we eventually obtain

$\begin{eqnarray}&&{\varepsilon }_{{\rm{obs}}}=\sqrt{{\psi }^{2}{\varepsilon }_{{\rm{H}}}^{2}+{m}^{2}{c}^{4}}-m{c}^{2}\end{eqnarray} \tag{ 5 }$

where ψ = c/V. The dashed curve in Figures 2 and 3 shows this condition with ψ = 1000. This corresponds to V = 300 km s⁻¹, which is close to the lower limit of the typical range of V in many plasma environments (see Table 1). In the limit of ε_H ≫ mc², Equation (5) reduces to Equation (1) or ε_obs = ψε_H = ε_limit.

It is evident from Figure 2 that the Hillas prediction generally agrees with the observed maximum energies of protons within an order of magnitude. An exception is the radiation belts at Earth and Saturn (orange), as their maximum energies substantially exceed the Hillas limit (solid line). We consider this excess is due to externally supplied CRAND protons. Their observed energies are still lower than or comparable to the hard limit (dashed line).

The Hillas prediction also agrees with the observed maximum energies of electrons within an order of magnitude, as shown in Figure 3, although there are a few large deviations. In Cygnus A, the observed maximum energy ε_obs,e is nearly 9 orders of magnitude lower than the prediction. As described in Section 3.3, we do not consider this a lower limit. For solar flares and the Crab Pulsar, the deviation is roughly ∼4 orders of magnitude, but the observed values can be regarded as a lower limit. Therefore, it is still possible that the particle energies do reach the Hillas limit in these environments.

To further evaluate the validity of the Hillas limit, we performed a multivariate analysis where the dependencies of the observed maximum energy ε_obs on the parameters V, B, and L are modeled as

$\begin{eqnarray}&&{\varepsilon }_{{\rm{obs}}}=D{V}^{a}{B}^{b}{L}^{c}\end{eqnarray} \tag{ 6 }$

where D is a constant, and the indices a, b, and c are free parameters to represent those dependencies. This is, in a sense, a purely empirical, bottom-up analysis with no prior knowledge of particle acceleration theories. Although some values of ε_obs are considered as lower limits, we took the logarithm of the above expression and performed a linear regression, treating all values equally, after removing obvious outliers: Earth's radiation belt and Saturn's radiation belt for protons and solar flares, the Crab Pulsar, and Cygnus A for electrons. The resultant best-fit model for protons is

$\begin{eqnarray}&&{\varepsilon }_{\mathrm{obs},{\rm{p}}}=(9.6\times 1{0}^{-4})\,{V}^{1.47\,\pm \,0.56}{B}^{0.79\,\pm \,0.10}{L}^{0.86\,\pm \,0.07}.\end{eqnarray} \tag{ 7 }$

It is evident that all indices are close to unity, indicating that the observed maximum energies of particles depend almost linearly on V, B, and L, consistent with Equation (2). For electrons, we performed the same analysis and obtained

$\begin{eqnarray}&&{\varepsilon }_{\mathrm{obs},{\rm{e}}}=(2.6\times 1{0}^{-2})\,{V}^{1.27\,\pm \,0.23}{B}^{0.78\,\pm \,0.09}{L}^{0.82\,\pm \,0.06}.\end{eqnarray} \tag{ 8 }$

Again, the indices are close to unity.

Therefore, we conclude that the Hillas limit actually holds in many particle acceleration environments, even when the parameters B and L vary greatly across many orders of magnitude. However, we emphasize again that we treated all values equally, including those considered lower limits. Thus, these multivariate analysis results should be validated in the future as improved measurements of the maximum energy of particles become available. Additionally, we excluded known outliers that cannot be explained by the Hillas limit alone. For such outliers, additional physical processes must be considered, such as an external supply for protons and energy loss and weak scattering conditions for electrons. These topics are discussed in the next section.

5. Discussion

While it has been established that accelerated nonthermal particles often exhibit a power-law spectrum, how far the power-law spectrum extends and what limits the particle energy have remained unclear. Here, in this study, we have confirmed earlier suggestions (K. Makishima 1999; T. Terasawa 2001) that the highest observed energy of particles ε_obs does reach the Hillas limit ε_H. In fact, they often match quite well, i.e., ε_obs ∼ ε_H, over a very wide range of parameters, despite the fact that many of the observed values should be regarded as lower limits. The multivariate analysis that treats all ε_obs values equally returned indices that are consistent with the Hillas limit expressed as Equation (2). This universal scaling law becomes even more impressive when we recall the earlier argument that the Hillas limit, which initially appears as a one-shot direction acceleration by a system-wide electric field (Equation (2)), can actually be interpreted within the scenario of stochastic and diffusive acceleration with strong scattering (T. Terasawa 2001), as expressed in Equation (3). Furthermore, the hard limit based on the condition of r_g < L also appears to nicely explain a few cases of planetary radiation belt that clearly exceeded the Hillas limit. In this section, we first explored the possible implications of this scaling law on the underlying mechanism of particle acceleration. We then discuss some of the deviations from the Hillas limit that we have identified in our systematic review of the maximum energy of particles.

5.1. Particle Acceleration Mechanisms

We emphasize that it is possible for a stochastic acceleration to work simultaneously with energization by the system-wide motional electric field. This is perhaps most clearly described by J. R. Jokipii (1987) who showed theoretically that, if the cross-field diffusion is sufficiently strong across a quasi-perpendicular shock, particles can be energized as they drift along the V × B electric field while experiencing diffusive shock acceleration, resulting in a higher rate of energy gain than in the case of a quasi-parallel shock.

It has also been suggested that both coherent and stochastic components of acceleration work for particle acceleration in MHD turbulence due to magnetic reconnection (e.g., J. Ambrosiano et al. 1988; P. Dmitruk et al. 2003). In their work, a normalized version of Equation (2) is used, following W. H. Matthaeus et al. (1984), and the system size L is replaced with the Alfvén transit time τ_A. Then, it was demonstrated with test particle simulations that, as τ_A increases, the maximum energy of particles also increases in MHD turbulence (P. Dmitruk et al. 2003).

Similarly, two-dimensional (2D), full particle simulations have shown that the highest energy of particles is matched with the prediction of qE_RΔZ where q is the particle charge, E_R is the reconnection electric field, and ΔZ is the displacement of particles in the direction of the reconnection electric field (see Figure 11 of M. Oka et al. 2010). The match was observed in simulations of multi-island coalescence where particles change their direction of motion (and thus experience scattering) at the merging point. More formal treatment of particle acceleration in quasi-2D interacting magnetic islands can be found in, for example, G. P. Zank et al. (2014). A key point here is that the same effect might take place in 3D magnetic reconnection with turbulence. Therefore, we speculate that, even in a turbulent and diffusive environment, a small number of fortuitous particles can still experience energization while moving across the system as implied by the Hillas limit.

5.2. Electrons in Radio Galaxy Lobes

While it is difficult to distinguish the two different acceleration mechanisms with our data set alone, more detailed investigations of radio galaxy lobes (including the hot-spots) might provide an important clue. This is because there is a strikingly large discrepancy between ε_obs and ε_H for electrons in Cygnus A. The gap is roughly 9 orders of magnitude. Such a large gap exists in other FR-II radio galaxies as well, and it has been an important subject of research (e.g., F. Casse & A. Marcowith 2005; A. T. Araudo et al. 2018, and references therein). It was argued that the observed maximum energy cannot be explained by synchrotron cooling even after considering magnetic field amplification at the shock front and that the hot-spots of FR-II radio galaxies are indeed poor accelerators of particles (A. T. Araudo et al. 2016, 2018).

This is in stark contrast to the heliosphere (ACR) and the Crab Nebula where particles are accelerated at a reverse shock, similar to hot-spots of radio galaxies, and yet, the maximum energies of particles reach very close to the Hillas limit. At and around the termination shock of the heliosphere and the Crab Nebula, magnetic reconnection could occur frequently due to the intrinsic presence of current sheets in the solar and pulsar winds. Such current sheets, when compressed by the termination shock, can lead to enhanced turbulence via magnetic reconnection and associated particle acceleration (e.g., J. F. Drake et al. 2010; G. P. Zank et al. 2014, 2015; Y. Lu et al. 2021). On the other hand, jets in radio galaxies do not necessarily carry current sheets, although the precise magnetic field structure within the jets is unknown. We therefore conjecture that magnetic reconnection is lacking at the termination shock of radio-galaxy hot-spots, thereby reducing the chance of particle scattering by reconnection-induced turbulence.

Further investigations of magnetic field structure in the hot-spots may lead to a better understanding of the gap between observations and the Hillas limit of the maximum energy of particles of radio lobes, and more generally, to a deeper insight into the particle acceleration mechanisms in various plasma environments.

5.3. Electrons in Solar Flares

The maximum energy of electrons in solar flares is also much smaller than what is predicted by the Hillas limit. This is particularly important because protons do reach the energy much closer to the Hillas limit prediction. It should be emphasized again that the highest-energy electrons of solar origin were observed by in-situ measurements that are limited by the sensitivity and energy coverage of the instrument. In general, measuring high-energy electrons in space is more difficult than measuring protons, due to various reasons including the effect of background noise (secondary electrons) caused by high-energy protons. Therefore, it is important to continue improving the diagnostics of high-energy electrons of solar origin. If we determine there really is a large discrepancy between the Hillas limit and the observed maximum energy of electrons, it indicates that a physical process that is unique to electrons, either acceleration or energy loss, is at play. Because electrons reach the Hillas limit energy levels in many space plasma environments and there is no reason to believe turbulence is weaker in solar plasmas than in space plasmas, we tentatively consider the possibility that the maximum energy of solar electrons is limited by radiative energy losses.

Figure 4 shows the timescale of different energy-loss processes, as a function of electron energy, obtained by using formulas summarized by M. S. Longair (2011). For Coulomb collisions, we used a formula provided by P. Saint-Hilaire & A. O. Benz (2005). It is evident that, while collisions dominate in the nonrelativistic regime, energy loss by synchrotron emission dominates in the relativistic regime (see also J. M. McTiernan & V. Petrosian 1990; and Appendix of S. Krucker et al. 2008).

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Energy-loss time for fast electrons in the solar corona, demonstrating the importance of synchrotron cooling in the relativistic regime. Four different processes (as annotated) are considered with densities n = 10⁸–10¹⁰ cm⁻³ and magnetic fields B = (5–50) × 10⁶ nT, or 50–500 G . In the case of inverse Compton scattering, we considered interactions with photons of solar origin instead of those from the cosmic microwave background.
Download figure:
Standard image High-resolution image

In fact, the timescale of synchrotron energy loss for a particle of energy ε is expressed as

$\begin{eqnarray}&&{\tau }_{\mathrm{synch}}\sim 38{\left(\displaystyle \frac{B}{100\,{\rm{G}}}\right)}^{-2}{\left(\displaystyle \frac{\varepsilon }{\mathrm{GeV}}\right)}^{-1}\,{\rm{s}}.\end{eqnarray} \tag{ 9 }$

This is shorter than the typical duration of large-scale impulsive flares (∼minutes) but longer than that of X-ray pulsations (∼0.5 s), which are thought to be individual sequences of energy-release (e.g., M. J. Aschwanden 2005).

This can also be compared to the energy-loss time due to inverse Compton scattering off photons of solar origin. With the solar luminosity of L_⊙ = 3.83 × 10³³erg s⁻¹, its energy density at the solar surface is ${L}_{\odot }/4\pi {R}_{\odot }^{2}\,c=1.31\times 1{0}^{18}\,\mathrm{eV}{m}^{-3}$ . Then, the energy-loss time is

$\begin{eqnarray}&&{\tau }_{{\rm{I}}C,\odot }\sim 7\times 1{0}^{3}{\left(\displaystyle \frac{\varepsilon }{\mathrm{GeV}}\right)}^{-1}\,{\rm{s}}.\end{eqnarray} \tag{ 10 }$

Similarly, if we consider inverse Compton scattering off CMB photons with an energy density of 2.65 × 10⁵eV m⁻³, the energy-loss time is

$\begin{eqnarray}&&{\tau }_{{\rm{I}}C,\mathrm{CMB}}\sim 4\times 1{0}^{16}{\left(\displaystyle \frac{\varepsilon }{\mathrm{GeV}}\right)}^{-1}\,{\rm{s}}.\end{eqnarray} \tag{ 11 }$

These timescales are much larger than that of energy loss by synchrotron emission (see also Figure 4), and electron energies could reach the Hillas limit if they were to lose energy only by the inverse Compton emission. Therefore, synchrotron emission remains as the possible mechanism of radiative energy loss during electron acceleration in the solar corona.

Now, if we consider the one-shot direct acceleration by the electric field E ∼ VB in the corona, where V ∼ 3000 km s⁻¹ and B ∼ 100 G, we can readily see that electrons would reach relativistic energies and travel over the distance of L ∼ 20 Mm almost instantly within ∼9 × 10⁻⁵s. Therefore, the energy loss by synchrotron is negligible in the direct acceleration scheme.

This is not necessarily the case for the diffusive acceleration scheme from which Equation (3) is derived. The acceleration timescale is ∼D/V², and the condition D/V² ≲ τ_synch gives us

$\begin{eqnarray}&&{\varepsilon }_{{\rm{diffusive}}}\lesssim \frac{100}{\sqrt{\eta }}{\left(\frac{B}{100\,{\rm{G}}}\right)}^{-\frac{1}{2}}\left(\frac{V}{3000\,{{\rm{km}}\,{\rm{s}}}^{-1}}\right)\,{\rm{GeV}}.\end{eqnarray} \tag{ 12 }$

This indicates that the electron energy will still have no problem reaching the Hillas limit, 0.05–5 TeV, depending on the parameters B and V, in the strong scattering limit of η ∼ 1. However, in a pessimistic case of B ∼ 500 G, V ∼ 1000 km s⁻¹, and η ∼ 10⁴, the electron energy might be limited to ∼150 MeV, only a few times higher than the highest energy currently observed, 45 MeV.

In summary, for solar flare electrons, the discrepancy between the observed maximum energy and the Hillas limit can be explained by the synchrotron energy loss, but the synchrotron limit depends on the parameters V, B, and L. In particular, the parameter values of V and B can vary depending on the observation methodology and the height of plasma from the solar surface. Based on Equation (12), there is still a good chance of detecting ultra-high-energy (up to ∼100 GeV) solar flare electrons that have not yet been detected. We call for further efforts to achieve better in-situ measurements of high-energy, solar flare electrons with higher sensitivity and wider energy coverage. We hope that future measurements can more clearly identify the higher-energy cutoff in the energy spectra, just like the roll-off energies seen in proton spectra during GLEs.

5.4. Protons in Planetary Radiation Belts

In contrast to the above cases of discrepancy where electron energies fall well below the Hillas limit, proton energies that substantially exceed the Hillas limit have been observed in the radiation belts of Earth and Saturn (Figure 2). As mentioned earlier, these protons can be explained by CRAND (e.g., S. F. Singer 1958a, 1958b; W. N. Hess 1959; Y. Li et al. 2023). The GCR-induced neutrons produce protons, electrons, and antineutrinos through β decay, but protons carry the largest energy and are trapped in the radiation belts. Interestingly, proton energies do not exceed the Hillas limit as much in Jupiter's radiation belt. In fact, it has already been argued that CRAND is weak in Jupiter's radiation belt (P. Kollmann et al. 2017). In Jupiter's plasma environment, particles are continuously supplied from the volcanic moon Io (e.g., C. M. Jackman et al. 2014; R. Guo & Z. Yao 2024, and references therein). Heavier ions rather than protons are the main constituent, while CRAND only provides protons and electrons. Therefore, it is not surprising that we have not seen an enhanced flux of highest-energy protons beyond the Hillas limit.

6. Conclusion

We tested the Hillas limit using recent results of observations of various plasma environments, combined with renewed, consistent definitions of the key parameters V, B, and L. We found that particles do reach the Hillas limit in many cases, confirming earlier findings (K. Makishima 1999; T. Terasawa 2001). However, there were some exceptions. Accelerated electrons in the FR-II radio galaxy Cygnus A are 9 orders of magnitude less energetic compared to the Hillas limit predictions. We speculate that this is due to possible absence of magnetic reconnection in the hot-spots and that scattering is too weak to achieve acceleration up to the Hillas limit. Similarly, in solar flares, electrons have been observed with energies up to tens of MeV, which is well below the energy predicted by the Hillas limit. We suggest that electrons with much higher energies could be detected with improved instrumentation, although there is a possibility that the energy is limited by synchrotron cooling. On the other hand, the maximum energy of protons in planetary radiation belts clearly exceeds the Hillas limit due to externally supplied CRAND protons.

It should be noted that A. Chien et al. (2023) already argued that there is a large gap between observations and the Hillas limit in many different plasma environments, where magnetic reconnection may play a major role (Type 2). In this study, we focused only on plasma environments from which we could confidently derive key parameters V, B, and L, as well as the observed maximum energy ε_obs. Consequently, our list contained many plasma environments where shocks are likely playing a major role in energy conversion (Type 1 in Figure 1 and Table 1). A more comprehensive test of the Hillas limit with a larger number of Type 2 plasma environments is left for future studies, and we anticipate more interdisciplinary discussions of the maximum energy of particles that use both theoretical and observational approaches.

Acknowledgments

We benefited from discussions with many experts from various fields. We particularly thank Hantao Ji for valuable discussions in the early stage of this study, Lyndsay Fletcher for drawing our attention to synchrotron cooling in solar plasma, and Shunsaku Nagasawa for his help in estimating the timescales of radiative losses. We also thank Drew L. Turner for providing relevant references on radiation belts, Christopher C. Chaston and Soléne Lejosne for their helpful discussions on the same topic, and Katsuaki Asano for his valuable input on the Crab Nebula flares. M.O. was supported by NASA grants 80NSSC18K1002 and 80NSSC22K0520 at UC Berkeley.