Relativistic electron acceleration at the bow shock of Jupiter and beyond

Data

Observations for this study are from the Juno spacecraft of NASA³². The energetic particle data are provided by the JEDI⁴⁷, which measures ions and electrons from about 30 keV to 1 MeV with an energy resolution of around 20%. JEDI consists of three identical sensor heads (JEDI90, JEDI180 and JEDI270) distributed around the spacecraft to optimize pitch angle coverage over a 160° × 12° field of view with an angular resolution of about 18°. The first two energy bins of JEDI used in this study are contaminated and are not included in the analysis, resulting in four energy bins covering approximately 100 keV to 1 MeV, as shown in Fig. 2. Lower-energy ion and electron observations are obtained from the JADE⁴⁸. JADE consists of two electron sensors (JADE-E) and an ion sensor (JADE-I), both measuring ions with energy per charge from 10 eV q⁻¹ to 46.5 keV q⁻¹ across 64 energy channels and electrons with energy per charge from 30 eV q⁻¹ to 32 keV q⁻¹, with a time resolution that is mode dependent and corresponds to about 2 min in the presented event. Magnetic field vector data are sourced from the Magnetic Field Investigation (MAG) instrument⁴⁹, which uses two fluxgate magnetometers to provide measurements with a temporal resolution of 1 s. All data are presented in the JSO coordinate system, a Jupiter-centred frame in which the x-axis points to the Sun, the y-axis is in the anti-direction of the orbital motion of Jupiter and the z-axis completes the right-handed system⁵⁰.

Data post-processing and density calculations

The raw instrument data were processed to generate the products used in this analysis. The JEDI energy-time spectrograms (Fig. 1c) were created by averaging data from all three sensors and all look directions. During the observation period, the instrument operated in a low-resolution mode, binning counts into six logarithmically spaced energy channels from about 30 keV to 1 MeV and into 300 s time bins. The count rates associated with the transient event, ranging from about 20 to 60 counts per second, are considered statistically significant. The electron energy efficiency correction detailed in ref. ⁵¹ was applied, although its effect is minimal in the low-radiation environment near the magnetospheric boundary of Jupiter. For JADE, proton densities were derived from JADE-I data using a numerical integration method on SPECIES=3 data⁵². Although JADE-I is not optimized for solar wind measurements⁵³, this method has been shown to be consistent with forward-modelled Maxwellian fits for similar events²⁸. The omnidirectional differential number intensities for JADE-E were calculated by averaging the observed intensities over 48 look directions, which are binned onboard in the low rate science mode of the instrument⁴⁸.

Bow shock and foreshock transient characterization

In Extended Data Fig. 1, a magnified timeseries of the foreshock transient interval (11:30–13:30 UTC) is shown. Energetic particle intensification and plasma density depletion begin at about 12:30 UTC, with a localized compression marking the trailing edge of the structure at approximately 12:50 UTC, typical features of foreshock transients^4,5,22,54,55.

To better characterize the global environment during this encounter, we use the local magnetic field conditions and the shock normal vector estimated in ref. ⁵⁶. Using this, we obtain a normal vector of [0.77, 0.45, −0.44], consistent with the duskward Juno location. The orientation of the magnetic field with respect to this normal is shown in Extended Data Fig. 1e, suggesting that the shock orientation transitions from an oblique or quasi-parallel to a quasi-perpendicular one. Specifically, during the formation and observation of the transient itself, the orientation becomes even more quasi-parallel. This shock geometry (with θ_Bn ≲ 60°) is expected to produce substantial populations of foreshock suprathermal particles^57,58,59 associated with the formation of foreshock transients^4,5,22,54.

Particle data further support this interpretation. The presence of diffuse, isotropic suprathermal ions and electrons indicates that the spacecraft is residing within the foreshock region⁶⁰. Specifically, the pitch angle distributions (PADs) of ions and especially electrons show a clear isotropic population of accelerated particles. These PADs demonstrate that particles are distributed across all pitch angles, a signature of well-scattered populations within the foreshock. This is in agreement with characteristics of accelerated electrons observed during foreshock transients at Earth^5,18,19. An illustration of the environment and associated transient is shown in Extended Data Fig. 2.

