A Self‐Consistent Model of Radial Transport in the Magnetodisks of Gas Giants Including Interhemispheric Asymmetries

Outward transport of plasma in the inner and middle magnetospheres of gas giants results from an interplay between mass loading from the inner dominant mass sources (volcanic moons), flux tube interchange in the centrifugally unstable magnetospheric plasma disk, turbulent heating of the plasma, and coupling between the equatorial plasma and the planetary upper atmosphere through magnetic field‐aligned current loops and/or Alfvén waves. We present a new analytical formalism describing large scale transport in gas giant systems, combining two historical approaches: radial diffusion of mass and energy through flux tube interchange, and angular momentum transport through corotation enforcement. Under the hypotheses of axisymmetry, steady‐state, and multi‐fluid plasma, we provide transport equations for total contents of flux tubes. They feature new transport parameters accounting for the latitudinal extent of the disk, and self‐consistently include field‐aligned potential drops in the magnetosphere‐ionosphere coupling. Our general formalism has a wealth of applications, two of which are presented, corresponding to the cases of the two gas giants: the effect of interhemispheric asymmetries in the resistive and magnetic properties between the northern and southern ionospheres on the transport of angular momentum at Jupiter, and the influence of the plasma disk thickness on transport at Saturn. We apply our formalism to derive ionospheric parameters and reproduce the Juno and Cassini data. Further work will allow for more complete numerical solutions of our equations, with the aim of capturing the broad complexity of fast rotating magnetospheric systems which can be found inside and outside the Solar System.


Introduction
The dynamics of gas giant magnetospheres in the Solar System is dominated by the interplay between the fast planetary rotation and internal plasma sources: Io in the case of Jupiter, Enceladus in the case of Saturn.Both generate a torus along and around their orbit where most of the magnetospheric plasma is produced by ionization of the neutral gas.The fast rotation of the planet drives the plasma outwards, resulting in the formation of a disk of plasma extending from the inner to the outermost regions of the magnetosphere: the magnetodisk.In such rotation-dominated magnetospheres, understanding the large-scale characteristics of plasma distributions first requires a correct description of the radial transport of the different charged particle species from their regions of generation (the moon tori) to the regions where they can be lost (generally in the outer magnetosphere and through the magnetospheric boundaries).
Building on the important findings of Voyager, Galileo, and Cassini, several lines of modeling provided important advances in this description of radial transport, mainly through two different approaches.These two approaches, denoted by "radial diffusion" and "corotation enforcement" models, are briefly described below.
Radial diffusion models focus on an accurate description of the radial transport of the plasma content in mass and energy of the flux tubes in the magnetosphere.Their description of radial transport is based on the flux tube interchange instability, studied for instance by Southwood and Kivelson (1987) and Ferrière et al. (2001).In this instability, a denser flux tube located closer to the planet is dragged outward under the influence of the centrifugal force and exchanges with a lighter flux tube located at a larger radius, experiencing a smaller centrifugal force, as detailed in the review papers of Achilleos et al. (2015) and Krupp et al. (2004).Multiple evidence of this kind of interchange motions have been reported by Bolton et al. (1997), Kivelson et al. (1997) at Jupiter and Azari et al. (2018), Rymer et al. (2009) at Saturn.It results in a net outward transport of plasma mass and energy in the magnetodisk.Radial diffusion models initially used diffusion equations to describe the radial transport of the integrated mass and energy contents of magnetic flux tubes (Saur, 2004;Siscoe & Summers, 1981).More recent models achieved a better agreement with observations by introducing some contribution of advective transport (Ng et al. (2018) at Jupiter and Neupane et al. (2021) at Saturn).The successful use of an advection-diffusion model by Ng et al. (2021) suggests that radial transport is more likely driven by advection rather than by diffusion at larger radii.Although they fully integrate the thickness of the disk and a physical description of flux tube plasma contents, radial diffusion models do not take into account exchanges of angular momentum between ionospheres and magnetodisks, as corotation enforcement models do.
Corotation enforcement models mainly focus on reproducing the rotation curve of the Jovian and Kronian magnetospheric ion populations as well as the main aurora generated in the upper atmosphere.They relate directly the location of the main auroral emissions to the regions of upward field-aligned currents corresponding to downward electron acceleration and precipitation.Following the early works of Hill (1979) and Pontius (1997), which successfully explained the rotation curves observed by the Voyager missions, Cowley and Bunce (2001), later followed by Cowley et al. (2005); Nichols and Cowley (2005); Nichols (2011b); Ray et al. (2010Ray et al. ( , 2012)), developed progressively more complex analytical models of angular momentum transport in gas giant systems.Electric current loops connecting latitudinal horizontal currents in the two conjugate ionospheres to radial currents in the equatorial magnetosphere play a key role in these models: in the quasi-static approximation they assume that the divergence of equatorial radial electric currents is balanced by currents flowing along magnetic field lines and closing in the planet's ionosphere.Observational evidence of such a large scale field-aligned current system have been found in both the ionosphere and the magnetosphere of Jupiter by, for example, Kotsiaros et al. (2019); Wang et al. (2021); Nichols (2011b); Al Saati et al. (2022); Kamran et al. (2022); Liu et al. (2023) based on Juno in-situ measurements.Mauk and Saur (2007) measured a spatial structure of less than 20 km at high latitude for those systems.The resulting DC current loops transfer momentum between the equatorial magnetosphere and the ionosphere-thermosphere via the Lorentz force.This family of models has been very successful in reproducing the rotation curves of the magnetospheric plasma as well as the average location of the main auroral emissions at Jupiter (e.g., Cowley and Bunce (2001)).While they assume interchange motions of magnetic flux tubes operating at small to medium scales in order to transport plasma radially while recirculating magnetic flux, these models do not explicitly describe these small-scale motions, contrary to radial diffusion models.
During the last 20 years, new multi-instrument spacecraft observations of the Jovian and Kronian magnetospheres became available.At Saturn, the Cassini mission performed a three-dimensional survey of the multi-phase nature of its magnetosphere and allowed a detailed description of the variation with radial distance of the thickness, composition and energetics of its magnetodisk.At Jupiter, detailed multi-instrument observations of the highlatitude and polar magnetosphere have been performed by the Juno mission, providing for the first time simultaneous diagnostics of magnetic fields, auroral emissions, particle precipitations, plasma flows, as well as fieldaligned and ionospheric current systems.This unique suite of observations provided ground truth to test the predictions of corotation enforcement models (e.g., Mauk et al. (2020), Al Saati et al. (2022), Wang et al. (2022)).But it also revealed the importance of significant interhemispheric asymmetries in all these quantities, for example, the field-aligned current systems and magnetic field at Jupiter (Al Saati et al., 2022;Connerney et al., 2022;Kotsiaros et al., 2019) and the rotation rate, periodicity, and asymmetrical coupling linked to seasonal Finally, Cassini and Juno observations provided quantitative estimates of the importance of departures from azimuthal symmetry of the magnetodisks induced by longitude and/or local time variations due to internal effects, such as the planetary magnetic field (Connerney et al., 2018(Connerney et al., , 2022)), and to external ones, mainly magnetospheric interactions with the solar wind (Arridge et al., 2011;Kita et al., 2019;Yao et al., 2022).
Building on the results of these missions, time is ripe to develop a generation of models to which can be assigned two objectives: (a) describe the large-scale transport of mass, angular momentum, and energy in one single formalism, and connect this large-scale transport to the smaller-scale turbulent transport of these same quantities; (b) integrate into these new models some of the lessons learned from Juno and Cassini observations that depart from pre-existing model predictions.
In the following, we first build on the two modeling approaches previously described to integrate them into a single theoretical framework.We then use the resulting integrated model to explore two categories of solutions: at Jupiter, we specifically address the effects of different sources of interhemispheric asymmetries on radial transport, current systems, and magnetodisk-upper atmosphere coupling; at Saturn, we include the effects of the latitudinal structure of the magnetodisk on the transport of mass and angular momentum.While we will keep the assumption of azimuthal symmetry in the new developments described in this study, we will discuss in our conclusion section the limitations that this simplifying assumption induces on the application domain of our model.
In Section 2, we introduce the different regions of the coupled planetary system described by our model, the different schemes and equation sets describing transport in each region, and the main underlying assumptions used for simplification.In Section 3, we use these schemes in order to provide explicit equations for the transport of conserved quantities in the magnetosphere, as well as for the calculation of current systems in each region.In Section 4, we couple transport and current systems between these regions in order to derive global solutions to the transport of mass, angular momentum, and energy in the system, and we explore these solutions at Jupiter and Saturn.In Section 5, we summarize the main features and outputs of the new model, discuss its limitations, and provide some perspectives for future works.

