Adiabatic Invariants Calculations for Cluster Mission: A Long‐Term Product for Radiation Belts Studies

The Cluster mission has produced a large data set of electron flux measurements in the Earth's magnetosphere since its launch in late 2000. Electron fluxes are measured using Research with Adaptive Particle Imaging Detector (RAPID)/Imaging Electron Spectrometer (IES) detector as a function of energy, pitch angle, spacecraft position, and time. However, no adiabatic invariants have been calculated for Cluster so far. In this paper we present a step‐by‐step guide to calculations of adiabatic invariants and conversion of the electron flux to phase space density (PSD) in these coordinates. The electron flux is measured in two RAPID/IES energy channels providing pitch angle distribution at energies 39.2–50.5 and 68.1–94.5 keV in nominal mode since 2004. A fitting method allows to expand the conversion of the differential fluxes to the range from 40 to 150 keV. Best data coverage for phase space density in adiabatic invariant coordinates can be obtained for values of second adiabatic invariant, K , ∼ 10 2 , and values of the first adiabatic invariant μ in the range ≈ 5–20 MeV/G. Furthermore, we describe the production of a new data product “LSTAR,” equivalent to the third adiabatic invariant, available through the Cluster Science Archive for years 2001–2018 with 1‐min resolution. The produced data set adds to the availability of observations in Earth's radiation belts region and can be used for long‐term statistical purposes.


Introduction
In 1958, U.S. spacecraft Explorer-1 (e.g., Williams, 1960) was launched into the orbit, carrying onboard the experiment originally intended for measuring the cosmic rays. The measurements were susceptible to the energetic electrons in the regions that later became known as Van Allen radiation belts (Hartley & Denton, 2014). There are proton and electron radiation belts. Electron radiation belts exhibit a two-zone structure: the inner belt lies from L = 1 to ∼2.5, the outer belt is located between L values from ∼3 to 7 with maximum electron fluxes at L ∼4-5 (Lyons et al., 1972). Two electron radiation belts are separated with a so-called slot region situated within L values from ∼2 to 3 (Baker et al., 2018). The slot region is formed by a balance between pitch angle scattering loss to the atmosphere and inward radial diffusion (Kavanagh et al., 2018;Lyons & Thorne, 1973) and is usually devoid of energetic electrons. Protons form a single radiation belt with maximum flux intensities between L values 3 to 4 (Ganushkina et al., 2011).
The radiation belts can pose significant hazard to satellites. Relativistic particles with energies >1 MeV can penetrate through satellite shielding and burn the equipment onboard down (e.g., Koller et al., 2009). Moreover, particles at smaller energies of ≈50-100 keV (a so-called seed population) can accumulate at the surface of the satellite thus generating internal currents that can also lead to satellite loss (e.g., Fennell et al., 2000). Although most of the spacecraft flying in the radiation belts are affected by the energetic particles, relatively few satellites provide measurements of the radiation belts population (Koller et al., 2009). Among them, one of the longest ongoing missions is the Cluster mission launched in 2000 (Escoubet et al., 1997). The mission consists of the four identical spacecraft, each carrying onboard the Research with Adaptive Particle Imaging Detector (RAPID)/Imaging Electron Spectrometer (IES) detector measuring electron fluxes at energies 30-400 keV. It should be noted that the IES detector was not specifically designed for radiation belts studies and as such is subject to background contamination by the high-energy electrons (>400 keV) and inner-zone protons. However, the recently developed background correction techniques In this paper we present a step-by-step guide for the calculation of the adiabatic invariants and conversion of the electron fluxes to PSD for Cluster electron data. Among these we describe the details of the new data product developed for the Cluster Science Archive-the LSTAR-available with 1-min resolution for more than 18 years of the Cluster measurements. The structure of the paper takes the form of four parts, including this introductory section. Section 2 is concerned with data set and the methodology of adiabatic invariants calculations and conversion of electron fluxes to PSD. Section 3 discusses the importance of the developed product and outlines general features of PSD behaviour in 2004-2018. The conclusions are drawn in the final section.