Focusing on the foreshock transient (12:30–12:50 UT), the electron PAD signature shows a progression as the transient passes through the spacecraft. This signature suggests that particles are accelerated in the approaching region, peaking within the transient and ceasing as the spacecraft exits the structure and the field rotates to a quasi-perpendicular regime after 12:50 UT. This strongly supports a local acceleration mechanism because if the source was external, energetic electrons would be observable over wider intervals. Instead, their strict localization to the transient structure implies they are generated in situ rather than being remote-sensed. Regarding the broader spatial context, based on spacecraft speed (about 4 km s⁻¹) and the interval duration, we estimate that Juno was residing approximately 1R_J upstream of the bow shock (Fig. 2, inset). This serves as an approximate estimate, as the bow shock at planetary flanks can change location rapidly. This estimate is consistent with observations at Earth, in which transients are observed at around 1−4R_E (refs. ^19,61,62).

To determine the exact geometry and scale of the observed foreshock transient, we first established its orientation using minimum variance analysis (MVA) on the magnetic field vector data in the JSO coordinate system⁶³. This technique identifies the principal axes of the variance of the magnetic field by finding the eigenvalues (λ_max ≥ λ_int ≥ λ_min) and the corresponding eigenvectors of the covariance matrix of the magnetic field over the interval containing the transient crossing. The eigenvector associated with the minimum eigenvalue (n_MVA) is interpreted as the normal direction to the boundary of the transient, assuming a quasi-planar structure. The validity of this normal was confirmed by ensuring a large ratio of the intermediate to minimum eigenvalues (λ_int/λ_min ≫ 1). With the boundary normal established, we then estimated the scale size of the transient, L, along this direction using the single-spacecraft timing method. The scale is calculated as L = |v_sw ⋅ n_MVA| × Δt, where v_sw is the upstream solar wind velocity and Δt is the measured duration of the passage of the spacecraft through the structure. Finally, the convection electric field −V × B points towards the transient sheet, which allows particles to concentrate and form the observed transient. The overall methodological approach we followed is a standardized process typically done when single spacecraft in situ observations are available^4,5,18,28. Specifically, for our case, we used a typical upstream solar wind velocity of v_sw = [400, 0, 0] km s⁻¹ in JSO coordinates, which is in agreement with estimations of velocity during that interval, and calculated the scale as L = |v_sw ⋅ n_MVA| × Δt, where Δt was taken as a 15-min duration of the passage of the spacecraft during the transient event. It should be noted that this 15-min interval, while relatively conservative, provides a realistic range of values for the spatial scale analysis (described below). The outcome of this analysis is provided in Extended Data Table 1.

Finally, a power-law fitted to the energetic tail of the JADE and JEDI data (≥10 keV) during the foreshock transient results (Fig. 2b) in a spectral index of P ≈ −1.85 ± 0.2 with the exact 95% confidence interval being determined using a non-parametric bootstrap analysis with 1,000 iterations⁶⁴. This value suggests an acceleration process with a signature and efficiency similar to that of DSA. The obtained index is well-bounded by the canonical DSA limit of P = −1.5, a feature of efficient acceleration, and is also consistent with the expected spectral softening from −1.5 (non-relativistic) towards −2 for electrons at relativistic energies at ≳1 MeV, as they are above the electron rest mass energy⁶⁵.

Spatial scales and energy limits

The maximum energy attainable by a charged particle is fundamentally constrained by the physical properties of the accelerating environment. This limit is known as the Hillas criterion³, which relates the maximum particle energy to the available potential drop across the system. For a characteristic magnetic field B and flow velocity V, the induced motional electric field creates a potential difference across a scale L that ultimately limits the maximum attainable energy of a particle, with charge q, to E_max = qBLV.