Description of the System
The system of interest for our study comprises the magnetized planet, its atmosphere, and its space environment, up to its magnetospheric boundaries with the solar wind.It can be separated into four main domains.Domain I, the planet's deep atmosphere and interior, is a key element in the system workings.Although not directly coupled to the other three domains by current systems, it fills two important functions.First, it includes the "planetary dynamo" region where the large-scale planetary magnetic field is generated and determines the distribution of magnetic flux at the base of the upper atmosphere.Second, it is through the deep atmosphere that rotation of the planet's interior is transmitted to the upper atmosphere, via a complex and still poorly understood interplay of viscous coupling, atmospheric waves, and circulation cells (Banks & Kockarts, 1973;Brooks et al., 2019).Thus Domain I is the one that imposes boundary conditions on magnetic flux and horizontal neutral atmosphere circulation at the base of the upper atmosphere, that is, Domain II.Domain II is classically divided into its neutral component, the thermosphere, and an ionized component, the ionosphere, generated by three main ionization sources: EUV and X-ray solar radiation, polar and auroral precipitating electrons and ions, and meteoric fluxes (see Nakamura et al. (2022) for a recent comprehensive description of these three sources at Jupiter).Domain II extends from the top of the mesosphere at a few microbars up to a few thousand kilometers above the 1-bar atmospheric level in the cases of Jupiter and Saturn.This domain is populated by a magnetized, weakly ionized gas in which plasma dynamics and energetics result from the interplay between collisions between neutral and charged particles, on one hand, and charged particle gyrating motions in the magnetic field, on the other hand.In the narrow altitude region (a few hundred kilometers) where collision frequencies between charged particles and neutrals approximately match charged particle gyrofrequencies, this interplay produces different mean motions for ions and electrons under the effect of electric fields and neutral winds (Banks & Kockarts, 1973;Millward et al., 2005).As a result, the ionosphere is an "Ohmic" medium, in which electric currents flow is proportional to the applied electric field, as will be described in Section 2.2 below.Domain III, immediately above, covers high latitude magnetic field lines, along which specific non-thermal transport processes take place.This domain is the critical link between the upper atmosphere (Domain II) and the magnetospheric domain (Domain IV).The complexity of its physics and transport processes, still poorly understood, has been illustrated at Jupiter by in situ observations by Juno (Lorch et al., 2022;Mauk et al., 2017Mauk et al., , 2020)).In this study we are interested in describing its electrodynamic properties, that is, the relationship existing between plasma acceleration and electric currents along field lines and field-aligned electrostatic potential differences.Our simplified approach to this difficult problem is described in Section 2.3.
Finally, Domain IV is the equatorial magnetosphere dominated by its magnetodisk, where plasma populations trapped on planetary magnetic field lines reside.For simplicity, we will describe the relationship between plasma Figure 1.A simplified sketch of plasma transport in the magnetospheric systems of Jupiter and Saturn.A neutral torus (in deep blue), present in the inner part of the magnetosphere along the orbits of Io and to a lesser extent of Europa for Jupiter, and Enceladus for Saturn, acts as a major source of plasma for their magnetospheres, mainly via electron-impact ionization.This plasma source feeds a radially extended magnetodisk (green area), through radial outward transport.The magnetosphere can thus be separated along radial distances into three different regions: a source region, a region of radial transport driven by electrodynamic coupling to the planetary fast rotation, and a loss region, in which the plasma is ejected from the system.All three sub-regions are coupled to the ionosphere through field-aligned current loops that transfer angular momentum between the equatorial magnetosphere and the planet's upper atmosphere.At high latitudes, field-aligned potentials and Alfvénic turbulence structures drive field-aligned currents and accelerate electrons, leading to the formation of auroras.
transport and electric currents in this domain (Section 2.4) by means of the magneto-hydrodynamic (MHD) approximation, in which magnetic field lines are assumed to be frozen into plasma flows (see Figure 1).Domains II to IV are coupled along and across magnetic field lines by the electric current loops illustrated in Figure 1, which flow parallel to the magnetic field between the ionosphere and the equatorial magnetosphere, and close across magnetic field lines at the two ends: latitudinally, inside the ionospheric conductor (Domain II), and radially, inside the magnetodisk (Domain IV).We are now going to describe the electrodynamic interactions between domains II to IV, assuming we know the magnetic field generated in Domain I. To this end, we determine in each domain the basic equations linking electric currents to the local properties of each medium, which include electrostatic potential, plasma populations, and their coupling to neutral gas domains (thermosphere and neutral tori).

Domain II: Ionosphere-Thermosphere
The ionosphere-thermosphere domain is modeled as a resistive medium in which electric currents are related to the electric field in the reference frame of the neutral air via a local 3-dimensional Ohm's law.In magnetosphereionosphere-thermosphere (MIT) coupling models, it is usual to note that the vertical extent of this medium is small compared to its horizontal extent, which is on the order of the planet's radius.Thus, Domain II is usually represented as an infinitely thin electric current layer, the ionospheric dynamo layer, connected to the magnetosphere (Domain III) by a system of magnetic-field-aligned currents.Even if this view has been qualified at Jupiter by recent Juno measurements (Mauk & Saur, 2007), it is not the aim of this study to extensively discuss the effect of the ionospheric vertical structure on global magnetospheric transport, and we neglect it in the rest of the paper.
Within this thin ionosphere, the curvature of the magnetic field lines through the ionosphere can be neglected, and the field lines can be approximated as straight lines (see inset in Figure 2).In addition, the electrostatic potential, magnetic field B, and ionospheric distance to the rotation axis are assumed constant over the ionosphere height.Using these approximations, vertical integration of the local 3-D ionospheric Ohm's law leads one to an altitudeintegrated 2-D ionospheric Ohm's law relating the electric field E to altitude-integrated horizontal electric currents J i : is the altitude-integrated 2-D ionospheric conductance tensor, expressed in terms of its Pedersen Σ P and Hall Σ H components, and v is the local plasma velocity.This Ohm's law is written in the rest frame of the neutral atmosphere, which moves at the local velocity v n = v cor + v N of the neutral wind in the planetocentric inertial frame.The local corotation velocity is given by v cor = Ω P × r, where Ω P is the planet's rotation vector, and r the position vector.The velocity of the neutral wind v N in this corotating frame is itself composed of an azimuthal component v N,ϕ and of a meridional component v N,θ , which we assume for simplicity to be constant with altitude over the vertical extent of the dynamo layer.However, there is no reason to assume that neutral winds blowing at two conjugate points in the northern and southern thermospheres are equal.On the contrary, available observations show that wind systems in the two hemispheres, particularly at high latitudes, may be very different.In order to specifically study the large-scale effects of these wind system asymmetries, we will use in the remainder of this study two different neutral wind prescriptions, namely v Nn and v Ns , in the northern and southern hemispheres, respectively.In general, in all the following, quantities that can have different values in the northern and southern hemispheres will be identified with a subscript n or s respectively, or k which can be either n or s.

Domain III: High Latitude Polar and Auroral Field Lines
Our description of the electrodynamics of high latitude field lines threading Domain III is based on the simplified assumption that field-aligned current flows, particle acceleration and precipitation, and Alfvénic turbulence generate a net field-aligned electrostatic potential drop, written Φ // , between the lower and upper altitude limits of each field line crossing Domain III.Since field lines are assumed to be equipotential in the magnetosphere itself (Domain IV, see next sub-section), one can write, that Φ // = Φ i Φ eq is the difference between ionospheric (Φ i ) and equatorial (Φ eq ) electrostatic potentials along each field line.This definition is valid in the frame co-rotating with the magnetic field and the plasma, where the magnetic field is static.
The relationship between field-aligned potential drops and electric current carried by precipitating electrons has been first described by Knight (1973), using a kinetic description of adiabatic electron motions under magnetic mirror forces.In this very useful simplified description, field-aligned currents carried by precipitating electrons correspond to the flux of electrons inside the loss cone at the top of the ionosphere.This flux is itself determined by the electrostatic potential drop Φ // .Using a Maxwellian distribution for the electrons leads to Knight's currentvoltage relationship.Our work uses the modified current-voltage relationship provided by Ray et al. (2009), which takes into account the additional effects of the centrifugal potential and non-monotonic behaviors along the field lines in the expression of field aligned ionospheric currents: where j x , T x , and R xk are the thermal current, the electronic temperature, and the mirror ratio at the location of the gravitational + centrifugal potential minimum.The thermal current writes , with n xk the electronic density at the same location and m e the electron mass.Typical values of those parameters are given by Ray et al. (2010).

Domain IV: Magnetosphere
In the magnetospheric domain, we aim at describing the global transport in the magnetosdisk.Hence, we limit ourselves to the modeling of the cold and thermal plasma.Although the chemistry and charge exchange processes play an important role in the neutral tori around Io at Jupiter, and over an extended radius range from 3 to 9 R S approximately at Saturn, we do not extensively study it, and simply include it as an input parameter in the form of plasma loading, to be determined from specific models of neutral regions.Turbulence is not specifically addressed either, and it simply acts as a heating source for the plasma in our model, as described by, for example, Saur (2004), Liu et al. (2021).The magnetospheric domain is therefore assumed to be a non-relativistic collisionless medium, to which the magneto-hydrodynamics equations apply, written in the planetocentric inertial frame as follows (the reader can refer to Tables A1 and A2 for a summary of all the notations introduced below).
The conservation of mass, where S m,pu is the source of mass due to pick-up of ions, and L m,r is the loss of mass due to recombination: (3) The transport of momentum, including pick-up processes and plasma mass loading through the addition of a Lorentz-like force due to a pick-up current j pu : The induction law for magneto-hydrodynamics: The transport of energy, including a heating rate S q and a loss term L q : Setting a pick-up conductivity σ pu = S m,pu B 2 , and defining E′ the electric fields in the local rest frame of the (nearly) Keplerian neutral cloud, the pick-up current in Equation 4follows an Ohm's law:

