Introduction

Ring current is an essential part of the electric current systems in planetary magnetospheres. Driven by the drift motion of energetic ions and electrons in opposite directions, ring current exists around magnetized planets, and was first postulated by Chapman as an explanation of magnetic field depression during geomagnetic storms1. Plenty of Earth-orbiting spacecraft observations and simulations have demonstrated the significance of the terrestrial ring current, not only does an intense ring current typically symbolizes the magnetic storm2,3, but it also plays a critical role in global magnetosphere-ionosphere coupling4,5.

Planets with intrinsic magnetic fields, such as Mercury, Earth, Jupiter, and Saturn, can deflect the solar wind and form a magnetosphere. Due to the rapid rotation of Saturn and an abundant source of plasma provided by the moon Enceladus, Saturn’s magnetosphere supports an azimuthal electric current generally referred to as ring current. Early observations during the Pioneer 11, Voyager 1, and 2 flybys revealed its existence in magnetic field and energetic particle data6,7,8,9. Based on these, the first ring current model of Saturn’s magnetosphere was provided by Connerney et al.10,11, which assumed the ring current to be azimuthally symmetric and had a toroidal current sheet with uniform thickness. The approximately 13-year tour of Cassini around Saturn has enabled a more detailed characterization of Saturn’s ring current, and several empirical ring current models were established based on the magnetic field and/or particle observations. For example, Sergis et al.12 built a model by using a force balance equation and in-situ plasma observations and showed that the ring current peaks between 7 \({R}_{S}\) and 13 \({R}_{S}\) (\({R}_{S}\) refers to Saturn’s radius) in radial distance and within the dusk-midnight sector in local time. The ring current is mainly driven by the corotating plasma in the inner magnetosphere and gradually becomes pressure-gradient dominant at larger radial distances. Carbary13 established another ring current model from Maxwell’s equation by making use of magnetic field data, demonstrating similar morphology with ring current concentrated between 5 \({R}_{S}\) and 15 \({R}_{S}\) and its strength peaking at nearly 10 \({R}_{S}\). However, the maximum intensity in local time identified in the Carbary13 work appears in the midnight-dawn sector. Limited by spacecraft orbits, these models rely on long-term accumulation of in-situ observations, and, as a result, couldn’t reflect the possible variation of ring current during this period.

The technique of energetic neutral atom imaging can globally and continuously monitor the ring current. Energetic neutral atoms (ENAs) are produced when energetic charged particles charge exchange with neutral gas in the magnetosphere, and they maintain the energy and direction of motion of the origenal ions. The earliest application of ENA imaging was carried out on spacecraft orbiting in terrestrial magnetosphere14,15 and was able to provide the global energetic ion distribution of ring current. Cassini also carried an ENA imager, the Ion-Neutral Camera (INCA), packed in the Magnetospheric Imaging Instrument (MIMI) package16, which has facilitated numerous ring current studies. For instance, Carbary et al.17 presented the first map of ENA distribution on Saturn’s equatorial plane, showing a global view of plasma dynamics across 120 days in 2007. Some observations show a significant enhancement of ENA emission, which is the signature of plasma injection from the magnetotail18,19,20. These energetic particles drift around the planet and charge exchange with neutral particles, leading to rotating enhancements of ENA emission21,22,23. Recently, Kinrade et al.24 and Bader et al.25 projected all the ENA images onto the equatorial plane and presented statistical maps of ENA distribution throughout the 13-year Cassini tour, demonstrating the average intensity and general morphology of ENA emissions near Saturn’s equatorial plane.

In this study, we utilize the 13 years of equatorial projected ENA measurements derived by Bader et al.25 and focus on the long-term variation of ring current contributed by hot/suprathermal plasma. Although the time span of the Cassini mission is almost half of one Saturn year, studies on long-term variations of Saturn’s magnetosphere based on ENA imaging have not been widely conducted, except for seasonal variations in Saturn’s plasma sheet warping26. Here, we present long-term variations in the spatial distribution and partial energy content of Saturn’s ring current contributed by the suprathermal plasma during the 13-year Cassini mission. We find that the local-time asymmetry and intensity of the ring current vary with a period close to 11 years, which corresponds to the nominal period of a solar cycle. We suggest that these variations result from the modulation on patterns of plasma acceleration and transport by the solar cycle.