Adiabatic Invariants
As defined in classical mechanics, under the changing conditions of motion one can define the quantities, called adiabatic invariants, which remain constant when changes are slow (Landau & Lifshitz, 1976). Since the timescales of charged particles motion in the Earth's magnetosphere are much smaller than those of the geomagnetic field change, adiabatic invariants can be introduced in order to understand a variety of physical processes. Charged particles in the inner magnetosphere undergo three types of periodic motion, and each type is associated with an adiabatic invariant. The first type, cyclotron motion around the magnetic field line, conserves the first invariant , which is expressed as follows: where p 2 ⟂ is the electron momentum orthogonal to the magnetic field direction and m 0 is the electron's rest mass and B is the magnitude of the magnetic field.
This invariant is associated with the magnetic flux through the surface covering the particle's Larmor orbit. In addition, converging of the geomagnetic field toward the poles makes the particles oscillate between two mirror points. This second type of motion is associated with an invariant J. In the general case it includes both parallel particle momentum and magnetic field geometry through the following equation: where the integral is taken along the field line between two mirror points s m and s ′ m , B(s) is the magnetic field at point s, and B m is the magnetic field at the mirror point. In the equation (2) the term outside of the integral is constant. Therefore, one can write another form of the second invariant, depending only on magnetic field and not the particle characteristics, as follows : Furthermore, an increase in distance between the mirror points leads to a decrease in particle momentum parallel to the magnetic field and particle energy, and vice versa, when magnetic mirrors converge the particles are accelerated to higher energies. The third invariant Φ is related to the azimuthal drift Journal of Geophysical Research: Space Physics 10.1029/2019JA027576 around the Earth and represents the conservation of magnetic flux through the closed drift shell of a particle. The third invariant is conserved on the timescales of hours in case of weakly changing magnetic field .
Adiabatic invariants are approximate integrals of motion (Somov B. V., 2013). In practice, it is often more convenient to use some equivalent proxies to estimate the values of the invariants. For instance, when no parallel forces act on a particle, the momentum is conserved along a bounce path (Green & Kivelson, 2004). One can therefore define an integral invariant coordinate I, which is a proxy of the second invariant (Roederer, 1970): I has the units of distance (km or R E ) (Konstantinidis & Sarris, 2015).
The Roederer L parameter, or L ⋆ , is a dimensionless proxy of the third adiabatic invariant and is defined as follows: where M is the Earth's dipole magnetic moment. The L ⋆ parameter shows the radial distance to the equatorial plane where the electron would remain if all external magnetic fields were removed (Green & Kivelson, 2004).
PSD as a function of adiabatic invariants can be obtained from satellite measurements of particle flux through the relation: In order to analyze PSD at fixed values of the adiabatic invariants, a number of steps involving interpolation should be performed. The technique used for such conversion applicable to the RAPID/IES electron observations is described in detail below.

Data Set
Launched in 2000, the Cluster mission with RAPID/IES detectors onboard has provided a large data set of electron flux observations in Earth's radiation belts and ring current. The Cluster constellation follows the polar elliptical orbit, systematically passing through radiation belts region (Escoubet et al., 1997). It should be noted that Cluster orbit has been evolving over time due to the lunisolar perturbations. Starting with the initial perigee at around 4 R E and apogee at 19.6 R E , the perigee was gradually declining up to 1.2 R E in 2010 and then rising again up to ≈6 R E in 2019. Therefore, Cluster has mainly been passing through the outer belt but traveled through the inner belt as well in 2007-2013 (for further details on the spatial coverage refer to Smirnov, Kronberg, Latallerie, et al., 2019).
Electrons with energies 30-400 keV are measured by the RAPID/IES instrument aboard each of the four satellites. IES is a solid-state silicon detector consisting of three pinhole systems. Three of these detectors are combined in a configuration providing electron flux observations over a 180 • range (Wilken et al., 1997). The energy channels of the IES detector are given in Table 1. Electron distributions are delivered in the nominal telemetry mode (NM, active most of the operation time) and burst mode (BM, allows 4 times more data but only for short periods of several hours per week) . The omnidirectional electron fluxes are provided in the nominal mode for all six energy channels ("ESPCT6" product) since the start of measurements, whereas the pitch angle distribution was provided at the beginning only in burst mode ("E3DD" product). Starting in April 2004 an additional electron product became available in the nominal mode. This product, called "PAD_L3DD," was added after reprogramming the onboard DPU and provides the electron pitch angle distribution for channels 1 and 3 for nine pitch angles, see details in Kronberg and Daly (2019). We note that for our purposes the "E3DD" product does not have enough temporal coverage, and thus, for conversion to PSD we will further after use the "L3DD" product.

