Trapped Particle Motion in Magnetodisk Fields

The spatial and temporal characterization of trapped charged particle trajectories in magnetospheres has been extensively studied in dipole magnetic field structures. Such studies have allowed the calculation of spatial quantities, such as equatorial loss cone size as a function of radial distance, the location of the mirror points along particular field lines (L‐shells) as a function of the particle's equatorial pitch angle, and temporal quantities such as the bounce period and drift period as a function of the radial distance and the particle's pitch angle at the equator. In this study, we present analogous calculations for the disk‐like field structure associated with the giant rotation‐dominated magnetospheres of Jupiter and Saturn as described by the University College London/Achilleos‐Guio‐Arridge (UCL/AGA) magnetodisk model. We discuss the effect of the magnetodisk field on various particle parameters and make a comparison with the analogous motion in a dipole field. The bounce period in a magnetodisk field is in general smaller the larger the equatorial distance and pitch angle, by a factor as large as ∼8 for Jupiter and ∼2.5 for Saturn. Similarly, the drift period is generally smaller, by a factor as large as ∼2.2 for equatorial distances ∼20–24 RJ at Jupiter and ∼1.5 for equatorial distances ∼7–11 RS at Saturn.


Introduction
The Earth's internal magnetic field is, to a good approximation, dipolar, and charged particles in the magnetosphere can remain trapped in this field, according to their kinetic energy, pitch angle, and equatorial distance. The motion of a trapped particle is characterized by three independent timescales. From fast to slow, these are the cyclotron (gyration) period, the meridional bouncing period, and the azimuthal drift period. Since the discovery of charged particles trapped in the Earth's magnetic field (Van Allen et al., 1959), such dynamics for a dipolar field have been extensively studied (e.g., Hamlin et al., 1961;Lew, 1961;Roederer & Zhang, 2014;Walt, 2005) and widely applied to, for example, the dynamics of high-energy electron and proton populations in the van Allen radiation belts. At the gas giant planets, Jupiter and Saturn, the magnetic field deviates substantially from a dipole configuration because of the internal source of plasma provided by the moons Io and Enceladus, respectively, and the fast planetary rotation period (∼10 hr). The magnetic field is stretched into a disk-like structure near the equator where centrifugal force is largest. This structure is often referred to as a magnetodisk (e.g., Gledhill, 1967;Kivelson, 2015). The characteristics of trapped charged particle dynamics in Saturn's inner magnetosphere have been studied using an approximate dipolar field (Thomsen & Van Allen, 1980). Later, Birmingham (1982) used the models of Connerney et al. (1981aConnerney et al. ( , 1981b of the Jovian and Kronian magnetospheric magnetic field based on Voyager magnetometer observations to analyze charged particle motion in the guiding center approximation. More recently, various studies involving charged particle dynamics such as ring current modeling (Brandt et al., 1981a;Carbary et al., 2009), energetic neutral atom (ENA) dynamics (Carbary & Mitchell, 2014), energetic particle injection dynamics (Mauk et al., 2005;Paranicas et al., 2007;, and weathering process by charged particle bombardment (Nordheim et al., 2017(Nordheim et al., , 2018, rely on these kinds of calculations assuming the dipolar approximation provided by Thomsen and Van Allen (1980). A notable exception is the study of Roussos et al. (2013) who compared energetic electron microsignature drifts observed by Cassini at Saturn with their model for bounce-averaged magnetic drift based on three different nondipolar magnetic field models of Saturn. However, observations show that the magnetic field increasingly deviates from a dipole field when moving out from the inner to the middle magnetosphere.
For time variations of the magnetic field that are slow compared to the corresponding timescale of each type of motion, an adiabatic invariant is defined (Öztürk, 2012). The first invariant, B , is associated with the cyclotron motion of the particle and expresses the conservation of the magnetic flux enclosed by the particle's gyromotion with cyclotron angular frequency Ω g = qB∕m, where q and m are the charge and mass of the particle. In the more general relativistic case the mass m is replaced by the relativistic mass m 0 , where is the Lorentz factor = 1∕ √ 1 − 2 and is the ratio v∕c of the particle speed v to the speed of light in vacuum c and m 0 the particle's rest mass. We will from now on consider the relativistic case for the sake of generality. The second invariant, J, is associated with the meridional component of motion along the magnetic field between the two mirror points in each hemisphere and implies that the particle moves so as to preserve the length of the particle trajectory between the two mirror points, even in the presence of electric fields or slow time-dependent fields compared to the bouncing period. The third invariant, Φ, is associated with the particle's azimuthal drift around the magnetized planet, and it represents the conservation of the magnetic flux encompassed by the guiding drift path (or drift shell) of a particle for magnetospheric changes slow compared to the drift period. For more details on the adiabatic invariants see, for instance, Northrop and Birmingham (1982), Öztürk (2012), and Roederer and Zhang (2014).
Conservation of the first adiabatic invariant B , defined as the magnetic moment of the current I generated by the charged particle moving on its circular path, I = qΩ g ∕(2 ), with velocity v ⟂ , and therefore gyroradius implies that the quantity sin 2 ∕B, where is the pitch angle of the particle with respect to the magnetic field, remains constant. As a consequence the pitch angle becomes larger for more intense magnetic field.
In the guiding center approximation, where the particle's geometric center of the gyration motion moves along the magnetic field line, the mirror point magnetic latitude, m , is defined implicitly through the expression of the magnetic field at the mirror point, B m = B(r m , m ), that is, the location where the particle bounces back (reverses its velocity component parallel to the guiding field line) where eq is the pitch angle of the particle at its equatorial location, with radial distance R eq , and magnetic field B eq = B(R eq , = 0).
For a dipole field in the guiding center approximation, m depends solely on eq and is the solution of the equation (Hamlin et al., 1961) cos 6 m − sin 2 eq √ 1 + 3 sin 2 m = 0.
The bounce period b , and the bounce-averaged azimuthal drift period d , related to the second and third adiabatic invariants, respectively, are then expressed as integrals of the motion of the guiding center particle along the field line (Baumjohann & Treumann, 1996) where ds is an arc element along the guiding field line, v || is the particle's velocity component along the magnetic field line, and the change of longitude Δ during one bounce period b is given by the following:

10.1029/2020JA027827
where r is the radial distance to the particle and the magnetic drift velocity v D is the sum of curvature drift (v c ) and gradient drift (v g ) velocities; that is, For a particle moving in an inhomogeneous magnetic field, keeping only the first-order term B in the Taylor expansion of B about the guiding center of the particle's motion, inserting in Newton's law, and averaging over a gyroperiod leads to the following expression for the magnetic gradient drift velocity (Baumjohann & Treumann, 1996) where v g is perpendicular to both B and B. Note that retaining only the first-order term B in the Taylor expansion of B about the guiding center requires the particle's motion to be helical in the smallest scale and that the magnetic field does not change significantly within a gyroradius, that is, that r g ≪ B∕∇B.
Similarly in a curved magnetic field, the guiding center of a particle will effectively experience a centrifugal force, associated with field-aligned component of motion, leading to a general force drift with velocity where R c is the radius of curvature vector of the guiding center trajectory; that is, R c points from the center of curvature to the field line. Similarly to the calculation of the magnetic gradient drift velocity, this expression requires the radius of curvature to be much larger that the gyroradius, that is, that r g ∕R c ≪ 1.
Thus, the longitudinal change Δ during one bounce period b can be split into two contributions, curvature (Δ c ) and magnetic gradient (Δ g ) components Note that in the case of a curl-free field, that is, in absence of any currents, such as a pure dipole field, the radius of curvature R c is antiparallel to ⟂ B (i.e., R c ∕R c 2 = − ⟂ B∕B), and v c reduces to but in the general case equation (9) has to be considered to compute v c .

Generalized Formulation of Particle Motion
For a parametrization of the magnetic field line in polar coordinates, r( ) (where r is the radial distance from the planet center and the magnetic latitude), the element of arc length ds along any magnetic field line is given by ds 2 = dr 2 + r 2 d 2 and by definition: Thus, For a pure magnetic motion, where only magnetic field B exerts a force perpendicular to v, the total kinetic energy is conserved. Assuming the adiabatic invariant B is also conserved, we can write the velocity components of the particle, parallel (v || ) and perpendicular (v ⟂ ) to the field, as a function of the constant total velocity v and the values of the magnetic field at the position of the particle, B, and at the mirror point, B m : 10 The bouncing period b can be rewritten as with the dimensionless function Φ defined as wherer = r∕R P andR eq = R eq ∕R P are lengths normalized to the planetary radius R P . For a purely dipolar field,R eq corresponds to the value of the classical McIlwain L parameter or L-shell; that is, R eq is the equatorial (maximum) radial distance to which field lines on the L-shell extend. It is worth noting that Φ depends solely on the values of the magnetic field along the field line.
Also along any field line parameterized in polar coordinates r( ), the radius of curvature vector R c is given by where n is the unit normal vector, lying orthogonal to B in the plane of the field line, and the second-order derivative d 2 r∕d 2 can be expressed as a function of B r , B and their first-order derivatives with respect to using equation (12). Finally, the curvature is defined as the inverse of the norm of R c , = 1∕R c .
In a similar way to the bouncing period, the bounce-averaged longitudinal drift period d can be rewritten as with the dimensionless function Γ defined as the sum where Γ c and Γ g correspond, respectively, to the contributions from the curvature, and gradient drift motions: whereB = B∕B P andB m = B m ∕B P are normalized field strength relative to the field at the planetary surface equator B P , and ∇ r and ∇ are gradient components in polar coordinates. It is worth noting that Γ∕Φ depends on the values of the magnetic field components along the field line and also on their steepness across the field line (through the field gradient terms) and the shape of the field line (through the field curvature).
In the case of a dipole field, both bounce and bounce-averaged drift periods have been approximated by various analytic expressions. Among the most commonly used are (Baumjohann & Treumann, 1996;Hamlin et al., 1961) 10.1029/2020JA027827 where the dimensionless functions, Φ in equation (16) and Γ∕Φ in equation (19), are approximated by first-order polynomials in sin eq , andR eq has been replaced by the dipole L-shell value.
We developed a numerical framework to compute the functions Φ ( R eq , eq ) and Γ ( R eq , eq ) for any arbitrary magnetic field structure and compute their best fit to bivariate polynomials in R eq and sin eq , in order to provide approximate expressions similar to equations (23) and (24) for any arbitrary magnetic field.