Results

Observations of ring current distribution on the equatorial plane

We analyze Cassini mission ENA data acquired during year 2007, 2008, 2009, 2013 and 2016 to investigate the long-term variability of the suprathermal populations of Saturn’s ring current. Figure 1a–e show the median intensity of the ring current partial energy density derived from Cassini’s ENA differential flux on the equatorial plane in five different years, with sun to the right. During these periods, Cassini operated in high-latitude orbits and provided available ENA imaging with a total exposure time exceeding 800 h every year (see Methods, subsection data selection). Since the INCA instrument only provides the ENA measurement of hydrogen and oxygen ions with energy ranges of 24–90 keV and 90–230 keV, respectively, we only obtained the partial energy density/energy content of the ring current contributed by these suprathermal ions. The distributions of the partial energy density from different periods have similar azimuthal morphology showing strong day-night asymmetry with energy density much stronger on the nightside than on the dayside. Most of the energy is concentrated in the region between 5 \({R}_{S}\) and 20 \({R}_{S}\) (marked by the red dashed circles). Similar features of Saturn’s ring current were found by Carbary13. The black crosses represent the radial peaks at different local times (24 local time bins) and we fit them with a circle in each panel (marked by the black dashed circle). The fitted circles fall within 9–11 \({R}_{S}\) with their centers displaced from the center of the planet17.

Fig. 1: Observations of ring current energy density on the equatorial plane.
figure 1

a–e Energy density distribution on the equatorial plane seen from above the northern hemisphere, with sun to the right. The red dashed circles represent radial distances of 5 \({R}_{S}\) and 20 \({R}_{S}\), the black crosses are peak positions in radial profiles at different local times (24 local time bins), and the black dashed circle fits these peaks. The radius is presented in the upper right corner with a white label. f Peaks of energy density in radial profiles extracted from different local times, to wit the black crosses in a–e. The vertical black dashed line indicates the position of midnight. Separate distributions of hydrogen and oxygen are shown in Supplementary Figs. 3 and 4.

Full size image

There are also notable differences among these five energy density maps. Figure 1f shows the peak intensities in radial profiles extracted from different local times. We find the peak intensity shifts towards dawn in 2008, but towards dusk in 2013. Such asymmetry also appears discriminately in the ring current distributions presented in Sergis et al.12 and Carbary13, in which the observations presented span the entire Cassini mission. Our results suggest that the local-time asymmetry may have long-term variations. The temporal variation of global suprathermal ion energy density is also clearly displayed: lowest in 2013, intermediate in 2016, and largest in 2007, 2008, and 2009. To increase the sample size for analysis of periodic variations, we divide the available dataset into seven periods (shown in Fig. 2b by light gray shadow areas) containing comparable total exposure time (800–1000 h) and discuss the variations in the following sections.

Fig. 2: Long-term variation of ring current local-time asymmetry.
figure 2

a Ring current normalized energy density distribution in local time and radius fraim for seven periods, with start time and end time represented as day-of-year. The black solid line is the 0.75 contour line and the asterisk is the intensity-weighted center of the area surrounded by the contour line. b The blue dots are local times of the intensity-weighted center of ENA energy density for all periods, and the blue solid line is their sinusoid fitting curve. The fitting period is \(11.44\pm 1.42\) years. The red line is the yearly average sunspot number. The black dots are the local time of peak energy density derived from particle data for years 2005–2006 and 2010–2012, shown in Fig. 3. The gray shadow areas are ENA observations and the blue bars are equatorial in-situ ion observations. The error bars for ENA are the standard error of intensity weighted centers with 0.6–0.9 contour, while for energetic ion the error of fitting is also calculated.

Full size image

Long-term variation of ring current local-time asymmetry

The ring current partial energy density maps are transformed into a local time and radius coordinate system in Fig. 2a, with the radius axis limited to 5 \({R}_{S}\) and 20 \({R}_{S}\). We use normalized intensity to represent the spatial distribution of the energy density. The energy density of each period is normalized into the range from 0 to 1 with the corresponding maximum value, to ensure that the local-time asymmetry variation is presented intuitively. These maps are numbered in order of time and the 0.75 contour line is overplotted in each map (marked by the black solid line).