Transport Equations
In this section, we aim at deriving transport equations for the mass, angular momentum, and energy contents of magnetic flux tubes in the magnetosphere (Domain IV), based on the magnetospheric equations of conservation of momentum, mass, and energy (Equations 3, 4, 6).
We translate the fundamental vectorial equations 3 to 7 into scalar transport equations, using the Euler potentials, which form the local frame (α, β, ζ) born by the local vectors (e α , e β , e ζ ) (Figure 2).The use of this frame requires the strong hypothesis of an axisymmetrical, fully meridian magnetic field which satisfies ∇ × B ⊥ B. In particular, we assume that the magnetic dipole axis and the rotation axis of the planet are aligned.Although this is a rather valid approximation at Saturn, it is arguable at Jupiter (see e.g.Phipps and Bagenal (2021), Section 5.2, and Appendix B for a discussion on the effect of dipole tilt and breakdown of axisymmetry at Jupiter and Saturn).The local frame is associated with scale factors h α , h β , h ζ , and allows for the use of a curvilinear abscissa s such that ds = h ζ ζ.This special frame and its properties are presented in details in Appendix A.
Limiting our derivation to the inner-and middle-magnetosphere, the transport mechanism in our model must conserve the magnetic flux and the axisymmetry.The interchange instability comes as the right conveyor since it conserves both invariants and it has been observed in the magnetospheres of Jupiter (Bolton et al., 1997;Kivelson et al., 1997) and Saturn (Azari et al., 2018;Rymer et al., 2009), and since most parts of both inner and middle magnetospheres have been shown to be unstable against interchange motions (Ferrière et al., 2001).
In order to derive the transport of integrated quantities from the most general equations, we use a quasi-linear approximation.In this approach, each physical quantity A is averaged at order 0, A 0 = 〈A〉, and the interchange motion acts as a first-order perturbation of null mean value, δA = A A 0 .Based on the axisymmetry of the system, we define averages as azimuthal averages over 2π; they correspond to the equilibrium state of the magnetosphere and the axisymmetrical solutions of our set of equations.The quasi-linear approach has been applied by e.g.Ferrière et al. (2001) to study the stability of plasma against interchange motion in fast rotating magnetospheres and describe flux tube interchange modes.To do so, they focused on zeroth and first order terms only.Based on their work, we derive the global transport of averaged quantities, driven by the combination of first order perturbations into second order terms, by averaging Equations 3, 4, and 6.

Equilibrium State of the Magnetodisk
The transport formalism is derived from the perturbation of a mean static state.In this first section, we focus on the mean static state, in which time variations as well as source and loss terms are neglected.The static state is characterized by the topology of the magnetic field and the distribution of plasma, modeled as a multi-species plasma in a rotation dominated magnetosphere.It is computed, based on the angular momentum transport Equation 4 in Domain IV.
The magnetic field topology can be derived assuming its axisymmetry, from the Grad-Shafranov equation, which is a projection of the equilibrium state of Equation 4 on a meridian plane (Caudal, 1986).
The source of magnetic field was expressed by Caudal (1986) in the case of a single ionic species, and extended by Achilleos et al. (2010) to the case of two ionic species.The latest approach can be generalized to an arbitrary number of ionic plasma species, although this requires that each species be an independent source of magnetic field, and be in equilibrium with its own associated electronic species.The validity of this approximation lowers with more numerous ionic species, and we advice that they be limited to the smallest possible number for practical applications.Outside the inner neutral tori, the magnetodisks were found to be dominated by two groups of ions with similar molecular masses: the "heavy ions," iogenic ions (sulfur, oxygen) at Jupiter and water ions at Saturn, and the "light ions," which are protons in both systems (Huscher et al., 2021;Kim et al., 2020;Nichols et al., 2015;Thomsen et al., 2014).Each group can be rather well approximated by a single fluid description, allowing for the use of the Achilleos et al. ( 2010) model with two ionic species.
The magnetospheric plasma follows a barotropic equilibrium along field lines, described for example by Maurice et al. (1997) from the field-aligned projection of Equation 4. For each plasma species, the field-aligned electrostatic force, due to the field-aligned magnetospheric electric field E 0/ / , balances with the centrifugal and gravitational forces, the potentials of which are added in r , and the pressure gradient.In the previous expression, M P is the mass of the planet, G the gravitational constant, Ω the rotation rate associated with the magnetic field, R and r respectively the cylindrical and spherical radii.Combining the barotropic equilibrium with the perfect gas law and the conservation of individual temperatures along field lines, predicted by Liouville's theorem, the variation of the number density n l0 for each species l (ionic or electronic) along field lines writes as follows, with T l the species temperature, m l its particle mass, Z l its charge number, and k the Boltzmann constant: Combining this equation with the quasi-neutrality of the plasma provides an expression of the magnetospheric parallel electric field as a function of the number densities of all plasma species: A useful approximation of Equations 8 and 9 is that of a single ionic-species plasma.In such a case, the plasma density distribution along a field line is a simple exponential profile: where ionic properties are indexed j, electronic properties e and equatorial values eq.We have neglected the gravitational contribution to the potential Ψ g , which falls below 15% of the centrifugal potential beyond 5.9 Jovian radii, the orbital location of Io at Jupiter, and 20% beyond 4 Kronian radii, the orbital location of Enceladus at Saturn.This expression will be used to derive simplified transport equations of integrated quantities in Section 4.2.3.Unless explicitly stated (cf.Section 4.2.3), this work does not require the specification of the number of plasma species, and we leave this choice to later applications of our model.The Grad-Shafranov equation, Equations 8, and 9 allow for the derivation of the plasma distribution along magnetic field lines and the integration of plasma properties over flux tubes.Departure from the mean static state will result in a transport of those integrated quantities (mass, energy and angular momentum) perpendicular to magnetic field lines.

Transport of Mass
Having set our description of the mean static state of the magnetosphere, we start the study of its perturbations with the transport of the total mass content of flux tubes in the magnetosphere (Domain IV), aiming at the derivation of a diffusion equation for the integrated mass content.However, before getting into the complex formalism of quasi-linear theory and diffusion processes, a first analysis of the MHD conservation of mass (Equation 3) provides a quantity which will be very useful in the following, namely the total mass outflow through a given flux shell Ṁ = ∫ s 2πR〈ρv α 〉ds.In the steady state, integrating Equation 3 over the entire volume inside a flux shell indeed leads to: B respectively the total source and loss of mass in a given flux tube per unit magnetic flux and R P the radius of the planet.Equation 11is not sufficient for a complete description of the transport of mass in our framework.Our model indeed requires the derivation of the equatorial distribution of plasma, for instance to compute the meridian distribution of plasma along magnetic field lines and the disk height.In addition, we aim at taking into account the temporal variability of the magnetosphere to some extent.We thus derive a full diffusion equation, starting from the equation of conservation of flux tubes total mass M = ∫ ρ ds B , with the help of the quasi-linear theory.Noting that in the 2-dimensional space (α, β), the conservation of the magnetic flux translate into a incompressible condition for the 2-dimensional flow ∇ ⋅ v = 0, the 2-dimensional conservation equation for M writes: Splitting all quantities between averaged and perturbed values, and averaging Equation 12 over 2π yields (neglecting terms of order higher than three): The corresponding first order perturbed equation is written, neglecting δS m,pu and δL m : Following Ferrière et al. (2001), we assume periodic temporal variations for the perturbations, with frequency ω.This allows us to write Equation 14 as an ordinate differential equation for δM with respect to β, the resolution of which yields the expression of the mass perturbation: This expression is then injected in Equation 13, resulting in: Journal of Geophysical Research: Space Physics where equatorial parameters are indexed eq.The first two terms of the right-hand side of Equation 16correspond to a diffusive transport.Since averages are performed azimuthally over 2π, the fourth and fifth terms are equal to zero.Neglecting the remaining third term allows us to write a diffusion equation for the mass content of flux tubes: where the full expression of the diffusion coefficient is

Transport of Energy
Under the same assumption used for the derivation of a diffusion equation for the mass in the previous section, and following exactly the same steps, we find a similar diffusion equation for the energy content of flux tubes:

Transport of Angular Momentum
The equation for the transport of momentum in the magnetosphere (Equation 4, Domain IV) can be translated into an averaged infinitesimal transport equation for the angular momentum relative to the rotation axis of the planet: where T L,z = R〈j α 〉B 0 and T pu,z = R〈j pu,α 〉B 0 are respectively the torques per unit volume due to the Lorentz force and the force associated to pick-up.As the pressure gradient force and the gravitational force do not have an azimuthal component, their associated torques do not have a vertical component.The sum of the Lorentz torque on the given flux shell, per unit flux, writes: It is thus directly linked to the flux of electromagnetic current through the flux shell.Therefore, the transport of angular momentum in the magnetosphere can be analyzed through the prism of the conservation of currents and the related theory of magnetosphere-ionosphere coupling.In the following, this latest approach will be used to derive a transport equation for the angular momentum.
In a steady-state, ∇ ⋅ j = 0 is integrated over the flux shell, yielding a relation between the magnetospheric transport of angular momentum through the equatorial currents in Domain IV, and the perpendicular component of the ionospheric integrated currents defined in Section 2.2 in Domain II, J iα,k with k = n or s depending on the northern or southern hemisphere: The left-hand side of Equation 22 is expressed in terms of integrated magnetospheric properties, using the transport Equation 20: The advection term separates into two gradient terms for Ω and R. From there, the relation between the rotation rate of the magnetospheric plasma and the electrostatic potential at the equator ( ∂Φ eq ∂α = R P Ω) allows us to write: In the above, Ṁ⊥ corresponds to the variation of angular momentum of the plasma linked to the transport of magnetic flux tubes to a flux shell with a different rotation rate.ṀR/α is directly linked to the transport of the internal angular momentum of the flux tubes through the disk and has the dimension of a mass changing time rate per unit magnetic flux, integrated over the disk thickness: The modifications of the internal angular momentum of a magnetic flux tube during its outward transport in Equation 26 is due to • the perpendicular deformations of the magnetic field lines during its outward transport (expansion, compression, …), as suggested by the scalar product e R ⋅ e α in the integrand first term, • and to the redistribution of plasma along the magnetic field lines, as suggested by the scalar product e R ⋅ e ζ in the integrand second term, e R being the radial vector of the cylindrical frame.Each of these specific transformations is illustrated in Figure 3.
Finally, the third term on the left-hand side of Equation 24 results from the integral pick-up current, and is computed, applying Ohm's law in the keplerian rotation frame of the neutral gas at the rate Ω K ≃ being the colatitude), with the pick-up conductance defined as Σ pu = ∫σ pu ds (see Equation 7).

Comparison to Previous Results
Equations 17, 19 and 24 generalize diffusion models of mass and energy (e.g., Saur (2004)) and transport models of angular momentum (e.g., Cowley and Bunce (2001)) established in the magnetospheres of gas giants.Under a certain number of assumptions, it is possible to recover the results of the studies mentioned above from our formalism.
We first simplify Equation 17 to recover the diffusion of mass presented for example, by Saur (2004) in the Jovian magnetosphere.The symmetry between the northern and southern ionospheres yields B 0,eq //e θ and allows us to write the variations in α in terms of the equatorial radius of the field line R eq = R P L: In addition, Saur (2004) neglect the time variations of the mass content M 0 as well as source and loss terms.Equation 17 then becomes: Similarly, we can write Equation 19 in terms of a radial diffusion of energy: Using Saur (2004)'s notations for the equatorial plasma density (n), the plasma disk thickness (Z 0 ), the turbulent heating rate Q turb , we write the total mass content of a flux tube M 0 = nZ 0 L 3 , the total energy content corrected by the adiabatic expansion W 0 = 3 2 kTnZ 0 L 3+γ , and the energy source term integrated over a flux tube S q = Q turb Z 0 L γ .Assuming that their normalized radial diffusion coefficient D and D α are linked through D = LD α , we recover their Equations 1 and 2 for the radial transport of mass and energy.
Finally, as detailed in Section 4.2, under simplifying approximations, the left-hand side of Equation 24 describes the radial outflow of angular momentum in the same form as the Hill-Pontius equation (Hill, 1979;Pontius, 1997).Indeed, in Section 4.2.1, we show that, assuming perfectly conducting field lines, the electric potential Φ e derivative is directly linked to the rotation rate of the plasma.In addition, under the hypotheses of hemispheric symmetry and an infinitely thin magnetodisk, the Euler potential α partial derivative translate into radial partial derivatives and the mass outflow terms ṀR/α and Ṁ⊥ write in terms of the mass loading rate Ṁ (Section 4.2.2).Equation 24 then becomes: The underlined terms are identical to the left-hand side of the Hill-Pontius equation used by Cowley et al. (2003) (Equation 15) and Nichols (2011b) (Equation 3).These authors did not take into account the source of plasma due to pick-up, which is equivalent to assuming Σ pu = 0 in our equation.As will be shown in the following, the righthand side of Equation 29is linked to the electrostatic potential of the field-line, hence to the rotation rate of the equatorial plasma in the case of perfectly conducting field lines.One then recovers the right-hand side of Equation 15 of Cowley et al. (2003) and Equation 3of Nichols (2011b).

Magnetosphere-Ionosphere Coupling
As already introduced in Section 2.1, Domains II, III and IV permanently exchange particles, via acceleration of magnetospheric particles and possibly extraction of ionospheric particles along field lines, as well as momentum and energy.Momentum exchange works via the current loops illustrated in Figure 1.In this view, first introduced by Hill (1979) and later developed by Cowley and Bunce (2001) and others, the j × B force associated with radial currents in the magnetodisk tends to accelerate the magnetospheric plasma toward corotation (Equation 20 in Domain IV), while its ionospheric counterpart at the other end of the current loop decelerates the ionospherethermosphere system (Equation 1in Domain II).The torques exerted at the two ends of the current loop exactly balance in steady state.
Having identified a set of equations relating current flows in each domain along this current loop to the electrostatic potential and to local plasma properties (Equations 1, 2 and 24), we can "close the circuit" and calculate all components of the current loop.This is what we do in the following sub-sections.First, we use the condition of non-divergence of the total current system, which holds in the approximation of slow electrodynamics, to derive a closed set of five equations prescribing five unknowns.Its solution determines the exchange of angular momentum between the magnetosdisk and the ionosphere, self-consistently taking into account the influence of field-aligned potential drops.We then introduce a few useful approximations, based on which we present full numerical solutions of our closed set of equations for the cases of the Jovian and Kronian magnetospheres with two different goals.The Jovian model is applied to study the effects of interhemispheric asymmetries on current systems and on the plasma rotation curve, focusing on asymmetries in Pedersen conductances and magnetic field strength and topology.The Kronian application is used as a toy model to investigate the influence of the thickness of the plasma disk on its rotation profile.Finally, having introduced in our formulation an inter-hemispheric asymmetry between the thermospheric neutral winds flowing in the two conjugate hemispheres, we study specifically the effects of this asymmetry on current systems and on the rotation curve of the magnetodisk.

A Closed Set of Equations for the Determination of M-I Coupling
In this section, we exhibit a set of five equations coupling five variables that determine together the state of M-I coupling: • the equatorial electric potential in the planetocentric inertial frame Φ eq , related to the angular rotation rate through R P Ω = ∂Φ eq ∂α ; • the ionospheric electric potentials in the planetocentric inertial frame in both hemispheres Φ in and Φ is ; • the ionospheric field aligned currents at the top of both hemispheres j //in and j //is .First, we use the non-divergence of currents in the thin ionospheric layer (Domain II) to relate horizontal ionospheric currents, described by the ionospheric Ohm's law (Equation 1), to field aligned currents flowing out of the ionosphere and to the ionospheric electric potential, taken in the inertial planetocentric frame, Φ i such that E = ∇Φ i .In an axisymmetric configuration, E β = 0, and the expression of the integrated meridional ionospheric currents along a field line is: In the above, Ω n is the rotation of the neutral atmosphere, I is the angle defined in the inset of Figure 2, and the index i refers to ionospheric values.This equation is valid in both the northern and southern ionospheres, in which k becomes either n or s.
Combining Equation 30 with the angular momentum transport Equation 24 results in a relation between all three electrostatic potentials: Figure 4 illustrates the physical interpretation of each term in the above equation.Its left-hand side is linked to the effective transport of angular momentum through the magnetosphere, while the right-hand side corresponds to the resistive influences of pick-up in the neutral torus and of both ionospheres.The reader can also check Figure 7 of Pontius (1997) and Figure 2 of Brooks et al. (2019) for an illustration of the mechanical analogy for the equation.
To complete the system, we introduce the remaining two unknowns, namely the field-aligned currents in both ionospheres.The null divergence of currents in Domain II, B ∂j ζ / B ∂s = ∇ ⊥ ⋅ j ⊥ , allows us to compute integrated field-aligned currents in each ionosphere, as a function of the integrated perpendicular currents: Equations 30 and 32 allow us to write a first relation between ionospheric field-aligned currents and the corresponding ionospheric electrostatic potentials:  (red), linked to the electrostatic potential Φ eq , results from the competing rotation of two to three resistive regions: two ionospheres (green) and the equatorial pick-up region if the field lines cross the neutral torus (blue).The coupling between those regions is enforced by radial diffusion of plasma through the plasma disk, and the associated field-aligned current loops.
The system is finally closed through two current-voltage relations (Equation 2) linking the field-aligned currents and the field-aligned potential drops in Domain III.
We illustrate the applicability of the general formalism we derived in the following section, addressing simplified cases under common assumption corresponding to the magnetospheres of Jupiter and Saturn.We focus mainly on the transport of angular momentum coupling magnetosphere and ionosphere (Equations 2, 31 and 33) and the transport of mass (Equation 17).