Trapped Motion Properties in Jovian Magnetodisk
Our University College London/Achilleos-Guio-Arridge (UCL/AGA) magnetodisk model  uses the formalism developed in Caudal (1986) to compute axisymmetric models of the rotating Jovian and Kronian plasmadisks in which magnetic, centrifugal, and plasma pressure forces are in equilibrium. The magnetodisk model computes by an iterative method the magnetic Euler potential , which contains all the information about the poloidal magnetic field of the axisymmetric magnetodisk and is constant along the field lines. A correction is added to at each iteration, starting from the Euler potential of the initial (plasma free) dipole field. The correction decreases as the algorithm converges toward a solution and stops when the solution does not change more than a prescribed tolerance. Our model does not account for current sheet distortion known as the warping (or hinging) of the magnetodisk structure when the dipole magnetic equator is tilted with respect to the solar wind direction (Arridge et al., 2008). However, it is important to note in this context that transformation-based methods have been developed in the literature, which allow axisymmetric "flat-magnetodisk" field models to be modified for purposes of modeling the fields of asymmetrically tilted/hinged current sheets (e.g., Arridge et al., 2008;Achilleos et al., 2014;Sorba et al., 2018;Tsyganenko, 1998).
Here we use the output of our magnetic field model for a standard dayside Jovian disk configuration where the subsolar magnetopause is located at R mp = 90R J , where Jupiter equatorial radius is R J = 71,492 km, and with a hot ion population characterized by the index K h = 3 × 10 7 Pa m T −1 (see Achilleos, Guio, & Arridge, 2010, for details). This index indicates the global level of hot plasma pressure in the outer magnetosphere (product of hot plasma pressure and unit magnetic flux tube volume).
In Figure 1, we compare and quantify the difference in the geometry of the dipole and magnetodisk fields in the inner and middle magnetosphere. In the upper panel, the Euler magnetic potential , associated with the poloidal field model, is color coded in cylindrical coordinates, and field lines (contours of constant ) are labeled with an "equivalent dipole" L * parameter.
For the dipole field, the parameter L * is equal to the equatorial distance R eq of the field line in R J units, (i.e., the L-shell value). For the magnetodisk field, it is equal to the equatorial distance to which a pure dipole field line, emanating from the same ionospheric foot point (at approximately the planet's surface, i.e., R = R J ), as the labeled magnetodisk field line, would extend. Hence pure dipole and magnetodisk field lines of equal equivalent dipole L * enclose equal magnetic flux. This definition is in complete agreement with the definition of the L * invariant coordinate, a dimensionless quantity introduced first by Roederer (1970) where Φ is the magnetic flux encompassed by the guiding drift shell considered. Thus, since the UCL/AGA magnetodisk and pure dipole field models are both centered and axisymmetric, the magnetic flux Φ i integrated over the polar cap region bounded by a given ionospheric colatitude i can be used to specify a flux shell of field lines, which extend from that colatitude to some characteristic equatorial distance R eq . If the field were purely a centered dipole, we would have L * = R eq . For a dipole-plus-disk field, we have L * < R eq , where L * now corresponds to the equatorial distance of a pure dipole field line emanating from the same colatitude i (and associated with the same bounded magnetic flux Φ i , since, at the ionosphere, the current sheet field is negligible compared to that of the planetary dipole see also Lejosne (2014), for instance, Figure 1.
The lower left panel shows the equatorial distance R eq (in units of R J ) of the magnetic shell of field lines as a function of the equivalent dipole L * , for the total range of the magnetodisk model output, for the dipole (blue solid line) and the magnetodisk (green solid line). For the dipole field this simply corresponds to the line with slope unity since L * =R eq = L. For the magnetodisk we can see that the field lines remain to a good approximation dipolar for equatorial distances corresponding to L * ≲ 4, that is, where the green line does not significantly deviate from the blue line.
The lower right panel shows the equatorial distance R eq of the magnetic shell normalized to the dipole L-shell as function of the equivalent dipole L * for a range covering the inner and well into the Jovian middle magnetosphere. We can see that the magnetodisk model field lines are stretched out from dipole configuration by a factor as large as ∼3.25 (right panel) and indicated by the green line deviating from and increasing faster than the blue line (left panel). The last closed field line in the magnetodisk model output, at R eq = 90R J , corresponds (i.e., has same ionospheric anchor point) to the dipole field line with L * ∼ 45.1. For R eq ≳ 30 R J , the field line stretching does not increase as rapidly, as seen by the inflection point at L * ∼ 13.6 indicated as a vertical dashed line in the panels. This behavior is an effect of the outer boundary imposed in the model at the magnetopause within which the magnetic field is confined. For that reason we will only consider equatorial distances R eq ≲ 30 R J , well into the middle magnetosphere, and including the orbit of Ganymede at ∼15 R J , to calculate the dimensionless functions Φ and Γ∕Φ that characterize the particle's bounce and bounce-averaged drift periods. This range of distances represents a regime of purer magnetodisk structure. We aim to study the near magnetopause field topology in a future investigation.
The calculations of the functions in equations (17)-(20) were carried over the intervals 2-30 R J for R eq , and 16-72 • for eq . The minimum pitch angle value 16 • corresponds to a particle mirroring at the planet's surface (loss cone angle) while the maximum value corresponds to particles mirroring at latitudes ≲5 • . Figure 2 presents the latitude of the mirror points m defined in equation (2) and computed for the equatorial range and for a wide range of pitch angle, for both the dipole and the magnetodisk fields, from the nominal Jovian model described above (as seen in Figure 1). For equatorial distances ≲5 R J , the mirror point latitudes for both dipole and magnetodisk fields are very similar, as could have been anticipated from the similarity of the magnetic fields in Figure 1. For the dipole field, the left panel in Figure 2, the latitude of mirror point d m does not depend on R eq , as expected from equation (3). This is essentially a consequence of the self-similarity of dipole field lines of different L. For the magnetodisk field (right panel), m m is decreasing substantially as R eq increases, reflecting the stretching and confinement toward the equator of the field lines, due to the corresponding equatorial confinement of the plasma (which carries current) due to centrifugal force. The small jump seen in m m at ∼7.6 R J is a minor artifact due to a discontinuity in the UCL/AGA magnetodisk model and corresponds to the inner edge of the hot plasma distribution, clearly visible in the modeled azimuthal current density (see for instance, Achilleos, Guio, Arridge, Sergis, et al., 2010;Achilleos, 2018).
Figures 3 and 4 present the dimensionless integrals Φ and Γ∕Φ computed using equations (17) and (20)-(22), mirror latitudes shown in Figure 2, and calculated for both the dipole and the magnetodisk magnetic fields.
For the dipole field (left panel in Figures 3 and 4), there is no dependency on R eq for either quantity, as expected from equations (23) and (24). Note how small the range of variations of these quantities for the dipole are compared to the magnetodisk case; only the largest isocontour Φ = 1 is seen in Φ d , while only the smallest isocontour Γ∕Φ = 0.45 is seen in Γ d ∕Φ d . For the magnetodisk case, on the other hand, note how Φ m (right panel in Figure 3), and thus the bounce period, drops for large values of both R eq and eq . Figure 3. From left to right, the dimensionless function Φ characterizing the bounce period defined in equation (17), as function of equatorial distance and pitch angle, for the dipole field and the magnetodisk. The same color range limit is used to facilitate the comparison. Black lines correspond to isocontours of the same value of Φ, separated by 0.25 units. Quantitatively, Φ m is smaller than Φ d by a factor as large as ∼8, and the average value for Φ d ∕Φ m is ∼2.5 for the data presented in Figure 3. This behavior is due to the strong decrease of m with increasing R eq , reflecting the equatorial confinement of the plasma. In the case of the magnetodisk integral Γ m ∕Φ m (right panel in Figure 4), inversely proportional to the bounce-averaged drift period as seen equation (19), a sharp increase can be noted for R eq in the range 19-25 R J and for eq ≳ 50 • . Quantitatively, the ratio Γ m ∕Φ m ∕(Γ d ∕Φ d ) is as large as ∼2.2; therefore, the drift period for the magnetodisk is smaller than for the dipole by up to the same factor. The average value of the factor Γ m ∕Φ m ∕(Γ d ∕Φ d ) for the data presented in Figure 4 is ∼1.6. Note that the dipole and magnetodisk drift shells of the same equivalent L * will enclose similar magnetic flux (as seen previously in the discussion about magnetic shell mapping in relation to Figure 1). The differences in drift period are the result of the different azimuthal drift velocities experienced by the particle due to different guiding line geometry as can been seen from equations (7)-(9). The difference in the curvature and magnetic gradient contributions is further discussed in section 6.
Similar to the jump in m m in Figure 2, the jumps seen at ∼7.6 R J on both Φ m and Γ m ∕Φ m (right panels in Figures 3 and 4) are artifacts due to the discontinuity introduced by the inner edge of the modeled hot plasma distribution. Note also the artifact visible mostly in Φ m but also faintly in Γ m ∕Φ m as a jump at large eq , just above ∼7.6 R J and moving toward larger R eq as eq decreases. This artifact corresponds to the field line, which is conjugate to the edge of the hot plasma distribution at the equator. One can also note a very faint jump at ∼5.5-6 R J corresponding to the position of the Io torus. These features in the plasma model conspire to create a total, superposed structure that retains a couple of distinctive sharp ledges in the profile of the relevant integrals. These features can be further understood by examining the signature of this discontinuity, seen as an arc about the equator at R eq ∼ 7.6 R J , in the magnetic field gradient ∇B∕B and field curvature maps in cylindrical coordinates, in, respectively, the middle left and right panels of Figure 11 in section 6.

Analytical Approximations of and ∕
In order to provide realistic and practical formulations for magnetodisk studies, we also computed best fits of our numerical results using bivariate polynomials inR eq and sin eq to account for the magnetodisk field structure and thus obtain analytic approximation formulae similar to equations (23) and (24) for the bounce and bounce-averaged drift periods of the Jovian magnetodisk studied here. We may express b and d as where the estimates  Φ and  Γ∕Φ of the integrals Φ and Γ∕Φ are bivariate polynomials of the form Note. Also shown are the value of R 2 , the coefficient of multiple determination, and RMSE, the root-mean-square residual (see text). The indicated equatorial range R eq is the one used for the fitting.
The fitting was first validated for the dipole field seen in the left panel of Figure 1. The fitted coefficients p 00 and p 01 of equation (27) for the estimates  Φ d and  Γ d ∕Φ d of the functions are summarized in Table 1, together with their uncertainties in parentheses and a measure of the goodness of fit. The polynomial coefficients for both the approximations are in very good agreement with the ones given by equations (23) and (24). The coefficient of multiple determination R 2 , defined by equation 1 in Kvålseth (1985), is a measure of goodness of fit for regression models. It can be interpreted as the proportion of the total variance in the model (i.e., the polynomial fits) that is able to explain the variance in the functions. For the dipole we can see that more than 99.9% of the fitted model reproduces the functional values.
The fitting was then carried out for the magnetodisk field in the upper right panel of Figure 1. We started by limiting our investigation to bivariate polynomials of degree two (linear combination of the six monomials forming its basis) and considering the yet unused four monomials, that is, the linear termR eq , the bilinear termR eq sin eq , and the second-order terms sin 2 eq andR 2 eq . We found that the third most significant term in the expansion is the bilinear termR eq sin eq with coefficient p 11 . The contributions of the other terms are much smaller and do not improve significantly the goodness of fit parameters (see the discussion regarding R 2 and root-mean-square error [RMSE] as in Tables 1 and 2 below).
The fitted coefficients p 00 , p 01 , and p 11 for the estimates  Φ m and  Γ m ∕Φ m of the functions are summarized in Table 2 for two different ranges in R eq and are now discussed further.
We can see that the estimate for Φ m performs very well for both ranges of R eq as indicated by the high 95% and 98% values of R 2 , and the small 6% and 3% values of the residual RMSE. The values for the coefficients p ij 's are consistent between the two ranges. The estimate for Γ m ∕Φ m , on the other hand, does not perform as well for the large range 2-30 R J . This can be understood by the structure of Γ m ∕Φ m , which exhibits a peak around 20 R J toward large pitch angles. This structure cannot be accounted for with a polynomial of degree 2, and this result is further confirmed by the good fit achieved for the subrange 2-22 R J where the peak structure is cut away.
We continued our investigation to improve the fit for Γ m ∕Φ m over the wider equatorial range R eq = 2-30 R J and considered all the terms in a bivariate polynomial of degree 3, that is, 10 terms, and investigated the polynomials with an extra fourth term. We found that the fourth most significant term in the expansion improving the coefficient of multiple determination is the term R eq 2 sin eq with coefficient p 21 . The resulting coefficients for the fit of Γ m ∕Φ m are given in Table 3. The fourth coefficient p 21 increases substantially the value of the coefficient of multiple determination R 2 from a value of 43.5% to 73.4% and decreases by the same factor the RMSE residuals.
The coefficients in Tables 2 and 3, together with equations (25) and (26), provide new approximate formulae, valid well into the typical Jovian middle magnetosphere and including the orbit of Ganymede, for the bounce and bounce-averaged drift periods.
For a charged particle of mass m and velocity v, or equivalently with kinetic energy E and rest energy E 0 = m 0 c 2 , we write the bouncing period Jup b at Jupiter, in a manner similar to Thomsen and Van Allen (1980), and in units of seconds, as where we substituted v by c in equation (25) and used the identity = . Note that the leading constant in equation (28) is in seconds, the kinetic and rest energies, E and E 0 have to be in the same units, and the other terms in parentheses are dimensionless. A note of caution is issued here when using this approximation, as the value of the polynomial in parentheses might become negative for sufficiently large equatorial distance R eq and large pitch angle eq , a clear limitation of the approximation. It is therefore important to apply the formula within its described region of validity in ( R eq , eq ) space.
Similarly, the bounce-averaged drift period Jup d in hour units is (0.55 − 055 sin eq + 0.10R eq sin eq − 2.54 · 10 −3R2 eq sin eq ) −1 , where we substituted m 0 v 2 in equation (26), using the identity , and where Z = q∕e is the charge number, for example, Z = −1 for electrons (drifting westward in the frame of the rotating planet), and Z = 1 for protons (drifting eastward in the frame of the rotating planet). Note that the leading constant in equation (29) is in hour MeV, and the kinetic and rest energies, E and E 0 have to be expressed in MeV in this case. The strength of Jupiter's equatorial magnetic field used is B J = 428,000 nT.
As pointed out in section 1, studies that involve charged particle dynamics calculation such as ring current modeling (Brandt et al., 2008;Carbary et al., 2009), ENA dynamics (Carbary & Mitchell, 2014), energetic particle injection dynamics (Mauk et al., 2005;Paranicas et al., 2007Paranicas et al., , 2010, and weathering processes by charged particle bombardment (Nordheim et al., 2017(Nordheim et al., , 2018 would definitely benefit from the expressions for the bounce and drift period presented here, since they reflect the significant influence of more realistic nondipolar field structure.
It is also important to note thatR eq denotes the true equatorial distance in the magnetodisk normalized to R J and can be mapped to the equivalent dipole L * -shell as shown in the lower panels of Figure 1.

Trapped Motion Properties in Kronian Magnetodisk
Here we use the output of our magnetic field model for a standard Kronian disk configuration where the magnetopause is located at R mp = 25R S , where Saturn equatorial radius is R S = 60,268 km, and with a hot ion population characterized by the index K h = 2 × 10 6 Pa m T −1 (Achilleos, Guio, & Arridge, 2010).  Figure 1 for the standard Kronian disk calculated with the UCL/AGA magnetodisk model. Figure 5 shows the differences in the geometry of the dipole and the magnetodisk fields for Saturn in a similar way to Jupiter presented in Figure 1. Note how the stretching of the magnetodisk is less pronounced for Saturn than Jupiter, a factor as large as ∼1.8 for Saturn versus ∼3 for Jupiter. The last closed field line in the magnetodisk model at R eq = 25R S corresponds to the dipole field line with L * ∼ 14.7.
For R eq ≳ 15 R S , the field line stretching does not increase as rapidly, as seen by the inflection point at L * ∼ 8.9 indicated as a vertical dashed line in the panels, and is an effect of the outer boundary imposed in the model at the magnetopause within which the magnetic field is confined. For that reason we will only consider equatorial distances R eq ≲ 16 R S , well into the middle magnetosphere, including the orbit of Rhea at ∼8.74 R S , to calculate the dimensionless functions Φ and Γ∕Φ. This range represents a regime of purer magnetodisk structure as previously considered for the case for Jupiter. We also aim to study the near magnetopause field topology of Saturn in a future investigation.   Figure 3 with the dimensionless bounce integral Φ but for the Kronian system. (17)-(20) were carried over the intervals 2-16 R S for R eq and 16-72 • for eq . Similarly to Jupiter, the minimum pitch angle value 16 • corresponds to a particle mirroring at the planet's surface (loss cone angle) while the maximum value corresponds to particles mirroring at latitudes ≲5 • . Figure 6 presents the latitude of the mirror points m for Saturn similarly to the case of Jupiter in Figure 2. Note how, like Jupiter, even for small distances ∼4 R S , the latitude of mirror point of the magnetodisk m m deviates significantly from the dipole case, reflecting the stretching and confinement toward the equator of the field lines.

The calculations of the functions in equations
The small jump seen in m m at ∼8 R S is a minor artifact due to a similar discontinuity in the UCL/AGA magnetodisk model as for Jupiter and corresponds to the inner edge of the hot plasma distribution, clearly visible in the modeled azimuthal current density (see, for instance, Achilleos, Guio, Arridge, Sergis, et al. 2010;Achilleos, 2018).
As in Jupiter's case, Figures 7 and 8 present the dimensionless integrals Φ and Γ∕Φ for the dipole and the magnetodisk fields calculated with equation (17) and equations (20)-(22) and the mirror latitudes shown in Figure 6 for both the dipole and the magnetodisk. Note how Φ m in Figure 7 presents very similar characteristics to the case of Jupiter in Figure 3. In particular, the value of Φ m drops for large values of R eq and eq . This is due again to the significant decrease of m with increasing R eq , reflecting the equatorial confinement of the plasma. In the case of Saturn, though, Φ m drops by a factor as large as ∼2.5 compared to Φ d , a moderate factor compared to the factor of ∼8 for Jupiter. The average value of the ratio Φ d ∕Φ m for the data presented in Figure 7 is ∼1.5, compared to ∼2.5 for Jupiter.  Figure 4 with the dimensionless bounce-averaged drift integral Γ∕Φ but for the Kronian system. Table 4  Same Table as Table 1 Figure 5 X p X 00 p X 01 R 2 × 100 RMSE R eq ∈ 2-16R S Φ d 1.27 (7 · 10 −5 ) −0.54 (10 · 10 −5 ) 99.9 0 .0035

but for the Dipole Field Simulation of the Kronian System in
Similarly to the Jupiter case, we note that the integral Γ m ∕Φ m for the magnetodisk (right panel in Figure 8) is larger than its dipole counterpart Γ d ∕Φ d (left in same figure), meaning smaller drift period for the magnetodisk than the dipole field. In the case of Saturn, the drift period is smaller by a factor as large as ∼1.5, moderate compared to the factor of ∼2.2 for Jupiter, and an average factor ∼1.2 is found for the data in Figure 8, compared to ∼1.6 for Jupiter. It is worth noting that the broad maximum in the integral Γ m ∕Φ m for of Jupiter around ∼20-24 R J (right panel in Figure 4), is not so clear in the case of Saturn (right panel in Figure 8), due to the discontinuity artifact in the magnetodisk model around 8R S . Nevertheless, a weak local maximum can be seen for large pitch angle and around equatorial distance ∼13 R S . This distance is close to the distance at which the north-south field ΔB z , produced by the magnetodisk current, changes sign (e.g, .
Finally, we followed the same methodology introduced in section 4 for Jupiter and computed analytic approximations of Φ and Γ∕Φ for the Saturn case for the equatorial range of distances indicated.
We first validated the dipole case at Saturn (Table 4) and note the complete agreement of the coefficients p 00 and p 01 , the coefficients of multiple determination and the RMSE residuals with the Jupiter case in Table 1.
The fitted coefficients p 00 and p 01 of equation (27) for the estimates  Φ m and  Γ m ∕Φ m for the magnetodisk case are then summarized in Tables 5 and 6.
As seen at Jupiter, the fit of  Γ m ∕Φ m is poor for the wide equatorial range considered, 2-16 R S and improves by reducing the upper boundary to 12 R S as seen in Table 5.
The same method used in section 4 to improve the fit of the bounce-averaged drift integral was carried out for Saturn, and the resulting coefficients are summarized in Table 6. Note the improvement reflected by a coefficient of multiple determination of 84.9% compared to 16.5% for the total range of equatorial distance and even 58% for the reduced range.  Similarly to Jupiter, the coefficients in Tables 5 and 6 together with equations (25) and (26) provide new approximate formulae, valid well into the typical Kronian middle magnetosphere and including the orbit of Enceladus, for the bounce and drift periods of a charged particle.
For a charged particle of mass m and velocity v, or equivalently with kinetic and rest energies, E and E 0 , we can write similarly to Thomsen and Van Allen (1980), the bouncing period Sat b expressed in second units: (0.44 − 0.25 sin eq + 0.12R eq sin eq − 7.18 · 10 −3R2 eq sin eq ) −1 .
As for the case of Jupiter in equation (29), kinetic and rest energies in equation (31) have to be in MeV units. The strength of Saturn's equatorial magnetic field used is B S = 21,160 nT.
Note the similarity of order for the bounce and drift periods at Jupiter (equations (28) and (29)) and Saturn (equations (30) and (31)), especially the cross term R eq sin eq . These magnetodisk formulae, however, compared to the reference values of the dipole case, indicate a stronger deviation from dipole field for Jupiter than Saturn. This comparison indicates the differences in the magnetodisk field geometry at these planets and therefore differences in their respective ring current densities. Such differences can be traced to the differences in plasma source rate (mass loading), an order of magnitude less for Enceladus in the Kronian system compared to Io for Jupiter (Vasyliñas, 2008). But, although the plasma source rate from Enceladus at Saturn is an order of magnitude smaller (in absolute terms) than that from Io at Jupiter, suggesting the current density and thus the magnetodisk field geometry should be very different, the values of the dimensionless mass input rate (scaled to relevant planetary parameters) are more comparable (Vasyliñas, 2008).

Curvature and Gradient Drift Contribution
Finally, we examine the respective contributions of the field curvature and the magnetic field strength gradient to the total longitude change over a bounce period Δ (proportional to Γ). These longitudinal changes are respectively denoted Δ c (proportional to the integral Γ c ) and Δ g (proportional to the integral Γ g ) and were introduced in equations (6), (10), and (20) in section 1.
In Figure 9 we compare the percentage of the total drift velocity due to curvature, as a function of R eq and eq , for the dipole case (left) and the magnetodisk (right) at Jupiter. For the dipole field (left panel), the drift contribution is not a function of R eq , as expected, and for eq ≪ 45 • the curvature drift dominates as m becomes larger, while for eq ≫ 45 • the gradient drift dominates as the motion becomes more confined to  Figure 9 but for the Kronian magnetodisk field presented in Figure 5. the equator. The magnetodisk field exhibits the same behavior as the dipole for R eq ≤ 7R J , as expected (see Figure 1), but for R eq ≥ 7.6R J the curvature drift largely dominates, even at large pitch angle. This behavior arises from the larger equatorial curvature of the magnetodisk. Note that, once again, the artifact seen at ∼7.6 R J , similar to the functions Φ m and Γ m ∕Φ m , is due to a discontinuity in the UCL/AGA Jovian magnetodisk model that corresponds to the inner edge of the hot plasma distribution as discussed in the previous section.
It is quite remarkable that for R eq ≥ 20R J , and independently of the pitch angle eq , the drift velocity v D is entirely due to the curvature of the field line, implying that the drift motion is entirely driven by the curvature of the magnetic field in this region of the Jovian magnetodisk. Figure 10 presents the same quantities as Figure 9 but for the case of Saturn. As pointed out for Jupiter, the artifact seen in this case at ∼8 R S is also due to a discontinuity in the UCL/AGA magnetodisk model that corresponds to the inner edge of the hot plasma distribution at Saturn. As in the case of Jupiter, the azimuthal drift at Saturn becomes dominated by curvature drift as R eq ≥ 8R S and for larger pitch angle. But unlike Jupiter, at Saturn the regime where the drift is entirely controlled by the curvature of the field is never reached.
The generally larger slopes of the Γ m c ∕Γ m isocontours in the Jovian model reflect the more intense ring current and larger field curvature in Jupiter's magnetosphere, compared to the Saturn system. This is further illustrated in Figure 11, which shows the much greater curvature m for the field lines of the outer equatorial Jovian model (right middle panel) compared to Saturn (right lower panel).
In order to further highlight the differences in the curvature and magnetic gradient contributions to the total drift in the Jovian and Kronian magnetospheres, Figure 11 presents the magnetic gradient inverse length scale ∇B∕B (left panels), and the curvature = 1∕R c (right panels), entering in equations (8) and (9), respectively. Note that all panels have the same color scales to facilitate the comparison. Superimposed on each panel are a selection of field lines (white solid lines), and the field line associated to the discontinuity corresponding to the inner edge of the hot plasma distribution (yellow dotted line).
The two upper panels present the dipole case at Jupiter. As explained in section 1, when deriving equation (11) from equation (9) in a curl-free field, the contributions to the azimuthal drift from the magnetic gradient ∇B∕B and the curvature = 1∕R c terms are identical by definition. This is confirmed in the two upper panels. The middle and lower panels present the Jovian and Kronian magnetodisks, respectively. Note how the structure of the magnetic gradient inverse length scale is similar at Jupiter and Saturn overall. The curvature for the field lines is also of the same order for Jupiter and Saturn for distances up to ∼16 R P . As described by  and Vasyliñas (2008), even though the absolute value of the ring current is much larger at Jupiter, the normalized ring current in both systems is comparable, even slightly larger at Saturn. The normalization factor for the current density is B P ∕(R P 0 ). In the outer equatorial Jovian model, that is, for R eq ≥ 25R J , on the other hand, the curvature is much more pronounced than for Saturn. Note that Figure 11 can also be used to check the validity of the guiding center approximation for a particle with given energy, by checking that its gyroradius is, at all times, smaller than both the radius of curvature 1∕ and the gradient length scale B∕∇B. This will be the object of a separate study. Figure 11. Comparison of the magnetic gradient inverse length scale (left panels) and curvature (right panels), in normalized unit R P −1 , in cylindrical coordinates. The upper panels are for a dipole field at Jupiter, the middle panels for the Jovian magnetodisk field, and the lower panels are for the Kronian magnetodisk field. The white contours represent field lines equidistant at the equator, while the yellow dotted line represents the field line at the discontinuity seen in Figures 9 and 10.

Conclusion
We have presented a formalism to calculate the bounce and the bounce-averaged azimuthal drift periods in the guiding center approximation for an arbitrary magnetic field, and we have applied the formalism to nominal models of Jupiter and Saturn's magnetodisks generated by the UCL/AGA magnetodisk model.
We have derived, for the first time, analytic expressions for the bounce and the bounce-averaged azimuthal drift periods for the average Jovian and Kronian magnetodisk structure, analogous to expressions for a dipole field, but where additional terms in the polynomial expansion inR eq and sin eq have been introduced to account for the disc structure. These expressions, valid well into the Jovian and Kronian middle 10.1029/2020JA027827 magnetosphere, represent an improvement over the global use of a pure dipole field, which has been extensively employed in previous literature.
Further studies would be needed to check the sensitivity of the coefficients of the polynomial expansion to different configurations of the Jovian and Kronian magnetospheres (compressed and expanded states) and thus how the solar wind and suprathermal population state influence the bounce and the bounce-averaged azimuthal drift periods. Even so, the formulae presented here are still applicable for a typical field configuration.
Other useful studies would include comparison of the results of the guiding center approximation calculation in this paper against results from particle tracing simulations. In particular, the investigation of the limits to which the adiabatic invariants are conserved and thus characterization of the range of validity (in terms of particle energy, for instance) of the approximate formulae presented in this paper. In a future extension of this work, we also aim to include the effects of centrifugal force on particle motion, which are expected to be more pronounced at particle kinetic energies comparable to or smaller than the change in centrifugal potential along their trajectories.