Axially Asymmetric Steady State Model of Jupiter's Magnetosphere‐Ionosphere Coupling

We present an axially asymmetric steady state model of Jupiter's magnetosphere‐ionosphere coupling with variable ionospheric conductivity dependent on the field‐aligned current density. We use Juno and Galileo data to construct a simple model of the equatorial magnetic field, and develop a method for solving the system of partial differential equations describing magnetosphere‐ionosphere coupling. Using this model, we study the behavior of the system with different radial mass transport rates of magnetospheric plasma and the effect of additional field‐aligned currents associated with Jupiter's nightside partial ring current. We compare the model magnetodisc current intensities with those determined directly from magnetic field measurements in various local time sectors, and find that the value of mass transport rate of 2,000 kg s−1 , larger than usually estimated, better accounts for the observed radial currents. We also find that the inclusion of field‐aligned currents associated with Jupiter's partial ring current helps to explain the local time variation of the radial currents, reducing the discrepancy between the model and the observations.

10.1029/2021JA029608 3 of 18 magnetic field and the temperature of the plasma correlates with enhanced brightness of the main auroral oval. They concluded that a transient enhancement can be caused by an increase in the hot magnetospheric plasma pressure and iogenic plasma outflow. Ray et al. (2014) studied local time (LT) asymmetries of M-I coupling at Jupiter using the Vogt et al. (2011) empirical magnetic field model. While the field model used was LT-dependent, they assumed that M-I coupling is locally axisymmetric. They used an effective ionospheric Pedersen conductivity of 0.1 mho, constant in latitude and LT, and a mass outflow rate of 1,000 kg s −1 . Lorch et al. (2020) used magnetic field measurements obtained from all the spacecraft that have visited Jupiter to map the average radial and azimuthal current intensities in the magnetodisc in radial distance and LT. Since the radial currents are a key part of the M-I coupling current system, these observations provide an important opportunity to study LT asymmetries at Jupiter.
In this paper, we develop an axially asymmetric variation of the Hill (1979) model, with variable ionospheric conductivity dependent on the field-aligned current density. In Section 2 we derive the model equations.
In Section 3 we describe the magnetic field model and our approach to solving the differential equations describing M-I coupling. We present the results of the modeling in Section 4 and discuss them in Section 5. Section 6 summarizes our conclusions.