Perfectly Conducting Field Lines
A first, common approximation consists in neglecting the field-aligned potential drop between the ionospheres and the equatorial plane.It allows us to reduce the five-equation M-I coupling system to a single equation for the transport of angular momentum through the magnetosphere.It is the approximation made for instance by Cowley et al. (2003) and Nichols (2011b) to study M-I coupling in Jupiter's magnetosphere and Pontius and Hill (2009) at Saturn.In this paragraph, we show a few results derived from this simplification.
Neglecting field-aligned potential drops, the electrostatic potential is constant along the entire field line: ∂Φ eq ∂α = ∂Φ ik ∂α , Equation 31 becomes a coupling relation between the magnetospheric rotation rate and ionospheric quantities: To clarify the ionospheric influence in the angular momentum transport equation, we simplify notations, introducing two generalized conductances, These terms carry the individual resistive properties of each of the two ionospheres, weighted by their magnetic field strength and their specific geometries.Together with the individual ionospheric rotation rates, they completely characterize the contribution of each ionosphere to the angular momentum transport.Similarly, we can define S pu = R 2 eq B 0,eq Σ pu in the pick-up region.
A straightforward application of Equation 34 is the derivation of the magnetospheric rotation rate without equatorial plasma disk ( Ṁ⊥ = ṀR/α = 0, Σ pu = 0).In such a case, as suggested by Figure 4, the rotation rate in the magnetosphere is driven by the interplay between the rotation of the two conjugate ionospheres, weighted by their properties: This expression is very similar to the plasma rotation rate derived by Brooks et al. (2019) at Saturn, based on a mechanical model taking into account conductances and atmospheric drag asymmetries in steady state.However, we extend their model with the use of more general parameters accounting for the magnetic properties of the two ionospheres.Including all new notations, the angular momentum transport equation then writes: Equation 37 is now formally identical to the Hill-Pontius equation (Hill, 1979;Pontius, 1997) and the Cowley and Bunce (2001) equation, but it generalizes their result to the broader family of magnetospheres with: • interhemispheric asymmetries through two ionospheric parameters: -S = S n +S s 2 , the average of the generalized conductances of the two ionospheres; -the ionospheric rotation forcing Ω*, which is the weighted average of the rotation rates of the northern and southern upper atmospheres, and also the rotation rate magnetic field lines would experience in the absence of any magnetospheric plasma, when the two conjugate ionospheres are directly coupled along field lines; • latitudinally extended magnetodisk through two transport terms for the magnetospheric angular momentum: -Ṁ⊥ , associated to the speeding up of flux tubes; -ṀR/α , associated to the deformation of flux tubes.• potentially extended source regions in the disk through the radial dependence of S pu (e.g., Saur et al. (2004) at Saturn).
For the sake of presenting simple numerical applications of our model, we use the formalism developed here in the following paragraphs.This formalism first allows us to simplify our five-equation system (Equations 2, 31 and 33), into a single equation for the transport of angular momentum (Equation 37).In addition, the above expressions yield a simple parameterization of the effect of ionospheric properties on the rotation profile in the magnetosphere, which is particularly useful for the study of the influence of the interhemispheric asymmetries on the rotation profile of the magnetosphere.