In the specific context of diffusive shock acceleration, this limit can be expressed as a confinement condition requiring that the upstream diffusion length of the particle, L_d, remains comparable to or smaller than the system size L (ref. ³⁴). The diffusion length (L_d) can be estimated through the expression L_d ≈ D/V_sh, where D is the spatial diffusion coefficient for the maximum energy, and V_sh is the velocity of the shock. For strong scattering, as is often assumed in the foreshock regions of planetary bow shocks^2,5, the scattering approaches the Bohm limit, at which the diffusion coefficient is \(D\approx \frac13vr_\rmg\), with v being the relativistic velocity of the particle and r_g its gyroradius. By equating the diffusion length to the system size (L ≈ L_d), we recover the velocity dependence inherent in the Hillas criterion: \(L\approx \frac13\fracvV_\rmshr_\rmg\).

To derive a quantitative expression for the maximum energy from this relationship, we express the velocity v of the particle and gyroradius r_g = p/(qB) in terms of its total energy E_total = E_max + mc², where p is the momentum of the particle, q is the particle charge and B is the magnetic field magnitude. The relativistic relations p = γmv and \(E_\,\rmtotal^2=(pc)^2+(mc^2)^2\) imply \(p=\sqrtE_\,\rmtotal^2-(mc^2)^2/c\) and v = pc²/E_total. Substituting these into the confinement condition gives

which leads to the quadratic equation for the particle’s total energy,

Solving this equation for the positive energy root provides the exact solution for the maximum total energy that a size-limited shock acceleration region can produce:

where the term A = 3qBLV_sh describes the properties of the accelerator. The maximum kinetic energy is then found by subtracting the rest mass energy of the particle, E_max = E_max,total − mc². In the ultrarelativistic limit (v → c and E_total ≫ mc²), this relation reduces to the simple linear form E_max ≃ αL, with α ≃ 3qBV_sh.

The calculation of energy limits in equation (3) requires three parameters: the local upstream magnetic field strength (B), the shock velocity (V_sh) and the characteristic system size (L). In this study, we use the upstream magnetic field for B and the relative velocity between the transient and the primary shock for V_sh. The acceleration region size (L) is described by the spatial extent of the transient, which is adequately equated to the precursor or foreshock region.

To model the relationship between the acceleration region size (L) and the characteristic system size (S) as shown in Fig. 3, we assumed a power-law dependency of the form L = k ⋅ S^m (ref. ²⁸). Two separate power-law models based on planetary observations are developed: a ‘typical’ model and an ‘extreme’ model. For the typical model, we used the standoff distance of each planet as the system size (S) and its typical observed acceleration region size as the associated value for L_typ. For the extreme model, we used the same standoff distances (S) but paired them with the maximum observed acceleration region sizes (L_ext). We performed a linear fit for each model in the log–log space using ordinary least squares as implemented by the statsmodels (v.0.14.4) of Python library. The resulting fits are shown in Fig. 3a.

The extension from planetary to astrophysical scales follows three points of reasoning: (1) We establish empirically that the acceleration region size (L) scales systematically with the global shock size (S) across planetary environments, where in situ observations confirm that large-scale foreshock transients are the primary acceleration sites. (2) We apply the Hillas criterion to relate this acceleration scale (L) to the maximum particle energy (E_max), a relationship validated through the presented Juno observations, previous research and recently shown to be applicable at planetary scales^5,34. (3) We combine these scaling relations under the premise that foreshock transient-driven acceleration operates universally across collisionless shocks, independent of whether the system is a planetary bow shock, protostellar jet or supernova remnant. This final step represents an extrapolation beyond what is testable with in situ experimentation and relies on theoretical expectations regarding the physics of collisionless shocks in different astrophysical contexts, further discussed below. To model the maximum particle energies for the planetary systems shown in Fig. 3(b), we set the acceleration size (L) to the maximum observed value (L_ext). For the astrophysical objects, where direct observations of L are unavailable, we used our extreme fit model to estimate a value. Specifically, we calculated the upper boundary of the 95% prediction interval for a new observation at the object’s system size (S). This interval accounts for both the uncertainty in the fitted model and the inherent variability in the data. The resulting energy range for each object, therefore, reflects the uncertainty propagated from its local shock speed and magnetic field parameters⁶⁶. This approach provides an estimate for the maximum achievable energy that is constrained by direct observation for the planets and by a robust statistical extrapolation for the astrophysical shocks.