Partial Differential Equation for Angular Velocity
In this section, we derive the partial differential equation for the plasma angular velocity profile, generalizing previous work to the case of axial asymmetry, thus forming a two dimensional extension of the Hill-Pontius equation (Cowley et al., 2002;Hill, 1979;Pontius, 1997). Calculation of the M-I coupling currents follows Cowley and Bunce (2001), with an angular momentum balance equation derived analogously to that given by Cowley et al. (2002), but now not assuming axial symmetry.
A simple way to map magnetically between the equatorial plane and the ionosphere is provided by Euler's potentials. A magnetic field can be expressed in terms of such potentials and as Both and are constant along field lines because the magnetic field vector is perpendicular to their gradients. In cylindrical coordinates, is usually chosen in such a way that its isosurfaces are meridians of constant magnetic longitude. We assume = and = ( ) . The model field is hence wholly poloidal with zero azimuthal component. Function is magnetic flux per unit azimuthal angle and is sometimes called the flux function. With this assumption For purposes of modeling we consider the internal magnetic field of the planet to be dipolar, for which function is given by where is the cylindrical distance from the magnetic axis, is the distance from the center of the planet, and is Jupiter's equatorial magnetic field strength ( = 4.17 × 10 5 nT in the JRM09 internal field model of Connerney et al. (2018)). Near the planetary surface the internal planetary field is dominant, so that assuming the planet is approximately spherical, the ionospheric is where is the perpendicular distance from the magnetic axis in the ionospheric layer. Since is constant along a field line, we can map between the magnetospheric equator and the ionosphere using 4 of 18 From current continuity, the structure of the current system shown in Figure 1, and the assumption of north-south symmetry it follows that on a given flux shell in a given azimuthal sector where is the radial current intensity in the equatorial current disc integrated through its north-south width, and is the ionospheric height-integrated Pedersen current intensity given by Here, we have assumed that the polar planetary field is near-vertical and of strength 2 (instead of the strict dipole formula = 2 cos ). Ω = 1.76 × 10 −4 rad s −1 is Jupiter's angular velocity, is the ionospheric plasma angular velocity, and Σ * is the effective height-integrated ionospheric Pedersen conductivity. The effective conductivity accounts for rotational lagging of the neutral atmosphere relative to rigid corotation due to ion-neutral collisions, and is reduced compared to the true value by an unknown factor 0 < (1 − ) < 1 , taken to be equal 0.5 following Achilleos et al. (1998) and previous related works (e.g., Cowley & Bunce, 2001;Cowley et al., 2002;Nichols & Cowley, 2004). From Equations 4-7 we obtain = 4Σ * (Ω − ).
We assume that the equatorial plasma is concentrated in a thin disc. The angular momentum per unit mass of the equatorial plasma is 2 ( ) , where is angular velocity. The flux of angular momentum is caused by radial transport of the plasma and its rotation around the planet. The change of angular momentum per unit time in the volume between and + and inside the sector centered at azimuthal angle with angular width is where ̇ is the radial mass transport rate per 2 radians of azimuth (full equatorial circle), and ̇ is the azimuthal mass transport rate per unit radial length. The Lorentz force torque per unit volume about Jupiter's center is × ( × ) , where is the position vector, is the current density, and is the magnetic field. If varies only slowly with on the scale of the sheet thickness, and varies only slowly with , then div = 0 guarantees that varies only slowly with z on the scale of the sheet thickness. Then the -component of the torque acting on the plasma inside the volume element considered is If the mass density of the plasma per unit area of the equatorial current sheet is ( Substitution of Equations 10 and 11 into Equation 9 gives We then substitute the width-integrated radial current given by Equation 8 into Equation 12 to obtain the following partial differential equation for the angular velocity Throughout the paper, we will assume ̇ to be constant at all local times (we discuss this assumption in Section 5), which means that (̇) = 0 . With this assumption from continuity of the mass flow div( ) = 0, PENSIONEROV ET AL.
10.1029/2021JA029608 5 of 18 where is plasma bulk velocity, we find We then can simplify Equation 13 to geṫ Physically correct solutions must converge to almost rigid corotation close to the planet, thus the boundary condition is where inner is any inside the region of nearly rigid corotation. In this paper, inner is taken to be 6 (near the orbit of Io).

Modulation of Ionospheric Conductivity by Field-Aligned Currents
The ionospheric Pedersen conductivity is modulated by precipitation of electrons accelerated by fieldaligned voltages in the auroral region. Nichols and Cowley (2004), using the results of Millward (2002), calculated how the height-integrated ionospheric Pedersen conductivity depends on the outward field-aligned currents density, and provided the following analytic approximation where ‖ is the field-aligned current density just above the ionospheric layer. The effective conductivity is then where Σ 0 is the background height-integrated ionospheric Pedersen conductivity taken to be 0.1 mho (Nichols & Cowley, 2004).
If we assume the absence of electric currents perpendicular to the magnetic field lines in the region between the current disc and the ionosphere, then ‖∕ is constant along the field lines. Following Cowley et al. (2002) we then find the equatorial field-aligned current density from the divergence of the radial currents If there exists a partial ring current in Jupiter's magnetosphere, we should add the divergence of the azimuthal magnetodisc currents ∇ to the equation to obtain We will later define these additional field-aligned currents explicitly as an input parameter for the model. The ionospheric field-aligned current density is then 6 of 18 Equations 16 and 22 constitute a system of partial differential equations (PDEs) for and ‖ . Their solution requires a second boundary condition at some distance 0 Ideally, we would like to set 0 = inner and solve system of Equations 16 and 22 radially outward from inner . But, this system is unstable in the near-rigid corotation region and is nearly impossible to solve this way.
In the next section, we will discuss the reasons for this in detail and describe our approach to obtaining approximate solutions.

