Three-dimensional MHD simulations of Ganymede’s magnetosphere

1. Introduction

[2] Ganymede, the largest moon of Jupiter, has a radius (R_G) of 2634 km, which is slightly larger than the planet Mercury (1 Mercury radius = 2439 km). Ganymede has a special place among moons because it is the only one known to have its own internal magnetic field [Kivelson et al., 1996, 1997]. Its permanent dipole moment has an equatorial strength of 719 nT that is ∼7 times larger than the ambient field and the dipole axis is tilted by 176° from the spin axis [Kivelson et al., 2002]. The internal field is strong enough to form a mini magnetosphere that shields the moon’s surface from the corotating plasma and the magnetic field of Jupiter. This makes Ganymede’s magnetosphere, a small magnetosphere embedded within another larger magnetosphere, unique in the solar system.

[3] The orbits of Galilean satellites lie very close to the rotation equator of Jupiter but they are tilted with respect to Jupiter’s plasma sheet because Jupiter’s magnetic moment tilts 10° from its spin axis. Therefore the external plasma and magnetic environment around the orbit of Ganymede varies with the nearly 10.5 h synodic period of Jupiter’s rotation depending on the moon’s location relative to the Jovian plasma sheet. At Ganymede’s orbit ∼15 R_J (R_J, radius of Jupiter = 71,400 km), the corotating plasma overtakes the moon from the moon’s trailing edge because its speed is greater than Ganymede’s Keplerian speed. The typical flow velocity of the ambient plasma relative to Ganymede is less than the magnetosonic speed in the flow. In such circumstances, a bow shock, which is a common feature of the planetary magnetospheres, is not needed to decelerate and deflect the background corotating plasma.

[4] The Galileo spacecraft flew close to Ganymede six times during the orbital mission. These six encounters are called G1, G2, G7, G8, G28, and G29, the numbers referring the orbits on which close passes of Ganymede (G) occurred. Table 1 summarizes some information regarding the geometry of these flybys. The six passes at significantly different locations both relative to the surface of the moon and relative to the magnetosphere provide us with a comprehensive sampling of Ganymede’s magnetospheric system. Let us briefly review the observations acquired by the Galileo spacecraft during its close encounters with Ganymede. During the first two flybys (G1 and G2), the measured magnetic field increased dramatically near Ganymede, which led investigators to infer the existence of an intrinsic magnetic field and to identify a magnetosphere associated with the moon [Kivelson et al., 1996]. Sharp rotations in the magnetic field were observed across the separatrices between the external and connected field lines indicating strong localized currents. Kivelson et al. [1996] first used a vacuum-superposition model to interpret the observed magnetic signatures in which they subtracted the background Jovian field inferred by Khurana [1997] from the measured field, then fitted the residual with a dipole field. In their continuing studies, Kivelson et al. [2002] used data acquired by the magnetometer on multiple passes to further characterize Ganymede’s internal moments by looking for a solution consistent with measurements on multiple passes and also allowing for a time variable inductive component. The discrepancies between the data and modeled field, such as the absence of sharp field rotations near the magnetopause in the model, were attributed to various effects which are neglected in their vacuum models, such as the currents associated with the interactions between the Jovian plasma and Ganymede.

[5] Changes in energetic ion anisotropy signatures detected by the Energetic Particle Detector (EPD) onboard Galileo showed that Ganymede’s magnetosphere significantly slows the ambient Jovian plasma in regions extending to several Ganymede radii away from Ganymede along the Jovian field lines connected to the moon [Williams et al., 1998]. The measured particle pitch angle distributions were used to estimate the surface field strength and to infer the field geometries, and the results are consistent with the magnetometer observations [Williams et al., 1997a, 1997b]. A variety of wave modes (whistler mode emissions, upper hybrid resonances, electrostatic electron cyclotron waves, radio emissions) like those found within other planetary magnetospheres were detected during the close encounters [Gurnett et al., 1996; Kurth et al., 1997] which also strongly suggested that Ganymede has an extended magnetosphere of its own. In particular, the upper hybrid resonance frequency observed within the magnetosphere was used to infer the electron density profile, which has a peak value of about 400 cm⁻³ above the surface (i. e., in the ionosphere) with a scale height of about 600 km (≈0.23 R_G) [Gurnett et al., 1996; Eviatar et al., 2001]. This information is used in our MHD simulations to define the inner boundary conditions that will be discussed later. Outflowing ions from Ganymede’s high latitude ionosphere detected by the Plasma Analyzers (PLS), indicate the existence of a polar wind similar to that observed in the terrestrial magnetosphere [Frank et al., 1997; Vasyliûnas and Eviatar, 2000]. In addition to the in situ observations, Hubble Space Telescope observed OI airglow features at Ganymede [Hall et al., 1998; Feldman et al., 2000] associated with a Ganymede aurora.

[6] All the observations mentioned above indicate that Ganymede’s magnetosphere is a complex system. To better understand the interaction of Ganymede with Jupiter’s flowing plasma we have used global magnetohydrodynamic (MHD) simulations with high spatial resolution to model the interaction. Global MHD simulations have been widely used for decades and been proven to be very useful in studying the interaction of the solar wind with planetary magnetospheres. Several numerical simulations of Ganymede’s magnetosphere have also been carried out. Kopp and Ip [2002] and Ip and Kopp [2002] used a resistive MHD simulation to study the variation of Ganymede’s magnetospheric configuration for different background field orientations. However, they compared the measured magnetic field with the simulation results for only one flyby (G2). Paty and Winglee [2004, 2006] used multifluid MHD simulations to describe the heating and interaction of different ion species with the magnetosphere. When multifluid treatments are included in their simulation, the magnetospheric size is consistent with that inferred from magnetic field observations. We will demonstrate in this paper that the magnetic field observations for all six close encounters can be reproduced by our single-fluid MHD model provided that the resolution of the computational grid is sufficiently high.