Physical justification for astrophysical applicability

The extension of the planetary-derived scaling to astrophysical environments necessarily involves an extrapolation beyond verifiable limits, as in situ measurements of foreshock localized processes at kiloparsec distances are not feasible for the foreseeable future. However, several physical arguments support this extrapolation as a reasonable theoretical expectation, grounded in our understanding of collisionless shock physics and the properties of astrophysical environments.

First, the astrophysical shocks that we consider (the protostellar jet HH 211 and the supernova remnants SN 1987A and SN 1006) represent curved, collisionless shocks, known from numerical simulations and indirect observational evidence to exhibit notable foreshock regions in their quasi-parallel geometries^2,8. The curved geometry of these shocks ensures that quasi-parallel configurations are naturally present across substantial portions of the shock surface, similar to planetary bow shocks.

Second, although the predicted acceleration region scales of 10⁸–10¹⁰ km for these astrophysical systems cannot be directly resolved, our estimates represent conservative lower bounds. These values are consistent with independent estimates of foreshock sizes^67,68. The consistency between our extrapolated predictions and these observational constraints provides support for the validity of our scaling approach, even in the absence of direct spatial resolution of the acceleration regions themselves.

Third, the formation of large-scale foreshock transient can arise spontaneously through plasma instabilities and wave-particle interactions, whereas transient formation is also facilitated by large-scale discontinuities in the upstream medium²². These variable upstream media are an expected feature of astrophysical quasi-parallel shock environments. The turbulent interstellar medium contains magnetized filaments and structures extending to parsec scales^69,70,71. Furthermore, supernova shock waves are expected to interact with complex circumstellar structures created by the progenitor star itself. These include stellar wind-blown bubbles and previously ejected dense shells of material, which can extend to parsec scales^72,73,74. When the expanding supernova shock encounters these pre-existing structures, it can generate large-scale magnetic field discontinuities. These encounters provide ideal conditions for seeding transient phenomena and driving particle acceleration through transient processes. Finally, the very high shock speeds observed at supernova remnants allow Alfvénic Mach numbers (M_A) to reach values of the order of 10²–10³, causing particularly strong foreshock regions upstream of their quasi-parallel geometries that are highly favourable for the formation and persistence of large-scale foreshock structures. We note, however, that the maximum extent of these structures may be regulated by the coherence length of the upstream magnetic field (L_B). For remnants in the galactic disk (for example, SN 1572), in which the shock radius can exceed L_B, the acceleration region might be affected by the coherence scale, possibly creating a plateau in the acceleration scale length. This is in contrast to high-latitude systems such as SN 1006 or young SRNs such as SN7D21987A, in which L_B > S. Overall, the interplay between tangled fields and foreshock development in the regime in which S > L_B remains an open question for future investigation.

Although these arguments establish the physical plausibility of extending our planetary framework to astrophysical scales, we emphasize that this extrapolation remains tentative in the absence of direct observational confirmation. Nevertheless, the internal consistency of our framework, validated at planetary scales and yielding predictions for SN 1006 that match observed maximum energies of about 100 TeV, provides encouraging support for the underlying physical picture.

Parameters of Solar System planets

To apply our acceleration model, we first defined the physical parameter space for each environment. This space consists of the characteristic size (S), the acceleration region scale size (L), the shock speed (V_sh) and the upstream magnetic field strength (B) of the system. For the planetary environments, the system size S is the standoff distance of the bow shock, and L is the typical scale of a foreshock transient. These values were adopted from established literature^28,75, with the dataset for Jupiter supplemented by the event presented in this work. Minor refinements to the values have been made on the basis of the most recent standoff distance statistics and observations from each planetary environment^{76,77,78,79,80}. A full range for the shock speed, V_sh, was selected for all environments to account for the variety of dynamic conditions observed at interplanetary shocks. This range is based on the typical relative speeds between the solar wind plasma and transient compressive structures, as documented in several previous works, with the upper limit essentially describing typical solar wind speed^4,5,28,81. Finally, the range for the upstream magnetic field, B, at each planet was estimated using the Parker spiral model, which describes the evolution of the heliospheric magnetic field strength with increasing distance from the Sun⁸².