Magnetic Field Model
The model equatorial magnetospheric magnetic field employed has been derived from Galileo and Juno magnetometer data, using data from all Galileo orbits and from Juno perijoves 0-22 that are currently available. As in Lorch et al. (2020), we split the data into eight 3 h wide LT sectors, yielding sufficient data to cover radial distances of interest, and allowing us to readily compare our results with those of Lorch et al. (2020). For our purposes, we are interested only in the equatorial magnetic field. To determine whether or not a data point is inside the current disc we used two conditions. The first one requires spacecraft to be closer than 4 to the center of the sheet according to the Khurana and Schwarzl (2005) model. We use a large half-width of 4 instead of the commonly used value of 2-2.5 because sheet crossings predicted by Khurana and Schwarzl (2005) are often shifted from the observed ones. The second condition requires the -component of the magnetic field to be smaller than 2 nT. Its purpose is to select the real crossings from the broad intervals picked by the first condition.
In each LT sector, we fitted the polynomial to the component of the residual magnetic field, that is normal to the local current sheet surface. The residual field was obtained by subtracting the JRM09 model of the internal planetary field (Connerney et al., 2018) from the data (We note that for the modeling we still used a dipolar internal magnetic field because of the assumed north-south symmetry and the neglect of dipole tilt). Because local time sector at 9 h is not covered by Juno and Galileo trajectories beyond 40 we were unable to get a valid fit for it. We instead used an average of two neighboring sectors (6 and 12 h) fits. Figure 2 shows the data and the fit for the midnight sector.
Because of the sparsity of data in the inner region and to avoid divergence for 10 , we used the current disc field model developed by Pensionerov et al. (2019) instead of the polynomial given by Equation 24, smoothing the transition between the two field regimes by linear interpolation. Table 1 lists the coefficients of the polynomial used in the LT sectors 00-21 (labeled by their central local time value), while Figure 3 shows the resulting approximation combined with the dipolar magnetic field and its magnetic colatitude mapping. Function for this field model was then obtained by integrating Equation 2 from the small non-zero value of = 0.01 to the local value. Because we use the magnetic field model as an input for the system of Equations 16 and 22, we solve it on a fine grid, but in 3 h wide LT sectors. The outer boundary for our solutions is set to be at 70 .