[7] In this paper we present the simulation results and describe how they can be used to further our understanding of the structure of Ganymede’s magnetospheric system. The simulation provides realistic description of the magnetic environment near Ganymede that includes the plasma currents generated by the plasma interactions self-consistently. A detailed description of our MHD model is given in 2. 3 describes the simulation results including the global structure and comparisons with the observations. We conclude with discussion and a summary.

2. MHD Model

[8] The simulation code we use is based on the one developed to model the interaction between the Jovian plasma and Io [Linker et al., 1998]. Details are given by Linker et al. [1998]. Here we briefly describe the original MHD model and indicate how we adapt it to the Ganymede simulation. This code solves the following form of the MHD equations with normalized units in spherical coordinates.

[9] Here is the vector potential, and the magnetic field is calculated as = ∇ × , = ∇ × is the vector current density, ρ is the mass density, is the flow velocity, P is the plasma scalar thermal pressure, γ is the ratio of specific heats taken as and η is the resistivity. Viscous terms that are necessary for capturing sharp fronts have been added into equations (3) and (4) through terms involving the viscous stress tensor W. Details of the viscous tensor are given by Linker [1987] and here we have left W in an unexpanded form for simplicity. The typical value of the fluid Reynolds number R_e (ratio of the viscous diffusion time to the advection time) used in the simulation is 20.

[10] We solve the coupled equations on a nonuniform spherical mesh. Temporal derivatives are advanced with leapfrog time differencing combined with a semi-implicit method while centered spatial differencing is used for the terms appearing on the right hand side. The semi-implicit scheme has been used in the simulation of laboratory plasmas [Schnack et al., 1987] and solar coronal plasmas [Mikić and Linker, 1994; Linker et al., 1999; Lionello et al., 1999]. The semi-implicit method introduces a term into the momentum equation (4) that effectively modifies the inertia of the short-wavelength modes, while accurately treating the long wavelengths, enabling the time step to exceed the Courant-Friedrichs-Lewy (CFL) limit for Alfvén and magnetosonic waves. The semi-implicit scheme achieves efficiency by removing short-wavelength high-frequency oscillations. These oscillations are on the scale of the grid spacing. Since we are primarily interested in the large scale configuration of Ganymede’s magnetosphere the semi-implicit scheme is well suited to this problem. The increased time step greatly reduces the overall cost of the calculation as has been found in computations of the solar corona [Mikić and Linker, 1994; Linker et al., 1999]. Preconditioned conjugate gradient methods are used to invert the matrices arising from the implicit solve. The conditions ∇ · = 0 and ∇ · = 0 are exactly preserved.

[11] This code has the option of including the effects of ionization and charge-exchange through source terms in the MHD equations. In this paper we have turned off the ionization and charge-exchange effects in the calculation since typically the ambient Jovian plasma density is low and thus these effects do not play important role in the interactions except when Ganymede is located very close to Jupiter’s center plasma sheet, as during the G8 encounter. The mass loading effects at Ganymede will be discussed in a future paper.

[12] In order to include as much of Ganymede’s magnetotail and Alfvén wing structure as possible we set the simulation domain from 0.5 R_G to 40 R_G. Figure 1a illustrates the coordinate system in which the simulation results are shown, a Ganymede-centered Cartesian system (referred to as “GphiO”) based on the direction of Jupiter’s convective flow. In the GphiO system, is along the incident flow direction, is along the Ganymede-Jupiter vector, positive toward Jupiter, and is parallel to Jupiter’s spin axis. However, the code is written in spherical coordinates, with θ measured from the + axis and ϕ measured from the XY plane (Y > 0) of the GphiO system. This implies that Ganymede’s equator is at θ = 0°, 180° while its poles are at θ = 90° and ϕ = 90°, 270° in the simulation coordinates. The code uses a structured and stretched mesh, which is similar to that used in the Io calculation by Linker et al. [1998]. The nonuniform mesh allows grid points to be distributed finely in regions of interest. In Ganymede’s case, we are particularly interested in the region of the magnetosphere within several R_G in radial distance. Simulations presented in this paper use 131 × 132 × 128 (r, θ, ϕ) grid points. Figure 1b shows a typical grid distribution in a (r, θ) cut through the 3-D mesh adopted in our simulations. The grid provides fine resolution within Ganymede’s magnetosphere and around the low latitude magnetopause (∼ 2 R_G to 3 R_G). The finest grid resolution in the radial direction is of order 0.01 R_G ≈ 26 km and the average spatial resolution within the magnetosphere is about 0.04 R_G ≈ 100 km. As we will demonstrate in the next section, high grid resolution is necessary to capture the signature of strong currents that produce sharp rotations in the observed magnetic field near the magnetospheric boundaries. We also have the freedom of setting a nonuniform grid in the direction. Normally we put finer grids around Ganymede’s equator and coarser grids around its pole in order to resolve the boundary currents and strong gradients at low Ganymede latitudes. The grid in the direction is uniform.

