Modeling the Jovian magnetosphere under an antiparallel interplanetary magnetic field from a global MHD simulation

We present preliminary results of a new global Magnetohydrodynamics (MHD) simulation model of the Jovian magnetosphere. The model incorporates mass loading from Jupiter's satellite Io, the planet's fast corotation, and electrostatic coupling between its magnetosphere and ionosphere (M‐I coupling). The basic configuration of the Jovian magnetosphere including the equatorial plasma flow pattern, the corotation enforcement current system, and the field aligned currents (FACs) in the ionosphere are presented under an antiparallel interplanetary magnetic field (IMF) condition. The simulation model results for equatorial density and pressure profiles are consistent with results from data‐based empirical models. It is also found that there are similarities between the FACs distribution in the ionosphere and the observed aurora features, showing the potential application of the simple ionospheric model to the complicated M‐I coupling. This model will help deepen our understanding of the global dynamics of the Jovian magnetosphere.


Introduction
The interaction between a magnetized planet and the supersonic solar wind plasma results in a cavity containing low-density hot plasma, known as a planetary magnetosphere. However, mechanisms that control the dynamics and configuration of a planetary magnetosphere vary from planet to planet. In the case of Earth, the solar wind driven convection and reconnection process involving open flux tubes ("Dungey-cycle") prevails over the whole magnetosphere. In contrast, the Jovian magnetosphere is controlled primarily by two mechanisms, one is the interaction of the solar wind with the planet's intrinsic magnetic field, and the other is associated with the internal plasma source combined with the planet's fast rotation (Vasyliunas, 1983). As one of Jupiter's Galilean moons, Io is the most geologically active object in the solar system. Io's volcanism accounts for the abundant sulfur dioxide (SO 2 ) in its atmosphere, from which the particles are constantly ejected into the Jovian magnetosphere, forming the neutral cloud around Io's orbit at about 5.9R J . The neutral cloud then experiences ionization processes and charge-exchange collisions, which become the dominant plasma source in the Jovian magnetosphere Krupp et al., 2004).
Jupiter has been visited by nine spacecraft so far (Voyager 1/2, Pioneer 10/11, Ulysses, Galileo, Cassini, New Horizons, and Juno). These missions have provided us much knowledge of the Jovian magnetosphere. Apart from in-situ observations, global MHD (Magnetohydrodynamics) simulation is an important method to investigate the planetary magnetosphere (see a review by Wang C et al. (2013)), which can reproduce many large-scale magnetospheric dynamics. However, to date only a few Jovian MHD simulation models have been developed. In the early days, the effects of the ionosphere and the Io plasma source were not seriously considered. For example, the locations of the inner boundary were set far away from the planetary surface: 30R J (1R J = 71492 km, the radius of Jupiter) in Miyoshi and Kusano (1997), Miyoshi and Kusano (2001), 21R J in Ogino et al. (1998), Walker et al. (2001 and Fukazawa et al. (2005). In these models, the plasma source and the ionosphere were not included. Then in Moriguchi et al. (2008), the ionosphere was introduced through an electrostatic model similar to Hu YQ et al. (2007). However, the inner boundary was located at 8R J and the plasma source Io was not added to the governing equations. Chané et al. (2013) proposed a model to include the mass loading of Io and the ionosphere. They implemented an extended ionosphere region between 4.5 and 8.5R J , where neutral-ion collisions were modeled. However, the Io plasma source was located at 10R J instead of its real location of 5.9R J . The model also enhanced the gravity artificially to 55 times the normal value for some practical and numerical reasons. Due to these unrealistic ionosphere and Io torus assumptions, the inner magnetosphere in their model was not realistic.
In this paper, we will present results of a newly developed Jovian MHD simulation model, which takes into account the effects of the plasma source Io and the ionosphere. The paper is organized as follows. A detailed introduction to our model is described in Section 2. Preliminary simulation results are presented in Section 3. Finally, we conclude this paper with a summary in Section 4.

Basic MHD Equations
We start with MHD equations that include the effect of the Io plasma source. In order to improve the accuracy of the calculation of Jupiter's magnetic field, we use the magnetic decomposition method, in which the magnetic field B is split into two components, the internal dipole field B 0 and the external field B 1 (Tanaka, 1994). Then the equations in a conservative form are expressed as where U is the vector of conservative variables, F(U) represents the corresponding flux, and S(U) expresses the source terms concerning Io's torus and planet's gravity, as follows: where ρ is mass density, and u=(u x , u y , u z ) is velocity, κ is the plasma loading rate from Io's torus, g is the Jovian gravitational acceleration, u n is the velocity of the neutral cloud in the Io torus, and p n (ρ n ) is the pressure (density) of the neutrals in the Io torus. The total pressure p, p 1 , and the total energy density e 1 are defined as where p th is the thermal pressure of the plasma, γ = 5/3 is the adiabatic index. The above equations are implemented in a Godunovtype numerical scheme (Florinski et al., 2013), in which the zoneaveraged MHD conserved values are updated from the numerical fluxes at the grid interfaces at each time step. We use the Riemann solvers to calculate the numerical flux at the interfaces. In our model, the extended HLLC (Harten-Lax-van Leer Contact) and HLLD (Harten-Lax-van Leer Discontinuities) approximate Riemann solvers (Guo XC, 2015;Guo XC et al., 2016) are adapted in order to meet the requirement of the above magnetic decomposition method. The solvers have the advantage of robustness at the low plasma β (the ratio between thermal and magnetic pressures) region, and can keep the same high resolution as their original versions (Li ST, 2005;Miyoshi and Kusano, 2005). To maintain zero divergence of the magnetic field, we use the source term cleaning method to reduce the numerical error caused by the divergence of B (Powell et al., 1999). It should be noted that only one kind of ion, the proton, is considered in our one-fluid MHD model.

Numerical Domain and Boundary Conditions
We take a Cartesian grid system with solar-magnetospheric coordinates, in which the x-axis, y-axis, and z-axis point to the Sun, the dusk, and the north, respectively. Our simulation domain is taken to be -300R J ≤ x ≤ 150R J , -250R J ≤ y, z ≤ 250 R J . In our model, the cells are smallest at the inner boundary of the domain, with a size of 0.5R J , and biggest, 4.1R J , at the outer boundary. The inner boundary is set at 5R J instead of 1R J to avoid the negative effects of the strong magnetic field. At the inner boundary, the radial component of the perturbation magnetic field B 1r is set to zero to satisfy the constraint of zero magnetic field divergence. The radial plasma velocity v r is set to zero so that the mass loading in the inner magnetosphere is entirely determined by the Io plasma source rather than the inner boundary. The transverse velocity of the plasma at the inner boundary is determined by M-I coupling superposed by the planet's rotation, i.e., , where E and B are electric and magnetic field, and Ω is the Jovian angular velocity. For the outer boundaries in the y(z) direction and the downstream of the x direction, open boundary conditions are applied: all quantities are calculated by equivalent extrapolation. On the dayside, the solar wind is prescribed as follows: The initial state in our simulation model is somewhat arbitrarily prescribed. Similar to the approach by Hu YQ et al. (2007), we divide the numerical domain into three sub-regions: the inner sphere of r < 5R J , the region of x < 100R J outside the sphere, and the region of x > 100R J . In the inner sphere, we set the B' = 0, V = 0, n = 4.0 cm -3 , p th = 0.07 nPa. In the region outside the sphere and on the left of x = 100R J , B' is assumed to be produced by the image of Jupiter's dipole located at (x, y, z) = (100, 0, 0)R J , V = 0, n = 20/r cm -3 , p th = 0.35/r nPa. Beyond x = 100R J , a uniform solar wind is prescribed. It usually takes up to 100 hours to reach a quasisteady state.

Magnetosphere-Ionosphere Coupling Method
We use the same M-I coupling approach as Hu YQ et al. (2007) and Moriguchi et al. (2008), i.e., the field aligned currents (FACs) are calculated through Ampere's law at around the inner boundary, the currents are then mapped down to the ionosphere along dipole magnetic field lines, where the electric potential is solved in a two-dimensional (2D) shell for prescribed Pedersen and Hall conductance. For simplicity, the ionospheric height-integrated Pedersen conductance in the present calculation is uniformly set to 0.5 S, and the Hall conductance is twice that value. These values are reasonable compared to values given by Strobel and Atreya (1983), who estimated that the ionospheric conductance should be between 0.2 S and 10 S. The ionospheric electric field is then mapped back to the inner boundary to specify the velocity perpendicular to the magnetic field. For the ionosphere at 1R J , a uniform (θ, φ) mesh is laid out in the domain with △θ = π/180, △φ = π/64.

The Internal Plasma Source Io
In our model, we intend to include the effects of the Io torus; accordingly, the plasma generated by Io is added in the MHD equa-tions by means of a spatially dependent mass density loading rate, κ (kg·m -3 ·s -1 ). Since the plasma loading rate is complicated and unpredicted in reality, we introduce a prescribed distribution function of κ so that the loading rate is strong in the center of the Io torus, decreases smoothly with distance from the center, and is zero outside the torus (Chané et al., 2013), where d is the distance from the center of the torus, and h is the radius of the torus. In this model, the center of the Io torus is at a distance of 7R J in the equatorial plane, and the radius h = 1R J . The plasma loading rate is quite uncertain. Bagenal and Delamere (2011) derived values between 260 and 1400 kg/s for the mass loading rate. Here we use a moderate loading rate of 800 kg/s. However, other plasma sources, including the ionosphere, Jupiter's rings, and other icy satellites, are neglected since their contributions are quite small compared to that of the Io torus .
Note that all the results presented in this paper are obtained when the Jovian magnetosphere reaches a quasi-steady state. The plasma flow around Jupiter is predominantly in the corotational direction. However, this azimuthal flow is also modulated by Vasyliunas cycling. As seen in Figure 1a, the azimuthal velocities around Jupiter are higher on the dawn side than on the dusk side, which is in agreement with observations (Krupp et al., 2001;Woch et al., 2004). This local time asymmetry in plasma azimuthal flow can be explained as follows: during a Vasyliunas cycle, the flux tubes on the dawn side are depleted by tail reconnection and have higher velocities as a result of the conservation of angular momentum, while on the dusk side, the flux tubes will be massloaded and have lower velocities (Palmaerts et al., 2017). The white solid lines show the magnetic neutral lines. It is found that the X-line at the magnetotail appears obliquely distributed at a radial distance between 60R J and 80R J in the regions near Jupiter, which is consistent with in-situ observations . The high-speed tail flow induced by magnetic reconnection can also be seen in Figure 1a.

Numerical Results
In Figure 1b, the magnetic field lines (plotted with white solid lines) are projected onto the noon-midnight meridian plane. As presented in this panel, the magnetic field lines are elongated in the tail and compressed in the dayside. In this case, the IMF is antiparallel to the Jovian magnetic field, which results in an open magnetosphere with a subsolar magnetopause located at 46.05R J and a bow shock at 56.90R J . For runs of a parallel IMF (not shown here), the subsolar positions of the magnetopause and bow shock are 54.45R J and 67.06R J , respectively. Our magnetopause and bow shock are closer to Jupiter than the average values proposed by Joy et al. (2002), who suggest that the range of the subsolar magnetopause positions is between 50R J and 100R J , while the range of the subsolar bow shock positions is between 55R J and 125R J . Since we use protons in our model, whereas the Jovian magnetosphere in reality is dominated by heavy ions (mainly S and O) derived from its plasma source Io , the centrifugal force acting on the rotating plasmas in our model is weaker than that in the real Jovian magnetosphere, which accounts for our model's closer magnetopause and bow shock positions. Figure 2 shows the simulated density and thermal pressure profiles along the tailward Sun-Jupiter line. For comparison, we plot the results of the data-based empirical models in the figure. The solid lines in Figure 2a show the density obtained from the empirical models, while the dashed lines indicate the results multiplied by 22, where we assumed a mass of ~22 amu for the ions similar to Chané et al. (2013). It should be noted that the number density obtained from the empirical models is for the observed ions, whose components are complicated and unknown. Since we consider only protons in our single-fluid MHD model, we cannot compare our results directly to those of these empirical models. Here the Frank et al. (2002) model is based on Galileo data in the magnetotail, and the Bagenal and Delamere (2011) model is based on Voyager 1 and Galileo data. Since we do not know the solar wind conditions and the initial state of the magnetosphere at the time of these observations, some discrepancy between them would be expected. In Figure 2a, the small peak at around 7R J in the density profile indicates the effect of Io plasma loading. It can be seen that the density and pressure profiles we obtain are reasonable, except in the inner magnetosphere (r<15R J ). This disagreement could be attributed to the MHD approximation in the inner magnetosphere, where the non-MHD processes, such as the chargeexchange process and the energetic ions, may need to be taken into account. The charge-exchange process between the energetic ions and neutral cloud that extends from Io inside of 10R J , together with the energetic ions with energy ranging from 20 keV to 50 MeV, which dominate (by about a factor of 10) the pressure in the Jovian plasmasheet Bagenal and Delamere, 2011), may greatly modify the magnetosphere since they do not add mass to the system but affect the momentum and energy transfer in plasma systems (Jia XZ et al., 2012). Figure 3 displays the azimuthal velocity in the equatorial plane along different local times. The black oblique dotted lines display the rigid and 75% of the rigid corotation velocity. Plasma in the inner magnetosphere corotates with Jupiter's rotation, and the angular momentum needed to sustain the corotation is from the M-I coupling. That is to say, in the Jovian atmosphere/ionosphere, plasma acquires momentum from the elastic collisions with neutrals, and the momentum is then transferred to the magnetosphere through FACs. As a result of the conservation of angular momentum, the plasma will lag the planet's rotation on its way moving outward. The corotation then breaks down at a certain distance where the M-I coupling cannot provide enough momentum (Hill, 1979). As shown in Figure 3, the average location of corotation breakdown (as defined by (Hill, 1979), 75% of the rigid corotation) is about 27R J . The azimuthal velocity in the equatorial plane is not uniform but depends on the local time; it is higher on the dawnside (06:00 LT) than on the duskside (18:00 LT). There is a small dip at around 7R J indicating the effect of the mass loading of Io, because the newly injected plasma moves with a lower velocity than the background plasma. Since our magnetosphere in this case is highly compressed, as mentioned earlier, the azimuthal velocity slows down more quickly. The plasma at a radial distance between 9R J and 20R J rotates at a speed slightly higher than the rigid corotation. This is because of the effects of the inner boundary. The plasma in the boundary layer at middle latitudes rotates a little faster than rigid corotation as a result of numerical effects, then the plasma with super-corotation velocity will provide superfluous angular momentum to outer regions. Similar numerical effects of the inner boundary can also be seen in Earth's MHD models, where radial mass flow from the inner boundary is generated even when the radial velocity of the inner boundary is set to zero (Welling and Liemohn, 2014).  Figure 4a, the current is directed inward at the pole and outward at around 15° colatitude. There is a segment of the outward current, located between 7:00 LT and 12:00 LT, where the current is significantly weaker than that of the other local times, forming a discontinuity for the ionospheric FACs. This is consistent with Hubble Space Telescope aurora observations by Radioti et al. (2008), who suggest that there is a discontinuity in the main emission located between 8:00 LT and 13:00 LT. As explained by Chané et al. (2013), this discontinuity is associated with a weaker equatorial radial current, which corresponds to a weaker decrease of the angular velocity with radial distance in the pre-noon sector (Figure 1a and Figure 3). Besides, the FACs are stronger at the nightside and post-noon sectors, which can also be seen in the pattern of Saturn's ionospheric FACs simulated by Jia XZ et al. (2012). Figure 4b displays the corotation enforcement current system on the nightside of Jupiter in the noon-midnight meridian plane. The radial current is shown with color contours, and the lengths of the vectors are defined by the amplitude of the current in the same plane (i.e., J φ is not shown). As seen in Figure 4b, the radially directed corotation enforcement current that flows in the equatorial plane is closed by direct Birkeland currents (outward FACs) and the return Birkeland currents (inward FACs) of both semispheres, which is consistent with the theories of Vasyliunas (1983).

Summary
We present a new global MHD simulation model of the Jovian magnetosphere. This model includes Jupiter's rotation, the plasma loading from its satellite Io, and a simplified planetary ionosphere. We ran a case under a typical solar wind condition with a due northward IMF. The numerical results show the basic configuration of the Jovian magnetosphere, such as its equatorial plasma flow pattern, the corotation enforcement current system, and the FACs in the ionosphere. The density and pressure profiles are consistent with data-based empirical models except in the inner magnetosphere at a radial distance less than 15R J . To avoid this discrepancy, non-MHD effects such as charge-exchange between plasma ions and neutrals and the pressure of energetic particles should be taken into account. The ionospheric FACs are also presented and found to be consistent with main features of the observed aurora, implying that the simple ionosphere model might be used to predict roughly the basic structures of aurora observed in the Jovian ionosphere. The maximum inward current value is 8.71 nA/m 2 ; the maximum outward current value is 4.49 nA/m 2 . The colatitude is plotted with grey dotted circles, and the local time is also shown in the figure. (b) The corotation enforcement current system on the nightside of Jupiter in the noonmidnight meridian plane. The radial current is shown with color contours, and the current in the same plane is represented by the black arrows. The filled green circle represents Jupiter's location and the inner boundary. The vector scaling is given in the upper right corner.