Method for Obtaining Numerical Solutions
Most of the Hill-Pontius differential equation solutions quickly diverge to very large positive or negative values of angular velocity in the inner magnetosphere. The physically correct solution converges to ∼Ω . In the one-dimensional case with an explicitly defined ionospheric conductivity function Σ * ( ) , the solution can easily be obtained by solving the equation radially outward with the boundary condition ( inner ) = Ω . This is because solutions with slightly different boundary conditions near the planet quickly converge to the one solution we are interested in as shown in the appendix of Cowley and Bunce (2003).
This approach cannot be applied to the Hill-Pontius equation combined with the equation for ionospheric conductivity modulation by field-aligned currents, because it becomes unstable in the inner region. Since we cannot start the solution from the rigid corotation region, we cannot ensure the fulfillment of ( inner ) = Ω by solving the equations radially outward. Nichols and Cowley (2004) deal with this problem by solving the equations radially inward. A solution obtained in this way eventually diverges to a large negative or positive angular velocity, with the physically correct solution sitting in-between. This family of solutions maps to a range of boundary conditions, where larger boundary angular velocities correspond to the solutions diverging upward and the smaller ones to the solutions diverging downward. It allows to use a binary search to find the boundary value that corresponds to rigid corotation near the planet. Nichols and Cowley (2004) fixed the field-aligned current at a distance of 100 and used a binary search to find the physically correct angular velocity. Tracing the solution deep inside the near-rigid corotation region requires the boundary value to be specified to a large degree of precision, quickly exceeding the 64-bit float accuracy. Nichols and Cowley (2004) traced the solution to 10-20 and used an approximate iterative solution in the inner region. The same can be done with a fixed angular velocity and binary-searched field-aligned current boundary condition. In Ray et al. (2010) and Ray et al. (2012), the authors describe two methods for solving a more complex system of equations, that takes into account ionospheric conductivity modulation as well as the effect of field-aligned voltages. The main method they use is the "critical current technique" (CCT) which is faster, but less accurate than the more precise, but far more computationally intensive "constrained predictor-corrector" (CPC), which they use for verification.
The method for solving the one dimensional Hill-Pontius equation with an explicit conductivity can easily be adapted to the two dimensional case. However, the Nichols and Cowley (2004) method for the equation with variable conductivity cannot because the binary search becomes impossible due to the influence of the azimuthal sectors on each other. The crux of the issue lies in the second term of Equation 16 that accounts for the net azimuthal transport of angular momentum 3 . It is useful, therefore, to estimate the significance of this term in comparison with  Table 1 (2011), which together with our magnetic field model yields the number density per unit area, and hence the mass density per unit area assuming an average ion mass of 20 amu.
To estimate the significance of azimuthal transport we used the Hill (1979) analytical angular velocity profile applicable in case of a purely dipolar magnetic field, because it is generally representative of the plasma angular velocity behavior at Jupiter. We also assumed an upper bound for the azimuthal derivative of to be 0.5 × (Ω − ) , consequent on the fact that as the angular velocity converges toward rigid corotation, its azimuthal derivative should converge toward zero. The choice of the coefficient 0.5 is based on estimates of the derivative obtained by solving the one-dimensional system of equations for each LT, separately. Using these assumptions, we compared the electromagnetic torque (the right hand side of Equation 16) with the azimuthal transport term. Figure 4 shows a comparison using Hill's solution with a characteristic distance = 25 (the distance at which the angular velocity starts to deviate significantly from rigid corotation). The azimuthal transport term becomes less significant the closer we get to the planet. This means that in with the coefficients used in the present model for LT sectors 00-21 h (Table 1). We also show the field profile derived by Lorch et al. (2020) 10.1029/2021JA029608 9 of 18 the region with less than some boundary we can neglect the second term in Equation 16 and solve the system of Equations 16 and 22 for each LT sector using the Nichols and Cowley (2004) method. Ideally we would like to be as small as possible to better justify the neglect of the azimuthal transport term. But we need to be large enough for the boundary conditions search to be possible. We picked = 25 as it was the smallest that worked reliably for the search. The result varies with the plasma density and , but in most cases the azimuthal angular momentum transport effect is at least several times less than the electromagnetic torque effect in the region 25 .
From the solutions are traced inwards to 15 . Tracing to the region closer to the planet becomes too computationally intensive. To obtain the solution in the region from 6 to 15 we interpolate the fieldaligned currents at the 15 boundary to zero at 10 using a cubic spline. We then calculate the corresponding ionospheric conductivity and use it explicitly to solve the simple Hill-Pontius equation radially outward in the 6-15 region. Because the field-aligned currents at 15 are comparatively small, the combined solution is practically continuous.
Finally, we use the values of and ‖ at which for each LT sector correspond to a solution that converges to approximately Ω in the inner region as boundary conditions to solve the system of Equations 16 and 22 at . We solve it radially outward on a fine grid in the same 3 h wide LT sectors. We no longer neglect the net azimuthal transport term, using the full left hand side of Equation 16. This approach usually yields a valid solution, though, depending on the magnetic field profile and the chosen fixed boundary condition, it can produce angular velocities and field-aligned currents diverging to large positive values.

Boundary Conditions
As mentioned above, the Nichols and Cowley (2004) method for selecting boundary conditions requires one of them to be fixed. The specific value for it is dictated by the observations. At the moment we do not have detailed time-averaged measurements of the angular velocity at = 25 . Recently, Lorch et al. (2020) derived the divergence of the observed magnetodisc currents, which we could use for our boundary conditions. However, it is calculated as a finite difference of the observed currents, and while the statistical errors of the currents themselves are relatively small, for the divergence they become significant. Instead of the divergence, we opted to use Lorch et al. (2020) radial current measurments directly. Using Equations 8 and 19 we can constrain the boundary conditions by the observed radial current The binary search works slightly differently than in the case when one of the values is fixed. We conducted the search for ‖ , while the corresponding was calculated from Equation 25 on each iteration. Not every value of the radial current has a solution converging to near-rigid corotation. We found that one of the important factors determining whether or not such a solution exists for a given boundary radial current is the value of ̇ . For example, with ̇= 1000 kg s −1 the observed radial currents are typically too large and the equation has no physically correct solutions. In such cases, we iteratively decreased the boundary radial currents from the observed value, each time conducting a new search, until a valid solution became possible. Figure 5 shows a diagram, outlining the algorithm for obtaining the solutions.