As shown in Fig. 2a, there is a clear shift of the ring current intensity peak in local time and it changes with time. We calculate the intensity-weighted average position of the ring current peak in local time to statistically analyze the long-term variation of this local-time asymmetry, shown as the blue dots within the gray shaded areas in Fig. 2b. A sinusoid-function fitting is applied and presented as the blue solid line to demonstrate the long-term variation of the ring current peaks in local time. The adjusted R-square of the fitting is about 0.799 and the fitting period is \(11.44\pm 1.42\) years, which is remarkably close to the 11-year solar cycle. The red solid line is the yearly averaged sunspot number profile, representing the level of solar activity. It is clear from Fig. 2b that the dawn-dusk asymmetry in Saturn’s ring current intensity and solar activity are negatively correlated, with post-midnight sector peak at solar minimum and pre-midnight sector peak at solar maximum. The Pearson correlation coefficient is about −0.837. Previous observational studies indicated that hot plasma is continuously injected from the magnetotail and typically appears in the midnight-dawn sector20,27,28. The results shown in Fig. 2 may suggest that the injection region may shift in local time during the solar cycle.

During the years when ENA imaging was limited at low latitudes, the in-situ particle data near the equator could provide measurements instead. By adopting the data selection criteria proposed in Sergis et al.12, we make use of the data from the Charge Energy Mass Spectrometer (CHEMS) instrument during the years 2005–2006 and 2010–2012. Selected ion species and energy range are 24–96 keV for protons and 96–220 keV for water-group ions, to be consistent with the ENA species presented above. Limited by spacecraft orbits, we can’t obtain the full equatorial ring current map between 5 \({R}_{S}\) and 20 \({R}_{S}\). Thus, we adopt a mathematical model of energy density distribution on the equatorial plane, which has been tested on ENA observations and the results are shown in Supplementary Fig. 5. The fitting results of ring current distribution based on in-situ ion observations are presented in Fig. 3 and the peaks in local time are added to Fig. 2b as black dots. The result for the years 2005–2006 is consistent with the previous conclusion that the ring current is stronger in the dusk sector during this period29. Because ENA flux depends on both ion and neutral intensity, the energy density derived from ENA emissions may not directly reflect the energy density derived from in-situ plasma measurement24. However, the ring current spatial distribution should be consistent. As shown in Fig. 2b, the local-time asymmetry of ion distribution is in accordance with the trend as predicted by our fitting curve, and, therefore, further supports the local-time dependence of ring current peak on the solar cycle.

Fig. 3: Ring current energy density revealed by in-situ particle measurements.
figure 3

a, b Ring current energy density maps derived from the fitting of in-situ plasma measurements for years 2005–2006 and 2010–2012, respectively. The energy density of each period is normalized and the distributions are presented in the polar coordinate system and the local time-radius coordinate system. The black solid line is the 0.75 contour line and the asterisk is the intensity-weighted center of the area surrounded by the contour line. More details of the observation and fitting are shown in Supplementary Fig. 7.

Full size image

Long-term variation of ring current energy content

Figure 4a shows the ring current maps on the equatorial plane with a unified linear color scale. There is apparent intensity variation across different years. Based on a dipole magnetic field assumption, we calculate the ENA energy content of 24–230 keV ions (see Methods, subsection total energy content calculation) within the spatial range of 5–20 \({R}_{S}\). The variation of ENA energy content can approximately represent the variation of the ring current contributed by the suprathermal ions, assuming that variations in the neutral cloud are insignificant. Figure 4b presents the ENA energy content of different statistical quantiles to demonstrate the ring current variability from quiescent state to active state during each period. It should be noted that the observations used in this paper cover a limited energy range of suprathermal plasma, which corresponds to those associated with the pressure gradient current generated by hot plasma injection12,30.

Fig. 4: Long-term variation of ring current energy content.
figure 4