[13] There are three boundaries appearing in the simulation: the core boundary, the inner boundary and the outer boundary. The magnetic flux corresponding to Ganymede’s internal dipole moment [Kivelson et al., 2002] is set at the core boundary (r = 0.5 R_G). The tangential electric field at this boundary is assumed to be zero in order to keep the magnetic flux in the core fixed at the value assumed for Ganymede’s intrinsic field. We advance the field in the interior of Ganymede between the inner boundary and the core boundary using equation (1) with = 0. This makes equation (1) into a pure diffusion equation. We make the diffusion time very long compared to the typical Alfvén timescale (τ_A = ≈ 12 s, where the Alfvén speed v_A = 220 km/s) by adjusting the resistivity η. The background value of the resistivity η inside of Ganymede and outside the ionosphere used in the calculation results in a Lundquist number (ratio of the diffusion time to the Alfvén traveltime) of 4000. At the inner boundary (r = 1.05 R_G), Ganymede’s ionosphere is treated as a finite-conductivity thin layer populated with cold dense plasma. At this boundary the plasma density and pressure are uniformly specified with fixed values. Flow stops at this boundary, i.e., = 0. The resistivity in this layer increases from the background value used far from Ganymede to the value assumed at the surface in order to mimic the finite electric conductance of the ionosphere. At Ganymede, the conductance of the ionosphere is quite uncertain. Kivelson et al. [2004] give a Pedersen conductance of 2 S while Eviatar et al. [2001] estimate the conductance to be ∼100 S due to ion pickup. In the present calculation, we set the height-integrated conductance in the ionospheric layer to be ∑ = 2 S. A mass density of 550 amu/cm³ with an average temperature of 20 eV are set as initial conditions and fixed at the inner boundary. The outer boundary (r = 40 R_G) is divided into the upstream (X ≤ 0) and downstream (X > 0) boundaries. We specify ρ, P, , and the tangential electric field _t on all points at the upstream outer boundary with constant values corresponding to those observed by Galileo for each individual pass modeled. Since Jupiter’s rotation period (∼10 h) and Ganymede’s orbital period (∼7.15 d) are much longer than the time it takes Jupiter’s corotating plasma to flow past Ganymede’s magnetospheric system, which is on the order of minutes, the background Jovian environment does not change significantly over the time interval required for the system to reach a quasi-steady state. Therefore using steady upstream physical parameters to label the boundary conditions is a reasonably good approximation. At the downstream outer boundary, a non-reflecting boundary condition based on the gas characteristic equation [Hedstrom, 1979] is applied and plasma is free to leave the simulation domain at this boundary. The gas characteristics give boundary conditions for the pressure, density and the radial (normal) velocity. Neglecting the magnetic field in the characteristic boundary conditions may cause some reflections at the outer boundary. In order to minimize the boundary effects, we place the outer boundary as far as possible (=40 R_G) from the region of interest at an acceptable computational cost and then cut the run off before the boundary effects affect the interior solution.

[14] We initially set up a uniform flow with the specified plasma density, pressure and flow velocity, and a uniform background field superposed on Ganymede’s internal field. Plasma with the assumed upstream properties flows into the simulation sphere from the upstream boundary (X < 0 hemisphere) and propagates in the direction parallel to the -axis. The physical parameters used in our simulations for each flyby are listed in Table 2. The background magnetic field is set in the YZ plane perpendicular to the flow (see Figure 1a) with specified values corresponding to the Jovian field inferred by Khurana [1997] around the encounter’s closest approach. In Figure 2, we show the model Jovian magnetic field [Khurana, 1997] including both components (B_r,III, Bθ_,III, and Bϕ_,III) and magnitude (in System III coordinates) at Ganymede’s orbit. The locations where Galileo’s six close encounters with Ganymede took place are indicated by vertical lines. The figure indicates that the external field conditions during six flybys include a large range of magnetic conditions at Ganymede’s orbit. In general, Bϕ_,III in System III coordinates, which is equivalent to B_x in the GphiO coordinates, is a minor component compared to the other two components B_r,III and Bθ_,III. As an approximation, we have neglected the B_x component and kept B_y and B_z (in the GphiO coordinates) of the background magnetic field in the simulation. For Ganymede’s magnetic field, the intrinsic dipole and the induced field from Kivelson et al. [2002] defined in the GphiO coordinates are added together to obtain the total dipole moment. As shown in Figure 2, the radial component of the background field (B_r,III) exhibits the largest variation. Given a conducting shell present within the moon, this time varying component of the external field can drive inductive currents within the shell and hence produce a time-varying magnetic moment that lies in Ganymede’s equatorial plane and points roughly toward and away from Jupiter [Kivelson et al., 2002]. Because the induction field changes for different flybys, the assumed X and Y components of the internal moment vary from pass to pass (see Table 2). Typically the Jovian plasma at Ganymede’s orbit is subcorotating [McNutt et al., 1979] and the convective flow speed inferred from the analysis of the anisotropies in ion distributions typically is about 80% of the rigid corotation speed, that is ∼190 km/s [Williams et al., 1997a, 1997b]. In most of our runs, we use a typical value of 140 km/s as the background flow speed in the reference frame of Ganymede. For one of the encounters, the G7 flyby, which occurred on the dusk side of the Jovian magnetosphere, the corotating flow speed is smaller than the typical value according to the EPD measurements [Krupp et al., 2001] and we use a value of 130 km/s for the background flow speed. Near Ganymede’s orbit, the plasma density is highly variable depending on Ganymede’s location relative to the Jovian central plasma sheet [Kivelson et al., 2004]. The ambient plasma density changes by a factor of 5 as the Jovian plasma sheet flaps over Ganymede. The average plasma number density is around 4 cm⁻³ with a range from 1 to 10 cm⁻³. We use the electron number density inferred from the upper hybrid resonance frequency measured by the Galileo Plasma Wave Subsystem (PWS) [Gurnett et al., 1996] assuming that the average mass per unit charge is 14 amu [Kivelson et al., 2004] and that the ions are singly charged. The total thermal pressure of the ambient plasma (including both thermal and non-thermal components) is set as P_t = 3.8 nPa for the G8 pass, the close encounter that occurred near the central plasma sheet, and P_t = 1.9 nPa for the other five passes which occurred away from the center of the plasma sheet. Some dimensionless parameters of the ambient flow including the sonic Mach number (ratio of the flow speed to the sound wave speed), the Alfvén Mach number (ratio of the flow speed to the Alfvén wave speed) and the plasma β (ratio of the thermal pressure to the magnetic pressure) calculated from the above physical parameters are shown in the last three columns of Table 2. The sonic Mach number M_s and the Alfvén Mach number M_A used in the simulations are less than one for all six passes. The plasma beta β typically is less than one. This is the case for five of the six passes but the G8 flyby occurred near the center of the Jovian plasma sheet and the plasma β exceeded unity.