Solutions Without Field-Aligned Currents From the Partial Ring Current
We now examine the solutions obtained using the method described in Section 3, and compare the model magnetodisc radial current intensities with those determined by Lorch et al. (2020). The key model parameter is the radial mass transport rate ̇ . For simplicity we assume it to be symmetrical in LT and constant with (we discuss these assumptions in Section 5). Here we compare solutions for the canonical ̇ value of 1,000 kg s −1 , and an increased value of 2,000 kg s −1 .
In Figure 6, we compare the width-integrated radial current intensities calculated from the model results for ̇= 1000 kg s −1 (blue lines) and 2,000 kg s −1 (orange) with the observed current intensities from Lorch et al. (2020) (black) in the LT sectors from 00 to 21 h. For each sector the root-mean-square deviation (RMSD) of the model currents from the observed ones is given. The major observation here is that the 1,000 kg s −1 model systematically underestimates the observed currents within ∼ 40 in all local times except noon. The 2,000 kg s −1 model, on the other hand, comes closer to the observed values in this region, especially in sectors 00/24, 03, 15, and 21 h. However it still underestimates the currents in the 06, 09, and 18 h sectors. The noon sector likely has a radial mass transport rate lower than 1,000 kg s −1 , as both models overestimate the currents there.
In the region closer than ∼15 in some sectors the observed radial currents become negative, most prominently at 21, 00/24, and 03 h. Radial currents in our model always converge to zero near the planet and so it cannot explain these observations. Beyond 40 , the behavior of the observed currents changes depending on LT. In the midnight-dawn sector, the radial current intensities tend to decrease with radial distance very slowly, while in the dusk sector and at 09 h they decrease significantly faster. The model currents are generally more similar in behavior to the observed currents in the dawn sector than in the dusk sector. In the LT sectors 00-09 h the 2,000 kg s −1 model has 40%-50% smaller RMSD than the 1,000 kg s −1 model. In sectors 12-21 h, the situation is opposite, though less pronounced: the 1,000 kg s −1 model has 10%-25% smaller RMSD than the 2,000 kg s −1 model.  Figure 7 demonstrates the angular velocities, ionospheric field-aligned currents, and effective conductivity for the 2,000 kg s −1 and 1,000 kg s −1 cases. Model field-aligned currents in LT sectors 15, 18, and 21 h rapidly increase beyond ∼50 , while the angular velocity falls very slowly or even trends back toward rigid corotation. This is not supported by observations and is likely an artifact of the model, caused by the assumption of constant radial mass outflow and the absence of additional field-aligned currents. These features disappear when we incorporate field-aligned currents from the partial ring current into the calculations in the next section. At some local times, there is a noticeable derivative discontinuity at in all of the variables. This is the result of neglecting azimuthal transport for . These discontinuities are not severe, which indicates that our approximation was justified.
Angular velocities in the case of 2,000 kg s −1 start to deviate significantly from rigid corotation slightly closer to the planet than for the case of 1,000 kg s −1 , thus producing stronger radial currents. The ionospheric fieldaligned currents and conductivity behave similarly for both radial mass transport values. Field-aligned currents are typically in the range 0.1-0.2 A m −2 , with peaks in the dawn sector at 20-30 reaching 0.4-0.5 A m −2 . The corresponding effective conductivities range from 0.1 to 1.3 mho.