Limit Case of a Thin Disk: Jupiter
Considering the limit of an infinitely thin plasma disk centrifugally confined at the equator allows us to greatly reduce the complexity of the global transport equations.It was used by Saur (2004) and Ng et al. (2018) to derive their mass and energy transport equations, and Hill (1979), Pontius (1997) and subsequent authors to derive their angular momentum transport equations in the magnetosphere of Jupiter.It is particularly well suited at Jupiter, as its magnetosphere is known to exhibit a thin, centrifugally confined equatorial disk of height 1-2 Jovian radii, denoted R J in the following (Liu et al., 2021).
Under this approximation, the vertical structure of the disk and magnetic field lines can be neglected, and all plasma moments can be represented by their integrals assigned to the equator, yielding expressions for the integrated quantities in terms of equatorial properties and of the thickness of the disk H.In particular, the radial velocity becomes v R = v α and the shell-integrated mass outflow Ṁ = 〈2πR eq Hρv α 〉.The mass outflow terms become: We use the formalism developed in Section 4.2.1 in the case of perfectly conducting field lines to present first applications of our model.Injecting the thin-disk expressions of the mass outflow terms in Equation 37 leads to: Note that writing the above equation requires that magnetic flux shells perpendicularly cross the magnetic equator, yielding ∂ ∂α = R J BR eq ∂ ∂R eq .Although it is highly unlikely if interhemispheric asymmetries are taken into account, it is commonly used to analytically model the transport of angular momentum (Hill, 1979;Nichols, 2011b).At the equator, the angle between the magnetic field and the equator is thought to be small, and so is the contribution of the vertical derivative in the total expression of the α derivative.
Two parameters carry the ionospheric influence on the magnetospheric rotation curve: Ω* and S Ṁ.Both parameters vary with the magnetic flux shell, due to the spatial variation of ionospheric conductances (Gérard et al., 2020;Nakamura et al., 2022), ionospheric magnetic field (Connerney et al., 2022), and neutral wind velocity (Smith & Aylward, 2009).However their spatial variations are poorly known (Stallard et al., 2018), and we thus approximate them as constants throughout the magnetosphere.Ω* simply acts as a (constant) scaling factor for the rotation profile; in the following, we normalize all rotation rates to Ω* and focus on finer effects, carried by the parameter S Ṁ.For an insight into the exploration of those specific asymmetries, the reader can refer to Brooks et al. (2019) modeling at Saturn.
For our exploration of the parameter space of interhemispheric asymmetries, we define the variation range of S Ṁ based on observations at Jupiter.Connerney et al. (2022) provided magnetic field amplitude and topology in the northern and southern hemispheres.Based on their Figures 5 and 6, we evaluated azimuthal mean, minimum and maximum values for the northern and southern ionospheric magnetic field intensity at the location of Io footprint, providing orders of magnitude for those parameters a few degrees poleward, in the main oval aurora.In addition, Al Saati et al. (2022) provided a comprehensive statistical study of the ionospheric properties in the northern and southern ionospheres, which the Pedersen conductances.We use the results from their Table 1, limiting ourselves to their "Trend A" corresponding to sub-corotating ionospheres.Mean and extreme values of each ionospheric quantity as well as S in the Jovian case are summarized in Table 1 Ṁ is varied in the range presented by Bagenal and Delamere (2011) for the Jovian plasma source.
Figure 5a shows the rotation curves obtained for various S Ṁ values, and the corresponding locations of the breakdown of corotation (vertical lines), computed as the location at which Ω Ω * falls below 75%, in accordance to the definition of Hill (1979) and following works.We use the magnetic field model developed by Ray et al. (2010), combining the CAN model (Connerney et al., 1981) bellow 21.78R J , and the Khurana and Kivelson (1993) model above 21.78RJ for our simulations.Figure 5b shows the evolution of the corotation breakdown radius R break with S Ṁ.These numerical results can be applied to any Jovian-like system, showing a strong, fast rotating magnetic field, bent by a thin magnetodisk at the dipolar-centrifugal equator.In Figure 5, the special case of Jupiter, computed for the canonical "Mean" value of S presented in Table 1 and another value of S Ṁ deduced from the observed location of the corotation breakdown (Figure 6 and following discussion), is shown amongst the broad family of rotating magnetosphere.
As a step forward toward completeness, we briefly present the resolution of the mass transport equation in a thin disk approximation.In the steady-state and in the absence of plasma source, D α B 0,eq ∂M 0 ∂α is constant through the disk.A dimensional analysis leads us to set: where α 0 = α(R eq0 ) is the inner limit of the integration range.In the limit of a thin disk, M 0 = ρ eq H B 0eq .Converting the α derivative into radial derivative in the above equation and integrating it yields the radial equatorial density profile: R′ eq D α dR′ eq (41) Following Saur (2004) and Ng et al. (2021), we set a power law for the diffusion coefficient: β .The equatorial density then becomes a power-law function of the equatorial radius: Our numerical model for the mass and angular momentum transport in the magnetosphere of Jupiter is compared to data in Figure 6.In Panel 6a, we compare two Jovian models with the Kim et al. ( 2020) radial velocity measurements.The first model assumes canonical values for S and Ṁ taken from the literature (0.7 R 2 J .B J .mho from Table 1 and 1 t.s 1 respectively), and strongly overestimates the data at high radii.The second model uses a value deduced from the corotation breakdown observed in the data: the azimuthal velocity measurements (Figure 6a) indicate that the breakdown of corotation is located at about 35 R J , in agreement with Cowley and Bunce (2001); from Figure 5b, this value corresponds to a S Ṁ ratio of approximately 5 × 10 5 R 2 J .B J .mho.s.kg 1 .Assuming a mass loading at Io of 1-1.5 t.s 1 (Liu et al., 2023), this would indicate a S value of 0.05-0.07R 2 J .B J .mho, much lower than what is shown in Table 1.This is likely due to an overestimation of the ionospheric Ṁ exploring the Jovian range shown in Table 1.The dashed gray curve corresponds to the canonical Jupiter configuration found in the literature with S = 0.7 R 2 J .B J .mho,Ṁ = 1 t.s 1 , while the dashed blue curve is computed for a value of S Ṁ best reproducing the observed location of corotation breakdown at Jupiter (Figure 6   Pedersen conductance by a factor ∼10, which can be explained by the fact that the data taken from Al Saati et al. (2022) are peak values of the conductances, whereas our approach should use mean values.Figure 6b compares our numerical solutions for the density profile at the equator to the Bagenal and Delamere (2011) empirical model below 20R J .The disk scale height is taken from Bagenal and Delamere (2011) and Liu et al. (2021).The parameters β and D α0 are fitted to the empirical model.The solutions satisfactorily fit the Bagenal and Delamere (2011) density measurements for β = 5.6, which is consistent with the value found by Saur ( 2004) (β = 6.6), considering the relation between our diffusion rate and theirs: D α = L 1 D (see Section 3.5).

Limit Case of a Dipolar Field: Saturn
Contrary to Jupiter, the Kronian interhemispheric asymmetries in terms of ionospheric properties, and especially conductances, are not yet well constrained.Hence, we do not focus on this aspect for our application at Saturn, even though it will be included in the S parameter in the angular momentum transport equation.Although the latitudinal extension of the Kronian magnetodisk is thought to be greater than the Jovian one due to the topology of its magnetic field, its effect on transport has not yet been explored.We thus rather choose to study this specific point for the present application.
At Saturn, the magnetic field can be fairly well approximated by a dipole up to 5-10 kronian radii (Achilleos et al., 2010).In the special case of a one-ion-species plasma, this allows us to express the integrated mass outflow terms similarly to the thin disk case, using correcting geometrical factors h d1 , h d2 , h d3 taking into account the field line topology (see Appendix C).The resulting transport equation for the angular momentum at Saturn is: For consistency with the Jovian application in the previous section, we solve the transport of mass for the Kronian extended-disk case.The equatorial density can easily be computed by replacing M 0 in Equation 40 by its expression given in equation C3.Assuming a power-law for the diffusion coefficient D α , the resulting density profile is then very similar to Equation 42, switching from R J to R S (Saturn's radius): Given a plasma temperature profile in the magnetosphere, the transport of angular momentum and mass is solved in the dipolar configuration with the same approach as the thin-disk case.The equations are integrated from 5 R S to limit the influence of the Enceladus torus and any plasma source is neglected in our simulation.
We first present the influence of the plasma disk thickness on the rotation and density curves in Figure 7.It is directly correlated to the plasma temperature (see equation C7): the hotter the plasma, the thicker the disk; a temperature of 0 keV thus corresponds to the thin disk approximation.Livi et al. (2014) showed that temperature profiles for the H + and OH + ions at Jupiter are within the range 1-400 eV.Temperature variations do not greatly  Note.For each physical quantity, we show a mean value and two extreme cases (E1 and E2), from which we compute representative S n , S s and S values, in units of R 2 J .B J .mho = 2.05 × 10 12 m 2 .T.mho, and the range likely to comprise these parameters at Jupiter.impact the modeled rotation and density profiles, hence we set the plasma temperature to be constant throughout the magnetosphere.We fix a value for S Ṁ for all profiles in Figure 7 accordingly to the literature.From Voyager I UVS and radio data, Eviatar et al. (1982) found values of the ionospheric conductance on the order of 0.02 mho, though those values were qualified by Eviatar et al. (1983), who measured conductances of 0.2 mho at Saturn, based on Voyager/PLS measurements.On the other hand, Saur et al. (2004) fitted rotation models to Saturn's magnetosphere and found Pedersen conductances of a few hundredth of a mho for a magnetospheric mass source of 40 kg.s 1 .The mass loading rate in the Kronian magnetosphere, dominated by Enceladus cryo-volcanism, has been shown to be highly variable in time (Bagenal & Delamere, 2011).The value from Saur et al. (2004) is in the lower part of the 12-300 kg.s 1 range presented by Bagenal and Delamere (2011) in their literature review.We set S .B S .mhoand Ṁ = 300 kg.s 1 for our exploration of the influence of the disk thickness.In Figure 7, three solutions are presented for three representative situations based on the range of temperatures measured by Livi et al. (2014) for the equatorial H + and OH + ions: an infinitely thin disk with kT = 0 eV, a strongly inflated disk with kT = 500 eV, and an intermediate case with kT = 30 eV.The thickness of the disk substantially impacts the transport of angular momentum:the plasma stays at corotation at higher radii for a more extended disk (Figure 7a).The impact is smaller on the equatorial density distribution (Figure 7b), and not as straightforward as for the rotation curve.
Secondly, we compare our model to the Livi et al. (2014) equatorial measurements in the Kronian magnetosphere (Figure 8a) to refine constraints on the value of S Ṁ.The influence of the plasma temperature on the rotation and density profiles is taken into account to compute a plausible range of values for S Ṁ, fitting our models to the data in the thinnest and thickest cases shown in Figure 7 (kT = 0 and 500 eV respectively).Figure 8a presents two types of fitting approaches for the rotation curve.The first model (full blue curve) assumes a constant value for S Ṁ throughout the entire resolution domain, and, in doing so, fails to reproduce the slope of the rotation curve at all radii.The second model (full orange curve) uses three values of S Ṁ to fit the measured profile in three designated regions.Those regions are defined, based on the observations of Sittler et al. (2008): the first region, between the Ṁ and β according to the discussion in the main text.We use a mean ion mass of 20 m P , corresponding to water ions.The thickness of the disk is varied through the mean plasma temperature, which values are set as constant throughout the magnetosphere in a range which corresponds to the Kronian magnetospheric environment (Livi et al., 2014).The profiles are normalized to 0.68 Ω S (7a Livi et al. (2014)) and 1 (7b) at 5 R S .
tori of Tethys and Dione, from 5 to 6.5 R S corresponds to a slowly increasing flux tube mass content with radii and rather constant ion production; the second region, between Dione and Rhea, from 6.5 to 10 R S , to slowly decreasing flux tube mass content and ion production; and the third one, outside the orbit of Rhea at 10 R S , to a stronger decrease in mass content and ion production (Figures 14 and 18 of Sittler et al. (2008)).Our second model fits the rotation curve well in all three regions with values of S Ṁ comprised in the ranges [6.5 × 10 5 , 1.1 × 10 4 ] [1.4 × 10 6 , 2.6 × 10 5 ] and [1.4 × 10 4 , 2.2 × 10 4 ] R 2 S .B S .mho.kg.s 1 respectively in the first, second and third regions.Those variations are likely due to either a change in the ionospheric properties at Saturn, especially the ionospheric Pedersen conductances, or to a modification of the mass flux rate Ṁ.The second option would suggest the presence of mass sources and losses in the regions from 5 R S to 15 R S , and any interpretation based on this option would be unreliable, as our model does not include them.Further work is needed to provide conclusive statements in this regard.As for the first option, it is hard to comment on variations of ionospheric properties at such spatial scales, because the mapping between equatorial and ionospheric locations is poorly known.For instance, Palmaerts et al. (2020) mapped Saturn's northern main aurora to an equatorial radii range of 6-17 R S , encompassing almost all our integration region.We finally highlight that solar wind interactions and convective transport start affecting the plasma rotation profile when getting close to 20 R S , the approximate location of the Kronian magnetopause on the day-side.This could explain the departure of the measured rotation profile from the modeled one around 15 R S , as those effects are not taken into account in our axisymmetrical formalism (we refer the reader to Section 5.2 for a discussion on this specific aspect).
Additionally, Panel 8b compares equatorial density profiles obtained through Equation 44with Livi et al. (2014) data for both H + and OH + ions near the equator.For each plasma species, we fix β = 5.6 and we fit our modeled density profile to the measured one for two cases, the thin disk approximation, and a thick disk configuration with  8a) and the mass (8b) to the Saturn data, taken from (Livi et al., 2014).Panel 8a: Fitting the value of S Ṁ to the rotation curve measured by Livi et al. (2014).Both blue and orange spans are delimited by two curves fitted to the data in the thin disk approximation or assuming a constant plasma temperature of 500 eV.The blue models are fitted to the data, holding S Ṁ over the entire radius range [5, 15] R S .To compute the orange curve, three regions are identified, each having a different, constant value for S Ṁ , independently fitted to the data.Panel 8b: Fitting the value of D α0 to the density profiles of H + (gray dots) and OH + (gray squares) taken from Livi et al. (2014), for a given value of β = 5.6.Similarly to Panel 8a, both spans are delimited by two curves fitted to the data in the thin disk approximation or assuming a constant plasma temperature of 500 eV.The thick disk case uses the modeled rotation profile obtained by the 3-region fit (orange span in Figure 8a) for the computation of the geometrical factor h d1 .
a constant kT = 500 eV, hence providing a range of probable profiles (blue filling).Here again, the latitudinal geometry of the disk plays a secondary role in the determination of the mass transport.Our simulations for the equatorial plasma density fail to capture the second-order variation of the data profiles, which can be due to the widely distributed nature of the plasma sources or the complex chemistry taking place in the disk, overlooked in our simplistic model.The last effect would greatly impact the formulation of the transport equation, as it would not allow the simplification of mass outflow terms (Equation C3 to C6).However, the slope is fairly well reproduced for both H + and OH + components beyond 9 R S where the disk is known to exhibit lower plasma sources.
The limitations of our model are clear, as is the direction further development should follow: accounting for the radially extended source regions in the Kronian magnetodisk, and the complex chemistry resulting from the rich neutral supplementation in the disk, as well as better constraining ionospheric parameters and magnetospheric plasma sources from Cassini data.The latest point would also be of major interest for the Jovian model since Juno is now probing the innermost regions of the magnetosphere and the Io neutral torus.Given the high temporal variability of the magnetosphere, for instance in terms of mass loading rate in the Enceladus torus, modeling departure from the steady state could also provide precious insight in the Kronian system.

Summary
This study presents a new derivation of radial transport equations for the conservative quantities, mass, energy and angular momentum, in the corotation-dominated regions of the magnetodisks of gas giant planets, unifying the two previous lines of models which separately described diffusive radial transport of mass and energy on one hand, and angular momentum exchange between magnetosphere and ionosphere/thermosphere on the other hand.Using the assumptions of axisymmetry and steady state or slow temporal evolution of the system, we derived multi-fluid plasma transport equations for these conservative quantities which take into account the finite thickness of the magnetodisk, the turbulence-driven radial transport, the coupling to the two conjugate ionospheres and thermospheres, and the possible existence of field-aligned electrostatic potential drops.Our formalism consists in 10 mutually coupled equations: three setting the equilibrium state of the system, two describing the transport of mass and energy content, and five remaining equations dedicated to the ionospheremagnetosphere coupling processes driving angular moment transport in the magnetosphere.
Inspired by observations carried out in the magnetospheres of Jupiter and Saturn by the Galileo, Cassini and Juno missions, we based our description of transport on the quasi-linear development of flux tube interchange, under which both magnetospheres happen to be unstable everywhere at radial distances larger than those of the main plasma generation regions (Io and Enceladus tori).Our results reasonably well reproduce the observed radial transport of the mass, energy, and angular momentum contents of flux tubes under the assumption of axisymmetry.
We self-consistently included the exchange of angular momentum between the two conjugate thermospheres/ ionospheres of the planet and the magnetodisk via field-aligned electric currents.This exchange tends to maintain magnetospheric plasma toward corotation and correspondingly slow-down thermospheric motions in the inertial frame.The possible generation of electric potential drops along field lines by the effect of field-aligned currents has been modeled, using the Ray et al. (2009) formula.
Unlike previous analytical transport modeling, our formalism is general enough to be applicable to the global transport of all three conservative quantity (mass, energy, and angular momentum) in the two gas giant systems in the Solar System, taking into account the particularities of each of them.It was especially built to include all the new observations of the Juno mission, ignored by most transport formalisms at Jupiter, which were developed before the mission was launched.We applied this general formulation to specific applications to the cases of Jupiter and Saturn.
At Jupiter, recent Juno and telescope observations showed significant inter-hemispheric asymmetries in planetary magnetic field, field-aligned currents, ionospheric conductances and auroral emissions which had not been taken into account in magnetospheric transport.Jovian applications of our model allowed for a first investigation of their influences on the magnetospheric rotation curve.We showed that, neglecting field-aligned potential drops, interhemispheric asymmetry effects could be described by the interplay of two parameters: the average of the generalized conductances of the two conjugate ionospheres given by Equation 35 and the average of northern and southern upper atmospheres' rotation rates weighted by the associated generalized conductances.The values of those parameters are poorly constrained, and even less are their variations with magnetic flux shell.Further observational determination of those parameters, combining both Juno measurements and remote observations, will be required to feed our model.In addition, in order to self-consistently account for interhemispheric asymmetries and the resistive properties of the two conjugate ionospheres, further developments of this formalism will need to include field-aligned currents and potential drops in the system, and their possible effects on ionospheric conductances (e.g., Nichols and Cowley (2004)), taking advantage of new in-situ Juno observations of the northern topside polar ionosphere by the Juno extended mission.
At Saturn, spacecraft observations have pointed to the presence of a vertically extended plasma disk in the magnetosphere (Kellett et al., 2009;Krimigis et al., 2007).However, most transport models at Saturn have approximated the disk as a thin plasma sheet (Pontius & Hill, 2009).The Kronian application of our model made it possible to study the disk thickness influence on the predicted rotation curve and equatorial density profile, using a simple, dipolar field and a single-ion magnetospheric plasma model.A first comparison to observations showed a moderate impact of the plasma disk thickness on the transport of angular momentum and mass.It also emphasized the complex radial structure of the Kronian magnetodisk and the need for radial variations of the ratio between ionospheric properties and mass outflow rate S Ṁ to properly account for the magnetospheric rotation profile.Further work will focus on reducing the strong hypotheses made at Saturn and numerically solving our set of transport equations in a more realistic configuration of the magnetosphere and the magnetodisk.

Application Domain and Perspectives
The main limiting hypothesis of the present formalism is the axisymmetry of the system, and particularly of its magnetic field, about its rotation axis.Abundant new observations at gas giant magnetospheres made it possible to assess the relative importance of actual departures from axisymmetry.To cite the most important ones: a) the importance of large-amplitude longitudinal and latitudinal variations in the magnetic field near the planet, more obvious at Jupiter (Connerney et al., 2018(Connerney et al., , 2022)); b) the importance of solar wind-magnetosphere coupling, which generates observable local time and temporal variations in auroral emissions, plasma flows, and current systems (Kita et al., 2019;Yao et al., 2022) at Jupiter, and is likely to cause a warping of the magnetodisk at Saturn (Arridge et al., 2011).
In the case of Jupiter, although Cowley and Bunce (2001) argued that the tilt of the magnetic field axis should not cause longitudinal variations in auroral properties, they did not take into account equatorial transport processes at large radial distances, where the flapping of the magnetodisk caused by this tilt is more important.In Appendix B, we show that the expected longitudinal variations of the dynamical environment caused by the magnetic field tilt at the magnetic equator are smaller than 2% and probably dominated by other types of asymmetries.For instance, Ray et al. (2015) showed that the thermosphere and the ionospheric conductances play a major role in angular momentum transport, so that asymmetries in zonal flows or conductive properties need to be accounted.At Saturn, Kaminker et al. (2017) showed important local time asymmetries in turbulent heating, which need to be taken into account in the transport of energy.While our formalism does not account for asymmetries in longitude and local time, it is a first step toward their future inclusion in global transport modeling.Indeed, keeping the perturbative approach of this work, and starting from its current description as a mean state, it will be possible to add small amplitude perturbations of the rotational symmetry at large spatial scales to represent azimuthal asymmetries.
In a near future, developments of this work for the Kronian and Jovian magnetospheres, combining simulations with existing observations, will make it possible to consistently take into account the influence of field-aligned potential drops and interhemispheric asymmetries on the transport of angular momentum, mass, and energy in the magnetodisk and to study the turbulent heating of magnetodisks.We will also include the plasma source terms of the two magnetospheres and extend our resolution domain to the Io and Enceladus tori.
Our general formalism is intrinsically applicable to exotic magnetospheric systems, and especially to the wealth of exo-Jupiter hosted by extrasolar planetary systems.In the past recent years, an increasing number of studies have investigated the detectability of Hot Jupiter' radio emissions generated by the interactions between stellar winds and planetary magnetospheres (Ashtari et al., 2022), as well as their internally powered emissions Journal of Geophysical Research: Space Physics 10.1029/2023JA032233 generated by the interaction between their ionospheres and their magnetospheres (Nichols, 2011a).Conversely, Weber et al. (2017) more recently showed that the specific conditions prevailing at the locations of Hot Jupiter could result in inflated ionospheres and dampened radio-wave emissions.The type of accurate physics-based modeling of the plasma environments of Jupiter-like exoplanets presented in this study will be crucial to predict their emissions and help design efficient detection strategies.