[15] For each set of upstream conditions, we let the system evolve in time until a quasi-steady magnetosphere develops. In simulations of the Earth’s magnetosphere, a quasi-steady state usually can be achieved on timescales of order of a low multiple of the time needed for the solar wind to propagate the length of the magnetosphere. In Ganymede’s case, we normally run the code for 5 ∼ 8 min or about a factor of four to six times the typical flow convection time across the magnetospheric region (τ ∼ ∼ 90 s, for v_flow = 150 km/s).

3. Simulation Results

3.1. Global Configuration

[16] We first show the overall configuration of the magnetosphere produced in the simulation. Figure 3a is a three-dimensional view of the magnetic field lines extracted for conditions of the G8 flyby at time t = 350 s. In this case the ambient magnetic field is mainly in the − direction according to Table 2. A well-shaped minimagnetosphere with standoff distance less than 2 R_G has developed in the simulation. The whole magnetosphere is in a roughly cylindrical shape differing from the bullet-like shape of planetary magnetospheres because in a sub-Alfvénic flow the dominant external magnetic pressure exerts force perpendicular to Ganymede’s magnetic spin axis whereas planetary magnetospheres are confined by the dynamic pressure of the super-Alfvénic solar wind flow on the upstream side and only on the downstream side by magnetic and thermal pressure. The magnetosphere has a small closed field line region near the equator and a large polar cap containing field lines that link to Jupiter. The polar cap is not symmetric with respect to the spin axis. The boundary of open and closed field lines is shifted toward higher latitudes by about 5° on the upstream side and toward lower latitudes by ∼10° on the downstream side due to the effect of the magnetopause currents. The shift of the polar cap relative to the dipole axis is almost same as the shift of the Earth’s polar cap [Holzworth and Meng, 1975]. Those field lines that connect Jupiter and Ganymede’s high latitude region are tilted with respect to the background field direction in both hemispheres. These field lines lie in the Alfvén wing. A cut through the Alfvén wing structure is shown in Figure 3b. This illustration is plotted in the XZ plane, which contains the projections of the ambient field lines and plasma flow. Color contours represent the x-component of the plasma flow velocity (V_x). Also plotted as the orange dashed lines are the Alfvén characteristics, which are tilted with respect to the background magnetic field at an angle given by θ = arctan(), where v_flow and v_A are the flow speed and the Alfvén speed of the unperturbed ambient plasma, respectively. Theoretically, these lines delineate the configuration of the Alfvén wing. Our simulation results clearly show that the open field lines, which have one end on Ganymede and the other end on Jupiter, are tilted relative to the background field orientation (in the − direction) at an angle consistent with the theoretical prediction.

[17] The incident plasma flow slows down as it approaches Ganymede’s magnetosphere. Because in this case the background plasma flow is sub-Alfvénic and subsonic, no bow shock forms in front of the magnetosphere. Instead, plasmas are slowed by the interaction with compressional magnetosonic waves that propagate upstream. Some of the incident flow diverts around the magnetosphere and is accelerated on the flank side. The ambient plasma has only limited access into the magnetosphere.

[18] Ganymede’s intrinsic field is nearly anti-parallel to the external field near the equator, and consequently magnetic reconnection is the major process for the plasma and energy entry into in the system. Plasma enters the Alfvén wing via reconnection and is convected across the polar cap toward the downstream region. Tail reconnection eventually returns part of the flow back toward the moon and ejects the rest down the tail. Within the Alfvén wing, which is confined by the yellow dashed lines, the plasma flow is significantly decelerated because the footprints of the open flux tubes in Ganymede’s ionosphere move at a much slower speed than does the background flow and the disturbances propagate along the magnetic field line via Alfvén waves. In the downstream region, the flow is accelerated gradually to the background flow speed by the magnetic tension force exerted on the plasma. Furthermore, a thin current sheet extends several Ganymede radii from the moon’s surface in the downstream region, separating magnetic field lines with opposite polarities. A detailed discussion of the associated current paths in the magnetosphere will be given in the next section.

[19] Because the external magnetic field varies slowly and remains in a favorable orientation (southward in this case) for reconnection, Ganymede’s magnetosphere enables us to investigate the reconnection process in a relatively stable external environment. We use Figure 4 to illustrate where reconnection takes place in the magnetosphere. Results shown in Figure 4 are extracted from the run for the G8 conditions presented above. In this case, the isosurface of B_z = 0, which surrounds the moon, roughly lies on the magnetopause boundary at low latitudes. The field lines traced near the magnetopause clearly are of three different types: external (Jovian), closed (shadowed by the isosurface) and connected field lines. The strongly bent field lines on the upstream flanks in Figure 4 show that reconnection is not confined to the upstream nose of the magnetopause. Fast northward and southward flows associated with reconnection are found on the boundary in regions where the field lines are reconnected. There is some evidence, such as the converging of field lines on the downstream side (near the blue and red regions in Figure 4), that show that reconnection also occurs on the downstream flanks. We also confirm that in the simulation fast horizontal flows (mainly in the ± direction) associated with reconnection (not shown in Figure 4) are found near the converging field regions.

3.2. Comparisons of Model Results With the Observations