Solutions With Field-Aligned Currents From the Partial Ring Current
While Lorch et al. (2020) provide the divergence of the observed azimuthal equatorial current, this divergence, as discussed above, comes with significant statistical errors. The solutions of the Hill-Pontius equation with variable conductivity are very sensitive not only to the magnitude of the ionospheric field-aligned 12 of 18 currents, associated with the divergence of the azimuthal currents ( (∇ ) for brevity), but also to their radial derivatives. The variation of the observed divergence from one bin to the next caused by the said errors, usually renders the equations unsolvable. Thus, instead of using the observed divergences to calculate (∇ ) directly, we employed a simple parametric equation. We found that the radial currents tend to increase when the radial derivative of (∇ ) is positive and to decrease when its negative. This is, to a certain degree, expected because the ionospheric conductivity increases with the density of the field-aligned currents, and the radial currents are proportional to the conductivity. According to Lorch et al. (2020), the azimuthal current is removed from the magnetodisc in the dawn sector and added back to it in the dusk sector. Removal and addition correspond to negative and positive (∇ ) , respectively. If we now assume that (∇ ) decreases in magnitude with distance from the planet, the sign of its derivative becomes consistent with the observed behavior of the radial currents at dawn and dusk. In the dusk sector the derivative of (∇ ) becomes negative, which can explain the faster decrease of the observed radial currents there. The opposite is true at dawn, where the positive derivative slows the decrease of radial currents with distance. On these grounds, we use the following equation for (∇ ) Parameter was set to be 5 , while the chosen values of and for each LT are presented in Table 2. This equation smoothly interpolates to zero at distances smaller than , with the smoothness of interpolation controlled by (as we assume the azimuthal currents become symmetrical in the inner magnetosphere) and its absolute value falls as 1∕ at greater distances. The specific values of and were chosen to best fit the observed radial currents with the model ones. The form of the approximation as well as its parameters are somewhat arbitrary, so we only aim to qualitatively study the effects of (∇ ) . Figure 8 shows an example of (∇ ) at 15 h LT calculated using Equation 26 (solid line), together with (∇ ) calculated from observations (dashed line). Figure 9 shows the radial currents in the same format as Figure 6, but with (∇ ) included. In sectors 12 and 00/24 h LT the addition of (∇ ) did not improve the fit, so we did not include (∇ ) in the final calculations. For the midnight and noon sectors such an assumption seems plausible because for a nightside partial ring current we expect (∇ ) to be present mostly at dawn and dusk. (∇ ) also did not improve the fit at 09 h LT. In this sector the observed (∇ ) is strongly negative, so if its magnitude decreases with distance its derivative is positive. As stated above, this leads to the model radial current falling slower than for the case without (∇ ) , while the observed current at 09 h, on the contrary, decreases very sharply. This behavior might be a result of a change in the radial mass transport rate with distance, unaccounted for in the present model. It also should be noted that, as indicated above, the 09 h sector has substantially less spacecraft coverage, which means that our magnetic field model and the Lorch et al. (2020) results might be inaccurate. Because of this we set (∇ ) to zero in this sector as well. In the rest of the sectors, (∇ ) significantly decreases the RMSD of the model with both radial mass transport values. For ̇= 2000 kg s −1 the RMSD is 40%-60% lower. For ̇= 1000 kg s −1 , the change is less dramatic with RMSD 10%-40% lower than in the case without (∇ ) . Although the addition of (∇ ) generally increases the model radial currents in the 1,000 kg s −1 case, this value still significantly underestimates the observations within ∼ 40 . For 2,000 kg s −1 the improvement in RMSD mostly comes from the region outside 40 . With

Table 2 Parameters and of the Approximate Form for the Ionospheric Field-Aligned Currents Associated With the Partial Ring Current Given by Equation 26, as Employed in LT Sectors 00-21 h for the Models With
Radial Mass Outflow Rates of 2,000 kg s −1 and 1,000 kg s −1 Figure 8. Field-aligned current from the partial ring current versus distance from the magnetic axis for 15 h LT. The solid line shows the parametric approximation given by Equation 26 with = 0.08 A m −2 and = 30 , corresponding to the 2,000 kg s −1 case in Table 2, while the dashed line shows the field-aligned current calculated from the observed currents from Lorch et al. (2020). Inside 20 the observed azimuthal current divergence was not taken into account due to relatively large errors, and was linearly interpolated to zero from the boundary value. Beyond 20 the observed divergence was interpolated by quadratic splines.
14 of 18 the RMSD twofold), the model still does not fit the observations beyond 40 well. Overall, the 2,000 kg s −1 model has significantly lower RMSD than the 1,000 kg s −1 model at most local times. Figure 10 demonstrates the angular velocities, ionospheric field-aligned currents, and effective conductivity in the same format as Figure 7, but for the model including (∇ ) . The anomalies in the dusk sector with diverging angular velocities are no longer present. Although the results changed for the individual sectors, the ranges and the behavior of all variables remain generally the same, with the dawn sector still having larger field-aligned currents peaks at 20-30 than the dusk sector.