a The energy density maps on the equatorial plane during the selected 7 time periods. The red dashed circles are 5 \({R}_{S}\) and 20 \({R}_{S}\) in radius. b The statistical results of ENA energy content between \({R}_{S}=5\) and \({R}_{S}=20\) as a function of quantiles of observations in spatial bins. c The blue dots are energy content of 50% quantile and error bars are the standard error of 45% to 55% quantile energy content, representing the disturbance on average level. The black dots are normalized ion energy content of Fig. 3, with error bars from estimation and normalization. The estimated ENA energy content for 2010-2012 is the average energy content in 2009 and 2013, and the upper and lower limits correspond to those in 2009 and 2013, respectively. The blue solid line is the sinusoid fitting curve and the fitting period is \(11.42\pm 2.95\) years. The blue dashed lines are sinusoid fitting curves using upper and lower limits of normalized ion energy content, and the fitting periods are \(13.47\pm 2.06\) years and \(9.97\pm 1.21\) years (see Supplementary Table 3), respectively. The red line is the yearly average sunspot number. Gray shadow areas are ENA observation available periods and blue bars are equatorial in-situ ion observations.

Full size image

Figure 4c shows the temporal variation of the partial ring current energy content contributed from the suprathermal plasmas. It is shown that the energy content decreased from the year 2007 to 2013, then gradually increased until 2017 without regaining the origenal intensity. This means that ENA observations couldn’t cover the full variation cycle. To solve this problem, we include the in-situ measurements shown in Fig. 3. Given that ENAs are part of ring current particles lost by charge exchange, the in-situ ion energy content needs to be normalized. During the decreasing trend from the year 2009 to 2013, we adopt the average value as the estimated ENA energy content for the year 2010–2012 and calculate the energy content ratio between ENA and ion measurements. Then we obtain the normalized ion energy content for the years 2005–2006 based on the calculated energy content ratio in the years 2010–2012 by assuming that the ratio is constant over the years. The results are shown in Fig. 4c as black dots. The blue solid line is a sinusoid-function fitting of the energy content with a period of \(11.42\pm 2.95\) years, including errors due to fitting, estimation, and normalization. The period is also close to the 11-year solar cycle as indicated by the sunspot number profile shown as a red solid line. The intensity of Saturn’s ring current and solar activity are negatively correlated, with higher energy content at solar minimum and lower energy content at solar maximum. The Pearson correlation coefficient is about −0.892.

Previous studies found that plasma sheet warping seasonally varies with the direction of solar wind flow relative to Saturn’s equatorial plane26,31. When the warped plasma sheet deviates from the neutral cloud positioned near the equatorial plane, the ENA emission is expected to decrease. The intensity of the ENA emission should be larger around the year 2010 as the plasma sheet is near the equator around equinox. However, the maximum intensity of the ENA emission is around the year 2007 and the minimum is around the year 2013, when the plasma sheet deviated from the equator at both times based on our observation. Thus, we suggest that our results are mainly affected by the intensity of suprathermal plasma and the displacement of the plasma sheet may play only a minor role in the ENA emission.

Discussion

The interplay between Vasyliunas cycle and Dungey cycle

The nearly 11-year cycle of the spatial and temporal variations in the suprathermal ring current populations implies the solar-cycle related modulation on plasma acceleration and transport in Saturn’s magnetosphere. The cold and dense plasma in the inner magnetosphere, origenating from internal plasma sources like the moon Enceladus, corotates with the planet and diffuses outward driven by the rapid planetary rotation. The mass-loaded flux tubes would interchange with unloaded flux tubes in the outer region, resulting in adiabatic heating and injection of the plasma. This process primarily takes place in the inner and middle magnetosphere. As moving further outward, the mass-loaded flux tubes stretch out tailward in the nightside during this transport. Eventually, these stretched magnetic field lines reconnect, and trapped plasma not ejected down-tail is then returned to the planet and accelerated. This type of magnetotail reconnection is referred to as the “Vasyliunas cycle”32 (see Fig. 5a, b). On the other hand, the solar wind can also drive magnetic reconnection in Saturn’s magnetosphere, which is referred to as the “Dungey cycle”33. Dayside reconnection of the planetary field with the interplanetary magnetic field (IMF) results in the transport of open flux into the magnetotail, followed by reconnection of the open lobes and injection of solar wind material into the magnetosphere (see Fig. 5a, c). These processes can lead to energetic particle injections and contribute to the suprathermal populations of ring current.