[20] In order to validate our MHD model, we compare our simulation results directly with the Galileo magnetometer observations for all six passes. Because the upstream boundary conditions are fixed throughout each calculation, the simulated magnetosphere does not vary significantly between successive time steps after the system attains a quasi-steady state. We find that the system reaches quasi-steady state after t = 270 s during which the ambient plasma has flowed a distance equal to 3 times the length of Ganymede’s magnetosphere. We then extract the components and magnitude of the magnetic field from the simulation along the Galileo trajectory to obtain the “simulated” flyby data. The results for all six close encounters are shown in Figure 5. In each subplot of Figure 5, the Galileo magnetometer observations and the simulation results are plotted together for each component and the magnitude of the magnetic field. In addition, the components and magnitude of the static superposition model results [Kivelson et al., 2002] are also plotted. For the G8 pass shown in Figure 5d, the simulated flyby data traced along the original Galileo trajectory is plotted and an additional trace shows the simulation results along a slightly shifted trajectory.

[21] From the plots, it is clear that the overall variations of the “simulated” magnetic field, including both the components and the magnitude, are in satisfactory agreement with the measured data. Changes in the background Jovian magnetic field and plasma conditions and spacecraft trajectories through different regions of the magnetosphere are well accounted for by our MHD simulations with a single set of inner boundary conditions. The magnetopause currents at low latitude are well resolved in the simulations that reproduce the rotations observed in the field orientation at the magnetopause crossings. In particular, the rapid rotations measured on inbound magnetopause crossings basically are well reproduced in the simulations despite small offsets of the locations. The discrepancy in the B_x component between observations and the vacuum superposition can be attributed largely to the field-aligned Alfvén wing current that also has been captured in our MHD simulations. For instance, in the inbound passes of the G1, G2, and G7 flybys, the observed B_x components show negative perturbations inside of the magnetosphere, which is consistent with the direction of the field bend back in the northern hemisphere. The simulation results produce consistent trends and amplitudes for the perturbations indicating that the intensities of the flowing currents in the Alfvén wing are well modeled. However, for some of the outbound magnetopause crossings, such as 07:36 UT in the G7 flyby and 08:43 UT in the G29 flyby, the field rotations are less sharp than the data indicate. The reason for this is that a fine spatial resolution in the MHD simulations is crucial in capturing intense localized currents near the boundaries. The spatial resolution of the simulation grid becomes less refined as the distance from Ganymede increases. The outbound passes occur at larger radial distances from Ganymede (around 5 R_G) than the inbound passes, which usually take place near 2 R_G. Referring to Figure 1b the grid separation at r = 5 R_G is about r = 0.26 R_G, where as at r = 2 R_G the grid spacing is r = 0.04 R_G. The grid resolution decreases with distance not only in the radial direction () but also in the other two directions ( and ) because the geometric factors in spherical coordinates (rθ and rsinϕ) depend on r. The effect of the changing grid spacing is to spread the currents at larger radial distances. Another discrepancy between the simulation and the observations is that the field magnitude is underestimated near closest approach for multiple passes, suggesting that the internal field model is probably inaccurate. For several flybys such as G1, G7, and G28, the virtual spacecraft seems to encounter the magnetosphere later and to exit the magnetosphere earlier than the observations show along the same trajectories. The time differences in the magnetopause crossings imply that the magnetosphere simulated in our MHD model is ∼10% smaller than that observed. We will return to this point later.

[22] During the G8 and G28 encounters, the Galileo spacecraft crossed the upstream side of Ganymede’s magnetosphere at low latitudes, G8 in the northern hemisphere and G28 in the southern hemisphere. The spacecraft passed closer to Ganymede (∼1.31 R_G) on the G28 flyby than on the G8 flyby (∼1.61 R_G). In the low-latitude upstream magnetosphere, the magnetopause currents enhance B_z and decrease B_x. Moreover, the magnetopause currents normally are sheet-like currents. Therefore the perturbations associated with these currents vary slowly with distance normal to the boundary. These features are clearly revealed in the G28 data in Figure 5e. Following an abrupt field rotation at 10:04 UT, the B_x component weakens and the B_z component increases dramatically as a result of magnetopause currents. The total field strength increases to nearly double that of the superposition model. After this crossing the B_y and B_z components vary very slowly for the next two minutes until the spacecraft moves closer to the moon and rapid changes are observed where the internal dipole dominates the field. These features are captured by the MHD simulations shown in Figure 5e. At the magnetopause, the field rotations and the width of the current layer produced by the simulation are consistent with the observations. However, the peak field strength near closest approach in the simulation is weaker than that observed.

[23] In the following we will focus on the G8 flyby because the field predicted for the G8 pass by the superposition model deviates from the observations considerably. In particular, the data show that the B_z component is positive inside the magnetosphere, suggesting that Galileo encountered closed field lines. Moreover, the pitch angle distributions measured by the EPD onboard Galileo show butterfly distributions in the magnetosphere, also consistent with a closed field line geometry. On the other hand, the B_z component in the vacuum superposition model remains negative and the trajectory crosses the lowest latitude open field lines. The difference between the data and the vacuum superposition model can be attributed to the compression and distortion of the magnetic field in the upstream magnetosphere [Kivelson et al., 1998]. The G8 pass occurred near the center of the Jovian plasma sheet where the plasma density is relatively high and the plasma β approaches unity. In this situation, the upstream side of the magnetosphere is highly compressed and the field lines are distorted in a manner similar to those in the Earth’s dayside magnetosphere. Strong currents at the magnetopause are required to produce the positive B_z perturbations where the spacecraft passes into the magnetosphere. Therefore the G8 flyby is a crucial pass for testing and validating our MHD model.