IRBEM Library
A variety of tools has been developed for the adiabatic invariants calculations. Among them is the physics-based LANLGeoMag providing L ⋆ values for the Van Allen Probes constellation in real time (Mauk et al., 2012). In addition, the LANLGeoMag computes several other magnetospheric parameters by obtaining a numerical solution to the equation of motion with high accuracy (Min et al., 2013). Koller et al. (2009) used a neural network to predict L ⋆ values based on the statistical training of the network. In the present paper we use the International Radiation Belt Environment Modeling (IRBEM/ONERA) library (Boscher et al., 2012), relying on the physics-based principles. The library is fully open access and allows computation of magnetic coordinates and fields as well as geographic and heliospheric coordinate transformations and atmospheric models runs. The IRBEM-lib is written in Fortran but can be called using the IDL and MATLAB interfaces. Recently, the Python interface has been developed allowing to employ the IRBEM-lib functionality through the open access high-level language (Shumko et al., 2018).
Several of the IRBEM functions can be used to compute the adiabatic invariants. Among others, there is a "make_lstar" function. It allows to calculate magnetic coordinates at any satellite position (Boscher et al., 2012). In Python interface, the function takes as an input: time in decimal format; specifications of the magnetic field model, output parameters that should be returned and coordinate system; satellite position and the magnetic field input depending on the selected model. We note that since the "make_lstar" function does not require pitch angles as input it neglects the slight dependency of L ⋆ on local pitch angle. In order to fully account for the latter, a "make_lstar_shell_splitting" function can be employed. It requires the same input as "make_lstar" with addition of local pitch angles. The outputs include the L and L * values, magnetic field strength at the mirror point (B mir ) and at position of satellite (B(s)), magnetic local time and the I value (see equation (4)). Therefore, by calling this function, one can obtain the values of the third invariant L ⋆ as a direct output and also the values of the second adiabatic invariant K through the relation

Conversion of Electron Fluxes to PSD
PSD analysis under fixed values of adiabatic invariants has notable advantages over the conventional electron flux studies. In particular, it helps to separate non-adiabatic from adiabatic effects and allows to employ observations from different magnetospheric regions . Conversion of the measured electron differential flux (E, , r, t) as a function of energy, pitch angle, position, and time to PSD ( , K, L ⋆ , t) as a function of three adiabatic invariants and time is done in the following steps described below. The described calculations were implemented as a set of Python and MATLAB routines and are available in open access at . Figure 1 gives the schematics of the PSD calculations. We first perform the L ⋆ calculations, as this parameter is necessary for the background contamination removal.