Parameters for HH 211

HH 211 lies at approximately 1,000 light years away and drives a narrow bipolar jet with an associated bow shock. James Webb Space Telescope measurements showed speeds of about 80–100 km s⁻¹ (figure 5 in ref. ⁷). Polarimetric ALMA/SMA observations have shown envelope magnetic fields of 40–100 μG (4–10 nT) within about 700 AU (refs. ^83,84). The projected separation from the protostar to the bow shock is about 1,000 AU (figure 1 in ref. ⁷), which we take as the effective standoff distance. To account for variability, we took 700 AU as the system size (S), and a shock speed between 50 km s⁻¹ and 150 km s⁻¹ based on the observed outer shock. The relatively cold environment upstream of these propagating jets allows the local Alfvén Mach number (M_A) to be relatively high, subsequently causing the formation of a strong foreshock precursor.

Parameters of SN 1987A

SN 1987A in the Large Magellanic Cloud has been extensively monitored through radio and X-ray observations. During the interaction with the dense equatorial ring, the forward shock decelerated to around 2,300 km s⁻¹, later re-accelerating to 3,500–3,600 km s⁻¹ (refs. ^85,86,87); for our analysis, we used a conservative range between 2,000 km s⁻¹ and 4,000 km s⁻¹. For the upstream interstellar magnetic field, we expect values of the order of 1–5 μG (0.1–0.5 nT)^88,89. Finally, for an equivalence of a standoff distance of the expanding shock (or arc), we used the radius of the blast wave of 0.1 pc, resulting in an approximately 200,000 AU system size (S).

Parameters of SN 1006

SN 1006 is a young, shell-type supernova remnant that provides compelling evidence for efficient particle acceleration to very high energies, similar to other historical remnants, such as Tycho⁹⁰, Cas A⁹¹ and RX J1713.7–3946⁹². Its shock velocity has been well-constrained by observations of its expansion, with estimates typically placing it in the range of 3,500–4,500 km s⁻¹. The strength of the upstream ambient magnetic field values is taken to be between 1 μG and 2 μG (0.1–0.2 nT). The system size of the remnant is approximately 9–10 parsecs, which translates to a system size (S) of roughly 2  × 10⁶ AU. A unique feature of SN 1006 is the strong observational evidence, derived from its X-ray synchrotron and TeV emission, that particles are being accelerated to energies of the order of 100 TeV within its shocks^8,36,38. This makes SN 1006 an ideal test case for evaluating acceleration frameworks and provides a direct benchmark for the predictions of our model.

The very high speeds observed at SN shocks allow Alfvénic Mach numbers (M_A) to reach values of the order of 10³−10⁴ causing a particularly strong foreshock upstream of their quasi-parallel geometry. This parameter range estimation is difficult to constrain as expanding SN shocks change during their lifetime. However, the parameter space above provides a reasonable order of magnitude estimate and uncertainty range for the maximum obtainable particle energy. The parameters for the planetary environments that were used for the fitting are shown in Extended Data Table 2, and all parameters used for our generalized maximum energy model are summarized in Extended Data Table 3. It should be noted that the obtained shock speed range for each planetary environment is predominantly driven by the relative speed of the solar wind with respect to the shock, whereas for astrophysical objects, it is primarily dictated by the outward shock expansion to the relatively stable interstellar medium. Finally, although the upstream magnetic field values used in this study may be elevated by compression effects within the foreshock or precursor region, the acceleration regions themselves are characterized by local magnetic depressions. Consequently, the magnetic field within these transient structures is expected to be notably weaker than in the surrounding environment.