[24] Our simulation results extracted along the original spacecraft trajectory (plotted in green dotted lines in Figure 5d) show a slightly positive B_z within the magnetosphere. However, the modeled B_z is not strong enough and the B_x component is weaker than observed. As a consequence, the simulated field magnitude fails to match the observations. Noting that the trajectory passes very close to the upstream polar cusp, we shifted the spacecraft trajectory in the + direction toward Ganymede’s surface by 0.05 R_G ≈ 130 km and in the − direction toward the equator by 0.08 R_G ≈ 210 km. The trajectory shift of 0.08 R_G is between one and two grid points in the and direction. The red solid lines in Figure 5d indicate the simulation results extracted from the same simulation along this shifted spacecraft trajectory. Now all three components and the field magnitude from the simulation fit the observations well. The self-consistently generated currents at the upstream magnetopause enhance the B_z component and decrease the B_x component relative to the superposition model. To better illustrate the geometry of this pass with respect to the magnetosphere we superpose the intersections between the original and shifted trajectories and the XZ plane onto the color contours of plasma thermal pressure in the XZ plane (Figure 6). Also plotted are the magnetic field lines. A close view of the trajectory intersections with the simulation grid superposed is inserted to show how many grids in the numerical mesh the shift distance is equivalent to. The spacecraft moves primarily in the + direction at nearly constant X and Z position and the closest approach occurs near Y = 0. We can see that the original trajectory lies in the cusp region where the field lines are open and the plasma thermal pressure is relatively high. On the contrary, the shifted trajectory lies inside of the magnetosphere where the field lines are closed and the plasma pressure is relatively low. Even though the distance shifted is relatively small compared to the size of the magnetosphere, the field and plasma properties in these two regions are dramatically different. In other words, our simulation results of the magnetic field at the location of the Galileo spacecraft are very sensitive to the size of the magnetosphere and to any possible asymmetry of the direction of the external flow or the internal magnetic moments. Other factors can possibly contribute to the need to shift the trajectory in this case. The upstream magnetopause location is determined by the balance of the total pressure between the external plasma, including the thermal, ram and magnetic pressures, and the magnetosphere. Reducing the external pressure will cause the size of the magnetosphere to increase. The ambient magnetic field is well constrained by the magnetometer measurements together with the Jovian field model [Khurana, 1997] but the thermal pressure and the dynamic pressure are not so well constrained. The ambient plasma thermal pressure used as the simulation input is 3.9 nPa, which is about the maximum pressure in the center of the Jovian plasma sheet near Ganymede’s orbit [Kivelson et al., 2004]. We have run other cases reducing either the plasma pressure or the dynamic pressure by a half. In these runs, the upstream magnetosphere is less compressed and distorted compared to the previous run and the B_z component near the closest approach is less positive at the location of the spacecraft trajectory, although the simulated magnetosphere expands slightly due to the reduction of the external pressure. On the other hand, within the upstream magnetosphere the magnetic pressure dominates the total pressure and essentially is determined by the strength of Ganymede’s intrinsic magnetic field.

[25] Another way of expanding the size of the magnetosphere is to increase Ganymede’s internal field strength but increasing the internal moment enough to move the magnetopause through the required distance would not be consistent with measurements closer to Ganymede. However, the upstream cusp, near which the spacecraft passed through during the G8 flyby, is affected not only by the field strength but also by the orientation of the internal field. The discrepancies between the simulation results and the observations may reflect some inaccuracies in the internal field model input in the simulation. We will consider this matter in detail in the discussion.

[26] In summary, the magnetic field results from our MHD simulations basically reproduce the signatures seen in the data. The discrepancies between the MHD model and the observations near the magnetopause boundaries may arise in part because of the inaccuracy of the internal field model. Thus encouraging agreement between our MHD model results and the observations suggests that our MHD simulations can be used to improve the determination of Ganymede’s intrinsic field.

3.3. Current Flows in the Magnetosphere

[27] In Figure 7 we have used the simulation results from the G8 flyby to show the current paths. The background magnetic field in this case is mainly in the direction and the configuration of the magnetosphere is more symmetric with respect to the axis than on other flybys. We have plotted color contours of the current density (J_y) together with projected magnetic field lines in the XZ plane in Figure 7a. There are two regions of intense currents flowing in the direction (−) parallel to the ambient electric field. On the upstream side is the magnetopause current sheet, responsible for shielding Ganymede’s intrinsic magnetic field in the ambient flow and also associated with the reconnection. In the equatorial plane (XY plane), this magnetopause current sheet is found to extend across most of the upstream side. It is this strong current sheet that produces the positive B_z signature in the upstream magnetosphere discussed in the previous section. Extending several R_G on the downstream side near Z = 0 is the tail current sheet that sustains the stretched magnetic field lines with opposite polarities above and below the magnetic equator.

[28] In Figure 7b we show the field-aligned current density (J_par) together with the field line projections in the YZ plane that is perpendicular to the incident flow. In this plot, the flow is into the plane and Jupiter lies to the left side (+). This plot shows the field-aligned currents, also called the Alfvén wing currents, that produce the bend-back of the field lines within the Alfvén wing structure. These currents flow toward Ganymede on the side closer to Jupiter and away from the moon on the side away from the planet. These intense currents concentrate near the open-closed field line boundaries and on the separatrix between the connected and external field lines because of the strong flow shear near these regions. In this simulation the total downward (or upward) field-aligned current integrated over the north Alfvén wing is approximately 2 × 10⁶ A. The sum over both hemispheres yields a total downward (or upward) field-aligned current of 4 × 10⁶ A.

[29] To better understand the current paths near Ganymede, we show in Figure 7c the current streamlines in three dimensions. Color contours of the current density J_y on the XY and XZ plane and J_par on the YZ plane are also plotted for reference. At high latitudes the currents primarily flow along the magnetic field. The currents divide into three major parts as they approach low latitudes near Ganymede. The first part (color coded as orange) closes near the equator via perpendicular currents (in +) along the upstream magnetopause. An intense J_y region clearly delineates the intersection of the upstream magnetopause with the equatorial region. Here the classic magnetopause currents, which close on themselves above the cusp at relatively low latitudes, and the Alfvén wing currents that close at large distances converge and cannot be distinguished. The total current in the + direction near the upstream equator reaches 0.6 × 10⁶ A, which is about 15% of the total current flowing toward Ganymede along the Alfvén wing flux tubes. Additional currents (also color coded orange) close on the downstream side of the magnetosphere via the thin and elongated tail current sheet. This tail current sheet carries a current of 1.4 × 10⁶ A that closes about 35% of the total current in the Alfvén wing. The rest of the Alfvén wing current, depicted by the green current streamlines in Figure 7c, closes through Ganymede and its ionosphere. In the interaction between a sub-Alfvénic flow and an obstacle without a strong intrinsic magnetic field, for example Io, the relative motion between a moon and the corotating plasma drives field-aligned Alfvén wing currents, which are thought to close through the moon and its ionosphere [Neubauer, 1980; Southwood et al., 1980]. In the case of Ganymede, which has a strong internal magnetic field, the Alfvén wing currents close on themselves not only through the moon and its ionosphere, in analogy to the interaction of Io with the Jovian plasma, but also through the magnetopause and tail current sheets, in analogy to the interaction between the solar wind and a planetary magnetosphere.