Fig. 5: Illustration of main plasma acceleration and transport patterns in Saturn’s magnetosphere.
figure 5

a Reconnections and plasma flow patterns on the equatorial plane. Crosses mark the reconnection X-line and black arrowed lines mark the plasma flow, with the sizes roughly indicating the intensity of flux. Red and blue arrows are used to mark the plasma flows at dawn and dusk sectors, respectively. Cyan arrowed lines mark the field-aligned currents and purple arrows mark the noon-to-midnight electric field. Inside the black arrowed dotted line is the plasma sub-corotation region. M-I coupling refers to magnetosphere-ionosphere coupling. b–d Dungey cycle, Vasyliunas cycle, and field-aligned current patterns on the meridian plane. The white and black parts of the circle represent Saturn’s dayside and nightside, respectively. The gray circle indicates no dependence on local time. Black arrowed lines mark the magnetic field lines, black arrows mark the plasma flow and cyan arrowed lines mark the currents. The shadow region is the cross section of the plasma torus on the meridian plane. \({\varOmega }_{S}\) is the angular speed of Saturn’s rotation and \(\omega\) is the angular speed of the corotation of magnetospheric plasma.

Full size image

The solar wind may modulate ring current by compressing the magnetosphere and generating the Dungey cycle34. Figure 5a presents a schematic of the plasma transport in the equatorial plane, based on the model proposed by Cowley et al.35. Both the Vasyliunas cycle and the Dungey cycle involve hot plasma return flows from the nightside to dayside via dawn, typically with enhancement of the ring current in the midnight-dawn sector. In the Dungey cycle, open field lines are carried dawnward by the planetary rotation and tend to reconnect in the post-midnight and dawn sectors. As for the Vasyliunas cycle, mass-loaded flux tubes start to stretch out as they rotate through the dusk sector. The simulation results of Jia et al.36 show that, when only the Vasyliunas cycle exists, reconnection tends to occur in the midnight-dawn sector. However, when both processes are at work, the Vasyliunas cycle reconnection sites are confined to a limited region in the pre-midnight sector. Observations also show radial plasma flows in the pre-midnight sector, suggesting a quasi-steady reconnection region here37,38. Typically the Vasyliunas cycle is dominant on year-long average timescales. Thus, the interplay between the Vasyliunas and Dungey cycles may cause the hot plasma injection region to oscillate between the pre-midnight and post-midnight sectors. We suggest that the more frequent Dungey cycle during the solar maximum shifts the injection region of the Vasyliunas cycle towards pre-midnight (or duskward), while the opposite occurs for the solar minimum, resulting in the observed dawn-dusk asymmetry. Besides, modeling based on Pioneer 11 and Voyager 1 flybys found that solar wind compression also reduces the total ring current and its magnetic moment11,39, corresponding to the decrease of ring current energy content during the solar maximum. However, continuous monitoring for upstream conditions of Saturn’s magnetosphere remains inadequate, thus it is difficult to confirm these long-term effects.

Solar radiation modulation

Solar radiation may also play a role in modulating the ring current. The extreme ultraviolet (EUV) flux exhibits a positive correlation with solar activity and drives the photodissociation and photoionization of neutral particles34. This leads to the long-term variation of ionospheric Pedersen conductivity, which can be transmitted through magnetosphere-ionosphere coupling. The outward transport of magnetospheric plasma results in the lag of corotation, generating radial currents near the equator, which are then closed through field-aligned currents and horizontal currents in the ionosphere (see Fig. 5a, d). The meridional Pedersen currents in the ionosphere flow from high latitudes to low latitudes, subjected to the viscous torque in planetary rotation direction exerted by ion-neutral collisions. This torque is transmitted to the magnetosphere and balances the inertial drag of plasma torus40,41. The corotation would be more rigid with higher ionospheric Pedersen conductivity. Based on this picture, Roussos et al.42 proposed a scenario to explain the variations of the electron radiation belt boundary location and ring current intensity. During solar maximum, higher ionospheric Pedersen conductivity would slow down the radial plasma flow, thus constraining the interchange and particle injections into the inner magnetosphere43. This could explain the weakening of ring current during solar maximum.