Preliminary Steps
(i) Calculations of L ⋆ parameter. L ⋆ coordinate is equivalent to the third adiabatic invariant. The values of L ⋆ with 1-min resolution can be obtained through the Cluster Science Archive (https://csa.esac.esa.int/ csa-web/; "LSTAR" product). For the "LSTAR" product, calculations were performed using the IRBEM library ('make_lstar' function). As mentioned above, for calculations of L ⋆ it is necessary to specify the satellite position, magnetic field model (since the fluxgate magnetometer onboard gives the magnetic field value at only one point) and geomagnetic/magnetospheric conditions. Tsyganenko-89 magnetospheric model, T89, (Tsyganenko, 1989) with Kp index as an input was employed. (ii) Background correction. Since the RAPID/IES electron flux measurements are subject to contamination in the radiation belts region, it is necessary to perform a background data correction to avoid  Kronberg et al., 2016) Channel erroneous values during the next steps. In the present paper we use the empirical contamination coefficients described in detail in Kronberg et al. (2016) and Smirnov, Kronberg, Latallerie, et al. (2019). Kronberg et al. (2016) derived empirical percentages of contamination for different L ⋆ values on all IES energy channels. The contamination coefficients for Channels 1 and 3 are shown in Table 2. In order to filter the IES electron fluxes, we subtract the part of intensity attributed to the contamination using the following equation: where I measured is the total measured electron flux intensity, p is the percent of contamination from the Table 2, and I electrons is the filtered electron intensity. (iii) Calculations of K. K is one of the forms of the second adiabatic invariant. It depends on the configuration of the magnetic field and is independent of the particle's energy, charge, and momentum (as can be seen from equation (3)). Cluster data resolves nine pitch angles from 10 • to 170 • with a 20 • step. Since K values are symmetric with respect to 90 • pitch angle, we compute K for five pitch angles from 10 • to 90 • via the "make_lstar_shell_splitting" subroutine implemented in the IRBEM/ONERA library. As input parameters we use local pitch angles at which to compute K, satellite position, and Kp index as a magnetic input parameter required by the T89 model. The output of the "make_lstar_shell_splitting" subroutine includes the I values (see equation (4)) and magnetic field value at the mirror point. We then obtain K through the equation (7). We note that there are different units commonly used during K calculations, such as G 1∕2 km and nT 1∕2 R e (Hartley & Denton, 2014). In this paper we use the latter units.  will be smaller and the particle motion will be closer to parallel to the magnetic field. Therefore, one needs to set up the specific values of K and before the following steps.
Step 1. Fit pitch angles to K to obtain k . The second adiabatic invariant K is independent of the momentum and depends only on the satellite position, local pitch angles, and magnetic field geometry. To each pitch angle there exists a corresponding value of K. As discussed above, the K values are symmetric with respect to 90 • pitch angle and for the fit one can use the pitch angles from 10 • to 90 • . We perform  a linear fit using the interp function implemented in the NumPy library (Oliphant, 2006). The illustration of such linear fit is given in Figure 1.
Step 2. Fit pitch angles to flux to obtain full pitch angle distribution. In order to retrieve flux values at a specific k it is necessary to compute the full pitch angle distribution (PAD). This is done by fitting a functional form, relating pitch angles to the observed flux. The standard approach uses the following relation: and seeks optimal values of A and n by the least squares fit. An example of such a fit is given in Figure 2. It should be mentioned that this method can be used efficiently to fit pancake, field-aligned and flattop distributions, but it cannot adequately reproduce butterfly PAD's. However, the electron PAD's observed by the Cluster/IES instrument in radiation belts region predominantly have the pancake shape. An example of the weekly averaged PAD for IES Channel 3 is shown in Figure 3.
Step 3. Perform the energy fit of differential fluxes using the relation̂(E) = a exp(bE). IES measurements are available for Channels 1 and 3; therefore, we only have two points for such a fit. However, for the IES electron fluxes we adopt the same distributions as in AE9 model (Ginet et al., 2013), where the a and b coefficients remain quasi-constant for all L * values in energy range ≈40-150 keV. In logarithmic scale it is equivalent to fitting a straight line to the data, which can be done in a unique way using two data points: wherêis the fitted value of flux, E i is the energy at which we want to estimate the flux, E 1 and E 3 are the effective energies of Channels 1 and 3, respectively (all terms are in logarithmic scale). With a step of 10 keV, we obtain the differential flux values using the equation (10). The example of such fit is shown in Figure 4 based on the full 3-D distribution taken from "E3DD" product in 2003. From Figure 4 it is evident that the proposed fit gives realistic estimates of flux for energies 40-150 keV.
Step 4. Calculate the first invariant and PSD values for all energies E i = {40, 50, … , 150} keV using the following formulas: where E i is energy in MeV, m 0 c 2 is the rest energy of the particle (for electrons equal to 0.511 MeV), B is the magnetic field measured by the satellite in nT, and is the local pitch angle; where (E i ) is PSD in Geospace Environment Modeling units (c/MeV/cm) 3 for the corresponding energy,̂(E i ) is electron differential flux from the equation (10) in units cm −2 ·sr −1 ·s −1 ·keV −1 . We note that since the Cluster mission has a polar orbit, it can travel through the dielectric cusp region, where the physical processes are totally different from those in radiation belts. The L-star approximation is based on the T89 model which is computationally efficient but may not account for local field irregularities. In order to avoid data from the cusp region, we put a limit on the magnetic field: It should be more than 80 nT.
Step 5. Interpolate the PSD values onto the target value of . Such interpolation can be performed, for example, using the Python function "numpy.interp" in the following way: PSD_target=numpy.interp (mu_target,mu,PSD) where lengths of mu and PSD arrays equal the number of E i . It is a common practice to use logarithmic values of PSD and during interpolation. A visual scheme of such interpolation is given in Figure 1. Fitting the differential flux versus energy should be performed in the energy range from 40 to 150 keV. Thus, one gets the value of PSD ( , K, L * , t) as a function of adiabatic invariants and time. The percent of data values for different in range 5-250 MeV/G and K in range 0.1-1258 nT 1∕2 R e for all Cluster satellites is shown in Figure 5. For all four Cluster spacecraft the most data values can be obtained when lies in range ≈5-20 MeV/G and for higher values of K (∼100-1,000). Examples of PSD merged over all four Cluster spacecraft for the most probable values of K and are given in Figure 6.

Discussion
The Cluster mission has been measuring the pitch angle resolved electron flux in the nominal mode since 2004, thus covering approximately 1.5 solar cycles: the declining phase of Solar Cycle 23 (2004-2008) and Solar Cycle 24 (2008. Therefore, the RAPID/IES electron measurements converted to the adiabatic invariants coordinates can be used for long-term statistical analysis of Earth's radiation belts and ring current dynamics. Figures 6a-6c show the PSD evolution along the declining phase of Solar Cycle 23, the historical minimum of solar activity in 2008-2009 and during the following solar cycle 24 for several values of the first and second adiabatic invariants. Figure 6d gives the monthly smoothed AE and F 10.7 indices. The F 10.7 index is a proxy of the solar emissions from chromosphere and corona, thermal gyroresonance over sunspots, and nonthermal emissions (Tapping, 2013). The F 10.7 index follows the general pattern of the solar cycle variation exhibiting peaks during solar cycle maxima. The auroral electrojet (AE) index is calculated as a difference between the upper and lower envelopes of the magnetic variations in H component from 12 geomagnetic observatories located at 61-70 • latitude in the Northern Hemisphere (Kamide & Akasofu, 1983). The AE index is a proxy of the substorm activity and has peaks during the declining phases of the Solar Cycles 23 and 24. We use the AE index due to the fact that the Cluster PSD coverage shown in Figure 6 is within L ⋆ from ∼4 to 9, which when mapped to the ground coincide with the latitudes at which the AE index is measured. Furthermore, Smirnov, Kronberg, Latallerie, et al. (2019) showed the correlation between the electron fluxes measured by the IES detector at energies ∼40-400 keV and the AE index and solar wind dynamic pressure, while no significant correlation was found between electron fluxes and the F 10.7 . This is due to the fact that the F 10.7 follows the sunspot cycle while the other parameters follow the occurrence of coronal holes and high-speed streams, which occur more often during the declining phase causing the radiation belts enhancements.
It has long been established that the substorm activity plays an important role in radiation belts electron transport and acceleration (e.g., Boyd et al., 2016;Jaynes et al., 2015;Zhao et al., 2017). The source population electrons with energies of tens of keV are injected into the inner magnetosphere during substorms and produce waves (e.g., the whistler mode chorus) (Boyd et al., 2016). The so-called seed population electrons with energies of up to a few hundred keV are also injected during the substorms and can be accelerated to relativistic energies by wave activity (e.g., Baker et al., 1998). Zhao et al. (2017) showed that the AL index had good correlation with PSD enhancements for energetic electrons, suggesting the direct transport of lowelectrons from substorm injections. From Figure 6 it is evident that at L ⋆ ∼4-6 the PSD follows the pattern close to that of the AE index. The PSD values are generally higher during the descending phase of the Solar Cycle 23 in 2004-2005 and exhibit minimum values during the time of the so-called "radiation belts desert" in 2009-2010, which also corresponds to the intense dip in the AE index. The largest PSD values at these L ⋆ can be observed during the descending phase of the Solar Cycle 24 in 2015-2017. Li et al. (2010) analyzed the global distribution of PSD for low-electrons and showed that higher PSD values correspond to the higher AE values (e.g., their Figure 9a), which goes well with our findings shown in Figure 6. Furthermore, Li et al. (2010) reported that the radial PSD gradient was positive for in range 0.05-2 MeV/G, suggesting the radial diffusion caused by the magnetospheric activity. Our results are in good agreement with findings of Li et al. (2010), since the PSD values in Figure 6 generally increase toward higher L ⋆ .

Conclusions
Using the step-by-step guide to the adiabatic invariants calculations presented in this paper, the large Cluster data set covering almost two solar cycles can be used for long-term statistical studies in terms of electron PSD in the corresponding coordinates. Furthermore, since the Cluster constellation consists of four spacecraft the data can be useful for PSD gradients analysis.
Using the IRBEM library, the L * coordinate was calculated with 1-min resolution for the Cluster mission in 2001-2018. The created product called "LSTAR" is available through the Cluster Science Archive (https://csa.esac.esa.int/csa-web/) as one of the "Auxiliary, MAARBLE and ECLAT support data" files. The developed product has numerous scientific applications in the field of radiation belts research, such as, for instance, the more precise satellite tracking in terms of L * coordinate compared to the simplified L-shell parameter. Cluster data converted to adiabatic invariants coordinates can be used for physical models validation with data, model set up in the form of initial or boundary conditions (e.g., for the nowcasting and forecasting models), and for data assimilative models. Furthermore, the Cluster data set can be used to study injections in terms of PSD. The developed product may be complementary to other data sets and may help more accurately reconstruct the radial profiles of PSD for different energies and locations. 501100001656). The authors are grateful to David Hartley (University of Iowa) and Mykhaylo Shumko (Montana State University) for valuable comments. The Cluster data can be found at CSA Archive (http:// www.cosmos.esa.int/web/csa/). The data of AE and F 10.7 indices were taken from OMNIWeb database (omniweb.gsfc.nasa.gov).