Appendix A: Notations
We provide a detailed summary of the notations used throughout the entire paper in Tables A1 and A2.In the following, we give some properties of the local reference frame (e α , e β , e ζ ) we use (see Figure 2).This frame involves both Euler potentials α and β such that Transforming from the spherical coordinates system to the local coordinates system is possible through the relation:   We choose the Eulerian system of coordinates due to its useful properties for the description of the magnetic field.First of all, the magnetic field is borne by the curvilinear vector: ∇ ζ//e ζ //B.Its norm is B = ‖B‖ = 1 h α h β .More importantly, the conservation of the magnetic flux can be simply expressed using Euler potentials.The surface vector on a magnetic flux shell writes dS ⊥ = h β dβdse α = r sin(θ)dϕdse α .And the surface vector parallel to the magnetic field lines is dS // = h α h β dαdβe ζ = dαdβ B e ζ , yielding the following expression for the magnetic flux: dΦ = B ⋅ dS // = dαdβ.The latest directly implies that the magnetic flux is conserved along magnetic field lines.
Vector operators on scalar or vectorial functions (resp.f or F) can be expressed in the local frame as follows (Ferrière et al., 2001) On the whole, it appears that the departure from the axisymmetry due to the magnetic field topology in the magnetospheric domain (Domain IV) shows similar contributions from the dipole component of the field and all other higher-order multipolar components.The resulting longitudinal variations in centrifugal potential, azimuthal velocity, magnetic field strength and magnetic pressure in the dipole equator are on the order of a few percents.