[30] At Ganymede’s orbit, Jupiter’s magnetospheric field varies in both orientation and field strength. We have shown the pattern of field-aligned currents for the G8 pass during which the ambient field was nearly symmetric around the axis. Next we show in Figure 8 the variations of the patterns of the field-aligned current systems in response to the background field orientations for another two flybys G2 and G28 during which the ambient field was tilted about 45° in the YZ plane toward (in the G28 pass) and away from (in the G2 pass) Jupiter. The G8 results are also shown for comparison. The top three panels in Figure 8 present the color-coded parallel current density and the projections of the magnetic field lines in the YZ plane perpendicular to the incident flow. Compared to the G8 case, the field-aligned currents in the G2 and G28 flybys show pronounced asymmetries in both the northern and the southern Alfvén wings. In the northern wing, the parallel currents appear to be more intense in the branch on the side toward which the background B_y points. The situation is reversed in the southern wing. For example, in the G2 case (Figure 8a) in which the background B_y points toward − direction, the parallel currents are stronger on the − side than that on the + side in the northern wing. Analogous patterns can be found in the G28 case (Figure 8c) where the background B_y is in the + direction.

[31] This feature is depicted more clearly in a cut above Ganymede’s surface. The field-aligned current density crossing a plane at Z = 1.5 R_G is shown in the bottom rows in Figure 8. For the cases with significant B_y in the external field (G2 and G28), the field-aligned currents are more intense and are concentrated in narrower or broader regions depending on the sign of B_y. That the parallel currents are organized by the background B_y at high altitude can be interpreted as follows. Let us take G2 as an example. Since there is a finite B_y component in the background field (for G2 the ambient B_y is negative) in the northern hemisphere, the shear in the magnetic field between the external field and Ganymede’s internal field is stronger on the side (−) that the ambient B_y points toward than on the other side (+) because the Y-components of these two fields are in opposite directions on the − side but in same direction on the + side. The situation in the southern hemisphere can be explained analogously. In the absence of strong B_y in the ambient field, such as in the G8 case, the parallel currents are nearly symmetric on the two sides of the -axis.

4. Discussion

[32] In these simulations, we have characterized Ganymede’s intrinsic magnetic field using the internal field model from Kivelson et al. [2002] in which Ganymede’s internal dipole field was fitted by using observations from three out of the six close encounters under the assumption that the magnetopause currents contributed a constant magnetic field. Considering the weaker field strength near the closest approach in multiple simulated flybys and the need to shift the trajectory in the G8 case in order to achieve better fits to the data, we think it is possible that Ganymede’s intrinsic magnetic field may be larger than the accepted value. In future work, we will use the more accurate description of field perturbations arising from magnetopause currents extracted from our MHD simulations to refine Ganymede’s internal field model as well as the inductive response. On the other hand, the simulation results show some offsets of the magnetopause crossing time, such as the outbound magnetopause crossing in the G28 pass and both inbound and outbound magnetopause crossings in the G29 pass. Because during these passes the magnetopause crossings occurred at relatively large distances from the moon (beyond 2 R_G), an increase of 10% in the size of the magnetosphere requires about 30% increase of the dipole strength. An increase of this size can not be reconciled with the magnetometer observations from the G2 flyby whose closest approach is only 0.1 R_G from the moon’s surface. The offsets of the magnetopause crossing may be caused by factors other than the internal field. The flow vectors inferred from the ions anisotropies measured by the Galileo EPD show that the plasma flow direction can deviate from the corotation direction near Ganymede’s orbit [Krupp et al., 2001]. In the present simulations we have assumed that the background plasma flow is strictly along the corotation direction (+ direction). Some departures from the corotation direction in the ambient flow could possibly twist the shape of the magnetosphere and displace the locations of magnetopause crossings along the original spacecraft trajectory. Exploration of the effects of changing external conditions appears worth undertaking but is beyond the scope of this paper.

[33] In the present work, we have neglected mass-loading effects due to electron impact ionization and ion-neutral charge-exchange. At Ganymede, water molecules are thought to be sputtered primarily from the polar regions by energetic ions from Jupiter’s magnetosphere. Paranicas et al. [1999] have estimated the total production rate to be ∼10²⁶ s⁻¹ which is less than the ionization source rate for the sulfur and oxygen ions (∼10²⁸ s⁻¹) at Io [Bagenal, 1997] by two orders of magnitude. Both the ion-neutral reaction rate and the ionization rate associated with electron impact depend on the ambient plasma density. When Ganymede is located outside of the Jovian central plasma sheet, where the ambient plasma density is only on the order of 1 or 2 particles per cm³, the total ionization rate is so small that neglecting the mass-loading effects at Ganymede should be a reasonable approximation. On the contrary, when Ganymede is close to the center of Jupiter’s plasma sheet, where the ambient plasma can be denser than elsewhere by a factor of ∼5, mass-loading effects can be important and can modify the interaction between Jupiter’s plasma and the moon. In fact, during the G8 encounter, the magnetic field magnitude starts to increase gradually above the background at an upstream distance of about 7 R_G as the spacecraft moves toward Ganymede’s magnetosphere from the upstream flank side. The pile-up of the magnetic field could result from mass-loading of exospheric neutrals in the vicinity of Ganymede. On the other hand, changes in plasma conditions due to mass-loading, such as changes in plasma temperature, velocity and density, can affect the the size of the magnetosphere. This could be another reason why the magnetosphere produced in current simulations is systematically smaller than the observed for several flybys. Therefore mass-loading effects should be included in the simulation to describe the interaction and this work will be carried out in our future studies. However, the general conclusions are not expected to differ greatly from those drawn from the present calculations.