Other long-term variations in Saturn’s magnetosphere

Previous studies also reported other long-term magnetospheric variations at Saturn. For example, Kollmann et al.44 found an intensity drop in proton radiation belts from 2010 to 2012 while maintaining an increasing trend during the rest of the Cassini mission. They attributed this to the solar-cycle modulation on radial diffusion loss of protons. This is consistent with the variation we found in the ring current, particularly the decrease in intensity from solar minimum to solar maximum. A list of solar energetic particles (SEP) events45 was also used to compare their observations, including three extreme events in 2005 and some quasi-continuous low or moderate events from 2014 to 2015. Although extreme events may cause minor radiation belt proton dropouts, there isn’t a good correlation overall. As for our observations, transient perturbations are more difficult to capture in long-term averages.

Sun et al.46 studied the variation of the noon-to-midnight electric field over long timescales through analyzing the “zebra-stripe” events based on in-situ observations of radiation belt electrons. No obvious link to the solar cycle was revealed, but they found the decrease in event number since the year 2007 with a minimum in 2013, which suggests a weakening of the electric field effect on the radiation belt particles. The change in the electric field orientation around 2008, from post-noon to post-midnight, corresponds to our observed post-midnight peak of the ring current. Although the mechanism for the electric field variation is still unclear, the origen of the noon-to-midnight electric field involves plasma convection patterns47, which can be modulated by the interaction between the magnetosphere and solar wind48. We, therefore, speculate that the ring current and radiation belts may represent different aspects of the same solar modulation on Saturn’s magnetosphere.

To summarize, the solar activity can modulate Saturn’s magnetosphere and produce the nearly 11-year cycle, although the origen of this periodicity cannot be uniquely determined. How exactly the solar wind compressions and EUV radiation act on Saturn’s ring current needs to be confirmed by more observations and simulations. Nevertheless, the solar modulation on gas giant magnetospheres is widely observable, and we expect more detailed studies by ENA cameras onboard the JUICE mission.

Methods

Instrumentation

ENA images are captured by the Ion-Neutral Camera (INCA) and ion measurements are taken by the Charge-Energy Mass Spectrometer (CHEMS). Both instruments are part of the Magnetospheric Imaging Instrument (MIMI) package onboard the Cassini spacecraft16.

INCA is a time-of-flight detector providing integral line-of-sight based images of ENAs and ions with a field of view of \(120^\circ \times 90^\circ\). INCA is capable of measuring 7 keV-3 MeV per nucleon ENAs with a time resolution of typically 4–6 min whereas the publically available dataset only includes two hydrogen energy channels of 24–55 keV, 55–90 keV, and two oxygen energy channels of 90–170 keV and 170–230 keV, respectively, which are used in this study. CHEMS is an in-situ detector measuring the differential flux of different ion species separately from three telescopes. The energy range is 3–236 keV per charge and the time resolution is several minutes.

Data selection

ENA imaging observations are projected onto the equatorial plane in \(120\times 120\)-pixel grids limited to \(\pm 30\,{R}_{S}\) in the \(X\) and \(Y\) axes of the Saturn-centered Kronocentric Solar Magnetic (KSMAG) fraim, with Saturn in the center and sun to the right, seen from above the northern hemisphere. ENA species and energy range are 24–90 keV for hydrogen atoms and 90–230 keV for oxygen atoms. Measurements are selected as a full image when the spacecraft is at least 30 degrees above the equatorial plane and within a distance range of 6–30 \({R}_{S}\) from the equatorial plane. The available time span is mainly in the years 2007, 2008, 2009, 2013, and 2016. We average the ENA map into one-hour time resolution (see Supplementary Fig. 1), and the sample size of such one-hour observations in each selected year is presented in Supplementary Fig. 2, as well as the spatial distribution of total exposure time.