Figure 2 .
Figure 2. Euler coordinates (α, β, ζ) are a useful set of coordinates to describe magnetized media, due to their properties in terms of the conservation of magnetic flux.They are born by the frame vectors (e α , e β , e ζ ) shown here in the planetocentric reference frame.The inset presents notations used for the ionospheric parameters.The ionosphere is assumed thin enough to neglect the curvature of the magnetic field lines through it; under this approximation, we define the inclination I of the field line as the angle between e ζ and e θ .

Figure 3 .
Figure 3. Origins of the angular momentum modifications in a magnetic flux tube as it moves through the magnetosphere, corresponding to different terms of the angular momentum transport Equation 24.Top: modifications of the azimuthal velocity of the magnetic flux tube when it moves to a region of different angular rate.Bottom: variations of the internal angular momentum due to deformations of the magnetic flux tube (stretching along field lines on the left, perpendicular compression/expansion on the right).

Figure 4 .
Figure 4. Model of angular momentum radial transport in the magnetosphere of a giant planet.The rotation of field lines(red), linked to the electrostatic potential Φ eq , results from the competing rotation of two to three resistive regions: two ionospheres (green) and the equatorial pick-up region if the field lines cross the neutral torus (blue).The coupling between those regions is enforced by radial diffusion of plasma through the plasma disk, and the associated field-aligned current loops.

Figure 5 .
Figure 5. Numerical simulation of the influence of interhemispheric asymmetries on the magnetospheric rotation profile, neglecting field-aligned potential drops and disk thickness.Panel 5a: Influence of S Ṁ on the magnetospheric rotation profile (in units of planetary rotation rate Ω J ) in the case of Jupiter.Five profiles are compared for five values of S Ṁ exploring the Jovian range shown in Table 1.The dashed gray curve corresponds to the canonical Jupiter configuration found in the literature with S = 0.7 R 2 J .B J .mho,Ṁ = 1 t.s 1 , while the dashed blue curve is computed for a value of S Ṁ best reproducing the observed location of corotation breakdown at Jupiter (Figure 6 and following discussion).Panel 5b: Evolution of the co-rotational breakdown radius with S Ṁ .The Jovian cases are shown in dashed lines, with colors corresponding to these shown in Panel 5a.
Figure 5. Numerical simulation of the influence of interhemispheric asymmetries on the magnetospheric rotation profile, neglecting field-aligned potential drops and disk thickness.Panel 5a: Influence of S Ṁ on the magnetospheric rotation profile (in units of planetary rotation rate Ω J ) in the case of Jupiter.Five profiles are compared for five values of S Ṁ exploring the Jovian range shown in Table 1.The dashed gray curve corresponds to the canonical Jupiter configuration found in the literature with S = 0.7 R 2 J .B J .mho,Ṁ = 1 t.s 1 , while the dashed blue curve is computed for a value of S Ṁ best reproducing the observed location of corotation breakdown at Jupiter (Figure 6 and following discussion).Panel 5b: Evolution of the co-rotational breakdown radius with S Ṁ .The Jovian cases are shown in dashed lines, with colors corresponding to these shown in Panel 5a.

Figure 6 .
Figure 6.Numerical solutions for the specific Jupiter case, neglecting the field-aligned potential drop and the disk thickness.Panel 6a (adapted from Kim et al. (2020)): Rotation velocity for the asymmetrical model of the Jovian magnetosphere, assuming the Ray et al. (2010) magnetic field model, full corotation at 6 R J , and S Ṁ = 7 × 10 5 (full gray line) or S Ṁ = 4.8 × 10 5 R 2 J .B J .mho (full blue line).Numerical results are compared to measurements from Kim et al. (2020) (black dots).Their full corotation profile (purple dashed line) do note exactly match ours (black dashed line).Panel 6b: Profile of the equatorial density, making the same assumptions on the magnetic field as in the top Panel (dotted and full lines).Our profile is plotted against the Bagenal and Delamere (2011) empirical model of the number density multiplied by an average ion mass of 64 m P (blue line), with associated uncertainty featured by the blue span, and is normalized to the empirical value at 6 R J .

Figure 7 .
Figure 7. Influence of the thickness of the disk on the transport of angular momentum (7a, rotation curve in units of planetary rotation rate Ω S ) and mass (7b).Models are computed, setting the values for SṀ and β according to the discussion in the main text.We use a mean ion mass of 20 m P , corresponding to water ions.The thickness of the disk is varied through the mean plasma temperature, which values are set as constant throughout the magnetosphere in a range which corresponds to the Kronian magnetospheric environment(Livi et al., 2014).The profiles are normalized to 0.68 Ω S (7aLivi et al. (2014)) and 1 (7b) at 5 R S .

Figure 8 .
Figure8.Comparison between our transport models for the angular momentum (8a) and the mass (8b) to the Saturn data, taken from(Livi et al., 2014).Panel 8a: Fitting the value of S Ṁ to the rotation curve measured byLivi et al. (2014).Both blue and orange spans are delimited by two curves fitted to the data in the thin disk approximation or assuming a constant plasma temperature of 500 eV.The blue models are fitted to the data, holding S Ṁ over the entire radius range [5, 15] R S .To compute the orange curve, three regions are identified, each having a different, constant value for S Ṁ , independently fitted to the data.Panel 8b: Fitting the value of D α0 to the density profiles of H + (gray dots) and OH + (gray squares) taken fromLivi et al. (2014), for a given value of β = 5.6.Similarly to Panel 8a, both spans are delimited by two curves fitted to the data in the thin disk approximation or assuming a constant plasma temperature of 500 eV.The thick disk case uses the modeled rotation profile obtained by the 3-region fit (orange span in Figure8a) for the computation of the geometrical factor h d1 .
Particle mass, charge number, temperature* m, Z, T

Figure B1 .
Figure B1.Longitudinal variations of the centrifugal properties of the magnetic equator due to the dipole tilt in the Jovian system.The distance of the field-line apex to the rotation axis of Jupiter is plotted against SIII longitude, in the case of an axisymmetric magnetic field, with its "dipole" axis tilted by 10.31°toward a SIII longitude of 196.61°(Connerney et al., 2022).The vertical dotted lines indicate a 2% fluctuation interval around the mean value for each radial sample.

Figure B2 .
Figure B2.Longitudinal variations of the magnetic equatorial properties in the Jovian system for the combination of the Connerney et al. (2020, 2022) magnetic field models.The magnetic field intensity is plotted in the dipole equatorial plane against longitude, considering a dipole axis tilted by 10.31°toward a SIII longitude of 196.61°.It is sampled at constant Jovicentric distances between 5 R J and 50 R J .The vertical dotted lines indicate a 5% interval range centered on the mean value for each radial sample.

Table A1
General Notations, Applicable to All Domains of the System.Quantities Are Defined in the Planetocentric Inertial Frame

Table A2
Specific Notations in Each Domain.Quantities Are Defined in the Planetocentric Inertial Frame Field-aligned current density at the top of the ionosphere j //i Horizontal current (integrated over the ionosphere thickness) J Pick-up conductance (integrated over the flux tube) Σ pu dr = dre r + rdθe θ + r sin(θ)dϕe ϕ = h α dαe α + h β dβe β + h ζ dζe ζ , (A2) where h α , h β and h ζ are scale factors, h β = R R p and ds = h ζ ζ, with s the curvilinear abscissa.
i Velocity of the neutral atmosphere in the planetocentric frame v n Corotation velocity in the planetocentric frame v cor = Ω P × r :