[34] Paty and Winglee [2006] used multifluid MHD rather than single fluid MHD to account for the magnetic field observations. They found that the upstream magnetopause location simulated by the single fluid model is ∼36% closer to the surface than that in a multifluid simulation whose dimensions correspond to the Galileo magnetic field observations. In contrast, the magnetosphere produced in our single fluid simulations with high spatial resolution is only ∼10% smaller than the observed one. It is not clear what aspect of the Paty and Winglee [2006] multifluid MHD treatment produced the ∼36% increase in the size of the magnetosphere. Finite gyroradius effects seem unlikely to cause such a large expansion of the magnetosphere. For the bulk of the ambient plasma near Ganymede with typical energy of ∼300 eV, the gyroradius of a heavy ion with mass number 14 in a 100 nT magnetic field is ∼100 km, which is less than 1% of the size of the magnetosphere. A single fluid MHD approximation should provide a reasonably good approximation for the description of the interaction of Ganymede’s system, although a multifluid MHD treatment can provide insight into the bulk motions of different ion species. In fact, we have demonstrated above that fairly good agreement can be achieved in our simulations by the use of a single-fluid MHD model using measured background properties with little adjustment. Although there are some discrepancies between our model results and the data near the magnetopause crossings, our current model is able to provide satisfactory representation of Ganymede’s magnetospheric field taking into account that the upstream conditions are not fully constrained observationally and mass-loading effects are not included. Furthermore, Paty and Winglee [2006] used supersonic and high Alfvén Mach number conditions for the ambient flow, which is different from the subsonic and low Alfvén Mach number flow conditions used in our simulations. The differences in the flow conditions could contribute significantly to the differences between Paty and Winglee [2006] single-fluid MHD and our simulation results.

[35] Stone and Armstrong [2001] modified empirically a model of Earth’s magnetosphere, with the magnetopause and tail currents included, to represent Ganymede’ magnetosphere. The fitting to the data is completed by using multiple parameters whose physical importance is not clarified in the model. The global configuration of the magnetosphere produced in our MHD simulation differs from that predicted by the Stone and Armstrong [2001] model in the sense that the magnetospheric currents, which separate the unperturbed external field and the magnetospheric field lines, close near Ganymede (within 3 R_G) in the Stone and Armstrong [2001] model. In contrast, the field lines appear to lie within the Alfvén wing in our model and the Alfvén wing currents which go through Ganymede and its magnetosphere close far away from Ganymede because the Aflvén waves can propagate a large distance in the time that the flux tube is in contact with Ganymede.

[36] In this paper, we have only shown the results extracted from specific snapshots in the simulations when the system is in a quasi-steady state. However, temporal variations, such as the strong magnetic field fluctuations observed during the G8 pass near the magnetopause crossings (around 15:50 UT and 16:02 UT), cannot be captured by a single snapshot. Although we have shown the variations of the global configuration of the magnetosphere under different background field conditions, we have ignored changes in the background properties over the time of a Galileo pass through Ganymede’s magnetosphere recognizing that changes occur on timescales of an hour. The temporal variability of the magnetosphere on short timescale (e.g., the time in which plasma convects across the magnetosphere) certainly is another important issue that needs to be investigated. Whether the reconnection is bursty or steady under very slowly varying external conditions is not clear yet and this matter needs to be investigated in great detail in the future. Because Ganymede is phase locked to Jupiter and hence there is no corotation electric field in the rest frame of plasma, the circulation of the plasma flow and the return of the magnetic flux to upstream within this tiny magnetosphere, whose space is largely filled by the moon, are expected to differ largely from the Earth’s case. We will address these issues in a future paper.

5. Summary

[37] In this paper, we have presented the results from our numerical simulations of the interaction between Ganymede and Jupiter’s magnetospheric plasma. A series of global MHD simulations with high grid resolution have been conducted by using upstream plasma and field conditions characteristic of each flyby and allowing Ganymede’s internal field to include both the steady and induced components. Our simulation results show that as a result of its interaction with a sub-Alfvénic and subsonic ambient flow, Ganymede’s magnetosphere exhibits a cylinder-shaped magnetospheric configuration with a large polar cap containing open field lines and a small closed field line region near the equator, which differs dramatically from the bullet-shape of planetary magnetospheres produced by interaction with the supermagnetosonic solar wind. In addition to the familiar structures such as the magnetopause and an equatorial current sheet, Ganymede’s magnetosphere generates an Alfvén wing that mediates the interaction of Ganymede with the plasma and ionosphere of Jupiter through field-aligned currents. During different encounters, the pattern of the field aligned currents varies according to the orientation of external magnetic field and asymmetries in the intensities of the parallel currents are organized by the ambient B_y. Direct comparisons of the magnetic field extracted from our simulations with Galileo observations yield satisfactory agreement for all six close encounters. Where our grid resolution is good, the magnetopause currents are well resolved in the simulations, such that the sharp rotations in the field orientation are consistent with the observations. The realistic magnetosphere generated in our simulations will not only help us to better understand the energetic particle behavior that has been observed by the Galileo Energetic Particle Detector but also allow us to refine Ganymede’s internal field model and to assess the inductive response.