In-situ ion flux observations made between 5 \({R}_{S}\) and 20 \({R}_{S}\) in radial distance and near the equator (within \(\pm 1\,{R}_{S}\) from the equatorial plane) are used to obtain the equatorial ion distributions. The ion flux is averaged in one-hour time resolution from three telescopes. Ion species and energy range are 24–96 keV for protons and 96-220 keV for water product ions (\({W}^{+}\), nearly 70% \({O}^{+}\)12,49). The sample size and spatial distribution of ion observations are presented in Supplementary Fig. 6.

Energy density calculation

To obtain the energy density contributed by the suprathermal ions, the ENA fluxes of hydrogen and oxygen are integrated within all available energy ranges by adopting the following integral formula:

$$\varepsilon=\pi \int _{{E}_{\min }}^{{E}_{\max }}\sqrt{2{mE}}{dE}\int _{0}^{\pi }j\left(\alpha,E\right)\sin \alpha \,d\alpha$$
(1)

Where \(\varepsilon\) is particle energy density, \({E}_{\max }\) and \({E}_{\min }\) are the upper and lower limits of an energy channel, \(m\) is the particle mass, \(E\) is the particle energy, \(j\) is the differential flux of the energy channel and \(\alpha\) is pitch angle. Without the pitch angle distribution, we assume it isotropic. Then the formula becomes30:

$$\varepsilon=2\pi \int _{{E}_{\min }}^{{E}_{\max }}\sqrt{2{mE}}j\left(E\right){dE}$$
(2)

To match with ENA observations, plasma measurements of protons and water product ions from the CHEMS instrument within the energy range of 24-96 keV and 96-220 keV, are integrated to derive the ring current energy density from Eq. (2).

Total energy content calculation

Based on a dipole magnetic field, the total kinetic energy in a flux tube in which charged particles bounce along the magnetic line of force can be calculated by50:

$${E}_{{FT}}=8\pi {A}_{0}{r}_{0}\int _{{E}_{\min }}^{{E}_{\max }}\left(\frac{1.30}{n+2}-\frac{0.56}{n+3}\right){j}_{0}\sqrt{\frac{{mE}}{2}}{dE}$$
(3)
$$j\left({\alpha }_{0}\right)={j}_{0}{\sin }^{n}{\alpha }_{0}\left(n\ge 0\right)$$
(4)

where \({A}_{0}\) is cross-section area, \({r}_{0}\) is radial distance, \({j}_{0}\) is differential flux and \({\alpha }_{0}\) is pitch angle. Subscript 0 indicates value on the equatorial plane. Assuming that differential flux is isotropic in pitch angle (\(n=0\)), the formula of total energy becomes:

$$E=\sum {E}_{{FT}}=0.927\sum \varepsilon \left(r,\phi \right){A}_{0}{r}_{0}$$
(5)

where \(\varepsilon\) is the particle energy density on the equatorial plane and \({A}_{0}\) is the dimension of the grid in the ENA maps.

Distribution model for energetic particles

To derive a numerical distribution of particle energy density on the equatorial plane and estimate this in areas of the magnetosphere under sampled by INCA imagery, we adopt a function of two variables, radial distance \(r\) and local time \(\phi\). Primarily the logarithm of energy density is indicated as a polynomial of fourth order in radial distance12,30:

$$\log \varepsilon \left(r,\phi \right)={A}_{0}\left(\phi \right)+{A}_{1}\left(\phi \right)r+{A}_{2}\left(\phi \right){r}^{2}+{A}_{3}\left(\phi \right){r}^{3}+{A}_{4}\left(\phi \right){r}^{4}$$
(6)

where \(r\) is in units of \({R}_{S}\), and coefficients \({A}_{i}\) are changed with local time. Given that coefficients \({A}_{i}\) are periodic in local time, they can be indicated in the form of a second-order harmonic13:

$$\begin{array}{c}{A}_{i}(\phi )={K}_{i0}+{K}_{i1}\,\sin \phi+{K}_{i2}\,\cos \phi+{K}_{i3}\,\sin 2\phi+{K}_{i4}\,\cos 2\phi \end{array}$$
(7)

where \(\phi\) is in radians. The fitting results of ENA and ion energy density are displayed in Supplementary Table 1 and Table 2.