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

1.   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₂) 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 (Khurana et al., 2004; 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.

2.   Model

2.1   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₀ and the external field B₁ (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₁, and the total energy density e₁ 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 Godunov-type numerical scheme (Florinski et al., 2013), in which the zone-averaged 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.

2.2   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., ${{v}} = \left( {{{E}} \times {{B}}} \right)/{B^2} + {{\varOmega}} \times {{r}}$ , 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 < 5 R_J, the region of x < 100 R_J outside the sphere, and the region of x > 100 R_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 quasi-steady state.

2.3   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.

2.4   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 equations 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 (Khurana et al., 2004).

Note that all the results presented in this paper are obtained when the Jovian magnetosphere reaches a quasi-steady state.

3.   Numerical Results

Figure 1 shows the configuration of the Jovian magnetosphere in the equatorial (Figure 1a) and noon-midnight meridian (Figure 1b) planes. Color contours of the plasma density in logarithmic form are displayed in both panels to delineate magnetospheric boundaries such as the magnetopause and the bow shock. In Figure 1a, the equatorial plasma velocities are shown with white arrows, the lengths of which are defined by the magnitudes of the velocities. A global view of the plasma flow pattern in the equatorial plane of the Jovian magnetosphere can be seen in this panel. 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 mass-loaded 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 (Krupp et al., 2004). The high-speed tail flow induced by magnetic reconnection can also be seen in Figure 1a.

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 (Khurana et al., 2004), 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  1.  Magnetospheric configuration under northward IMF condition. The filled white circle at the origin of coordinates represents the inner boundary of our global MHD model. The color contours of density are shown in each panel. (a) Equatorial plane. The plasma velocities are illustrated by white arrows, whose lengths are defined by the amplitudes of the velocities. The vector scaling is given in the upper right corner. The white solid lines show the positions where B_z = 0. (b) Noon-midnight meridian plane. The magnetic field lines are shown in white. The solar wind comes from the right.

DownLoad: Full-Size Img PowerPoint

Figure  2.  The plasma parameters as a function of radial distance. Comparison between our results (shown in black): number density n (a) and thermal pressure p_th (b), the fits to model based on observations by Frank et al. (2002), and Bagenal and Delamere (2011). The solid lines in panel (a) indicate the density obtained from the empirical models, and the dashed lines display the results multiplied by 22 (the assumed average ion mass). Our results are calculated along the tailward Sun-Jupiter line.

DownLoad: Full-Size Img PowerPoint

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 charge-exchange 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 (Mauk et al., 2004; 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.  Azimuthal velocity in the equatorial plane as a function of radial distance. The black oblique dotted lines display the rigid and 75% of the rigid corotation velocity to indicate the corotation breakdown (Hill, 1979). Our results calculated along different local times (00:00, 06:00, 12:00 and 18:00) are coded with different colors, and the local-time averaged velocity is in blue.

DownLoad: Full-Size Img PowerPoint

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  4.  The current system of Jovian magnetosphere. (a) FACs in the northern ionosphere. The FACs are shown with color contours and line contours, where the outward currents (warm colors) are shown with solid lines and the inward currents (cool colors) are shown with dashed lines. The maximum inward current value is 8.71 nA/m²; the maximum outward current value is 4.49 nA/m². 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 noon-midnight 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.

DownLoad: Full-Size Img PowerPoint

Figure 4 shows the current system of the Jovian magnetosphere. The FACs in the northern ionosphere are presented in Figure 4a. Shown in the background are color contours and line contours of the FACs, where the solid contour lines display the outward currents (warm colors), and the dashed contour lines display the inward currents (cool colors). The local time and the colatitude are also shown in the figure. As shown in 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).

4.   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.

Acknowledgments

This work was supported by grants from Chinese Academy of Sciences (QYZDJ-SSW-JSC028, XDA15052500) and NNSFC grants (41731070, 41574159, 41674146), and in part by the Specialized Research Fund for State Key Laboratories of China.