Discussion
A key feature of steady state M-I coupling models is their relative simplicity, compared to full 2D or 3D MHD modeling. This simplicity allows one to test the response of the system to different values of various parameters with low iteration time and computing power requirements, and to compare the results with observations. Here, we have developed a variation of this model, which is asymmetrical in LT. This allows us to compare the radial equatorial current intensities calculated using the model with those determined from magnetic field measurements by Lorch et al. (2020) in eight 3 h wide LT sectors centered on 00-21 h.
In this work, we compared the equatorial radial currents produced by the model using radial mass outflow rates of 1,000 and 2,000 kg s −1 , and found that the model currents are in significantly better agreement with observations when a transport rate of 2,000 kg s −1 is used. Currents produced in the 1,000 kg s −1 case are systematically lower than those observed. This result is unexpected, with most estimates of radial mass outflow rate being lower. In Jupiter's magnetosphere, the transport rate is generally assumed to be equal to the Io plasma production rate. Various empirical estimates of the plasma production rate have been made, ranging from 150 to 2,000 kg s −1 (Bagenal, 1997;Bagenal & Delamere, 2011;Broadfoot et al., 1981;Figure 9. Same as Figure 6, but for the models including field-aligned currents from the partial ring current.

10.1029/2021JA029608
15 of 18 Vasyliunas, 1983), and the canonical value of 1,000 kg s −1 has been used in many previous related works (Cowley & Bunce, 2001;Cowley et al., 2002;Nichols, 2011;Nichols et al., 2015;Ray et al., 2014). Nichols et al. (2020) used the canonical value as representative of the typical outflow rate, while using 2,350 kg s −1 for an enhanced plasma production case. Hill (2001) used the value of 2,000 kg s −1 in his calculations, while Nichols and Cowley (2003) studied solutions of the Hill-Pontius equation for various outflow rates from 100 to 10,000 kg s −1 .
The magnetic field observations used by Lorch et al. (2020) to calculate the currents, and those employed by us to construct our model of the magnetospheric equatorial field, are taken from observations made on the trajectories of several spacecraft, thus representing a time-averaged picture. We then take the mass transport rates used in our modeling to correspond to the average mass transport rate in the system and hence cannot explain the 2,000 kg s −1 value as temporarily enhanced. Another potential explanation comes from our neglect of the changes of mass outflow rate with distance from the planet. This can change the results in many different ways, as both the radial derivative of the outflow and the extra azimuthal flow created are a part of the differential equations. Our model fits the data better in the nightside magnetosphere, where the approximation of constant outflow is probably closer to reality. The inclusion of a variable outflow rate is one of the prime avenues for further research. However, the model with the canonical outflow rate of 1,000 kg s −1 underestimates the observations even within 20 , where the variability of radial outflow is probably much less pronounced, than in the outer regions. Thus, we find the outflow variability with distance unlikely to remove the discrepancy.
We also neglected in our calculations. It does not directly affect the angular momentum balance, but can "bend" the azimuthal sectors, changing the effective magnetic field profile. Because the model with 1,000 kg s −1 underestimates the observed currents across almost all of the local times, we find this change unlikely to affect our conclusions and leave it for future work.  Figure 7, but for the models including field-aligned currents from the partial ring current.
16 of 18 Nichols and Cowley (2004) compared their calculated radial currents with those derived from Galileo azimuthal magnetic field data obtained in the midnight LT sector, and found that the model currents better fit the observed values with ̇ set to a larger value of 2,000 or 3,000 kg s −1 . We used the same approximation for the conductivity dependence on the field-aligned ionospheric current density as Nichols and Cowley (2004), as well as the same atmosphere slippage coefficient of 0.5, which affects the effective conductivity. Preliminary tests with a lower slippage coefficient and hence higher effective conductivity did not show an increase in the model radial currents, while tests with a higher coefficient and consequent lower effective conductivity showed a decrease in currents. From this it follows that changes in the slippage coefficient do not alter the systematic underestimation of the observed currents by the 1,000 kg s −1 model. However, more rigorous study is needed on the behavior of the solutions with different approximations for the conductivity dependence on the field-aligned current.
We also considered the effect of field-aligned currents from the partial ring current on the solutions. Because they are sensitive to the radial derivatives of (∇ ) , we were unable to use the observed divergences directly. Instead, we used a simple analytic form for the resulting ionospheric field-aligned currents with parameters individually selected for each LT sector. The direction of the field-aligned currents we used is in agreement with the Lorch et al. (2020) observations. The inclusion of these currents in the model allowed us to significantly improve the agreement between the observed and model equatorial radial currents, reducing the RMSD by 40%-60% in most local times. However, both the form of the approximation and the specific parameters are to an extent arbitrary, with only the sign of the currents being directly tied to the observations. So, while the inclusion of such currents can help to explain the variation in behavior of the radial currents between the dawn and dusk sector, it does not necessarily do so. At local times 09, 15, and 18 h the inclusion of extra field-aligned currents was not sufficient to explain the sharp decrease in observed radial currents. As stated above, a possible way to improve the model in this regard would be to account for mass transport rate variations with the distance from the planet. Bonfond et al. (2015) used HST images to estimate the brightness asymmetry of the main oval of Jupiter's UV aurora. They found that in the southern hemisphere the dusk sector emission is on average ∼3 times brighter than in the dawn sector, while in the northern hemisphere the dusk sector is only ∼1.1 times brighter (possibly due to the northern magnetic anomaly complicating the analysis). As a possible explanation of this asymmetry, Bonfond et al. (2015) suggested the presence of a partial ring current in the nightside magnetosphere, whose field-aligned currents would strengthen the main oval aurora at dusk while weakening it at dawn. The calculations by Ray et al. (2014) are inconsistent with the Bonfond et al. (2015) results, since they predict stronger field-aligned currents in the dawn sector. Our calculations also show stronger field-aligned currents in the dawn sector. However, the total upward field-aligned current is not larger in the dawn sector, because in the dusk sector it covers a significantly wider latitude range. Figure 11 shows the upward field-aligned current integrated over the ionosphere in one hemisphere for each of the 3 h wide LT sectors. This figure shows our results for cases with and without (∇ ) as well as the Lorch et al. (2020) observations with and without azimuthal current divergence. We integrated over the ionospheric colatitudes that correspond to the equatorial radial distances range 6 < < 55 , to avoid the diverging model currents in the case without (∇ ) . The model current is 20%-50% smaller than the observed current at all local times other than noon. For the model currents in both cases and for the full observed current there is no strong dawn-dusk asymmetry. From these results, it follows that the additional field-aligned currents from the partial ring current do not necessarily affect the aurora in the simple way suggested by Bonfond et al. (2015). Additional field-aligned currents change the conductivity of the ionosphere, which in turn changes the angular velocity profile, and hence the field-aligned currents from the divergence of the radial currents. In the LT sectors 15, 18, and 21 h, which have positive (∇ ) , the total positive field-aligned current is less than in the case without (∇ ) , while in the 03 and 06 h LT sectors, which have negative (∇ ) , the total is larger. The resulting field-aligned currents depend strongly on the magnitude and radial derivative of the additional currents.

Conclusions
We have presented an axially asymmetrical variant of the steady state M-I coupling model for the Jovian magnetosphere. We have compared the radial magnetodisc currents calculated using this model with those derived by Lorch et al. (2020) from in situ magnetic field observations. 1. We found that the observed radial current magnitudes require an average radial mass transport rate of 2,000 kg s −1 , significantly higher than the value typically used of 1,000 kg s −1 2. We considered the effect of field-aligned currents associated with the nightside partial ring current on the M-I coupling system. We found that their inclusion allows a partial explanation of the diurnal variations of the magnetodisc radial currents, reducing the discrepancy between the model and observations by 10%-60%, depending on the local time A notable simplification in the present model is the assumption of a constant radial mass transport rate with the distance from the planet. Accounting for changes in the mass transport rate with distance is one of the directions for future work.

Data Availability Statement
Magnetometer data used in this study was retrieved from the Planetary Data System database (https:// pds-ppi.igpp.ucla.edu/) using (