Geophysicae Plasma and fields in the wake of Rhea : 3-D hybrid simulation and comparison with Cassini data

Rhea’s magnetospheric interaction is simulated using a three-dimensional, hybrid plasma simulation code, where ions are treated as particles and electrons as a massless, charge-neutralizing fluid. In consistency with Cassini observations, Rhea is modeled as a plasma absorbing obstacle. This leads to the formation of a plasma wake (cavity) behind the moon. We find that this cavity expands with the ion sound speed along the magnetic field lines, resulting in an extended depletion region north and south of the moon, just a few Rhea radii ( RRh) downstream. This is a direct consequence of the comparable thermal and bulk plasma velocities at Rhea. Perpendicular to the magnetic field lines the wake’s extension is constrained by the magnetic field. A magnetic field compression in the wake and the rarefaction in the wake sides is also observed in our results. This configuration reproduces well the signature in the Cassini magnetometer data, acquired during the close flyby to Rhea on November 2005. Almost all plasma and field parameters show an asymmetric distribution along the plane where the corotational electric field is contained. A diamagnetic current system is found running parallel to the wake boundaries. The presence of this current system shows a direct corelation with the magnetic field configuration downstream of Rhea, while the resultingj×B forces on the ions are responsible for the asymmetric structures seen in the velocity and electric field vector fields in the equatorial plane. As Rhea is one of the many plasma absorbing moons of Saturn, we expect that this case study should be relevant for most lunar-type interactions at Saturn.


Introduction
Since the Saturn Orbit Insertion (SOI) of the Cassini/Huygens spacecraft (July 2004), moonmagnetosphere studies have been primarily focused on two saturnian moons, Titan and Enceladus.Titan possesses a substantial ionosphere which forms a Venus or Mars-like interaction within a submagnetosonic flow.On the other hand, Enceladus's south polar source of dust and heavy gas molecules forms a large scale electromagnetic obstacle compared to the moon's size and results in an interaction which resembles that of a comet with the solar wind (Dougherty et al., 2006).Both moons are considered important plasma sources and potential drivers of the global dynamics of Saturn's magnetosphere (e.g.Gurnett et al., 2007).
Still, within the region occupied by Saturn's magnetosphere, at least 18 more natural satellites orbit the planet.Most of these moons, which sizes range from that of an asteroid (e.g.Telesto, Calypso) to that of a large moon (e.g.Mimas, Tethys, Rhea), probably have weak exospheres produced from sputtering of their surfaces by energetic particles.Such exospheres are not sufficient to slow down the incoming plasma flow or to cause a significant pile-up of the upstream magnetic field.To a good approximation, these moons are considered to be inert and are expected to have an "absorbing-body type" or "lunar type" interaction with Saturn's plasma.Such interactions have been studied extensively for supersonic plasma flows, using spacecraft observations, theory, numerical and laboratory simulations (e.g.Halekas et al., 2005;Kallio, 2005;Travnicek et al., 2005;Farrell et al., 1997;Ogilvie et al., 1996;Samir et al., 1983;Podgorny et al., 1982;Podgorny and Andrijanov, 1978;Dubinin et al., 1977;Whang and Ness, 1970;Whang, 1969).Khurana et al. (2008) and Roussos et al. (2007) have studied such interactions in the the subsonic plasma flows in Saturn's magnetosphere, using magnetic field and energetic electron Published by Copernicus Publications on behalf of the European Geosciences Union.E. Roussos et al.: Hybrid simulation of Rhea's magnetospheric interaction data, respectively.The latter study, however, focused mainly on the association of the effects of the plasma absorbing interactions with the global magnetospheric dynamics at Saturn.
Saturn's magnetosphere is an ideal environment for the study of lunar type interactions.Since Saturn's spin and magnetic axes are aligned and almost colocated, the plasma bulk velocity vector is continuously almost perpendicular to the magnetic field in the equatorial plane, where most moons orbit the planet.As a result, the bulk velocity vector, u, the magnetic field B and the resulting corotational (or motional) electric field, E=−u×B always define an orthocanonical system upstream of a moon.In addition, since many of Saturn's inert moons orbit within the plasmasphere of the planet, the upstream plasma moments can be considered quasi-stable for large time scales compared to the duration of a Cassini flyby or a simulation run.Such a configuration eliminates many complications which can otherwise mask important physical processes associated with the dynamics of a lunar wake.Moreover it makes the set up of a realistic simulation of the interaction more feasible.
For example, the solar wind impinging on the Earth's Moon can experience notable short term variations of its properties, e.g. in the interplanetary magnetic field (IMF) magnitude or orientation.Besides that, the angle between u and B is in most cases oblique.Furthermore, the macroscopic parameters of the Moon's space environment can vary substantially, as this can be located in the solar wind, in the Earth's magnetosheath, the current sheet or the tail lobes.The Jovian moons experience a similar variability as they can be positioned in or out of the plasma disk (since the axis of the Jupiter's dipole moment is tilted by 9.6 • with respect to the spin axis).In this case, the sampling of the interaction during only several flybys (as it is possible with most planetary missions) cannot be termed as representative and the comparison with simulations can be difficult.Such complications are less pronounced for Saturn's moons.
There are additional advantages of studying lunar type interactions at Saturn.The inert moons of Saturn are orbiting the planet in various distances within the magnetosphere and are exposed to different plasma environments.Therefore a comparative study can reveal the determining factors of the interaction (e.g.plasma betas, sonic and Alfvénic numbers etc.).Additionally, as within Saturn's plasmasphere the ion and electron thermal velocities are comparable or much higher than the plasma bulk velocity, the refilling of the wake should happen in a very close distance behind the moon (Samir et al., 1983).That is important as with a close flyby, both the structure of the wake and the processes refilling it can be monitored.This also requires smaller simulation boxes (or equally less convection time through the box) compared to the size of the boxes needed for simulating the wake of the Earth's Moon.Various simulations predict that the lunar wake is replenished by plasma at distances greater than 20-30 R L , where R L =1737 km, while for most moons of Saturn this should occur within a few moon radii.In this way, 3-D MHD or hybrid simulations can be set-up with increased spatial resolution.Fully kinetic 1-D or 2-D codes can also benefit in a similar way.
In the present study we utilize some of the advantages mentioned above to simulate the interaction of Saturn's moon Rhea with the magnetospheric plasma.For this purpose we use a 3-D hybrid code and a high-spatial resolution set-up.The scope of this simulation is to identify and understand the basic, macroscopic features and processes that occur in the three-dimensional wake of an absorbing body at Saturn and to exploit the importance of ion kinetic effects in the interaction.We also attempt to reproduce the signature seen with the magnetometer of Cassini during a November 2005 Cassini flyby of Rhea.Moreover, we review fundamental issues which are important in such simulations, such as the choice of boundary conditions at the surface and the interior of Rhea.
The paper is structured in the following way: in Sect. 2 we present several basic facts about Rhea and its surrounding environment and explain why we choose this moon for the simulation, instead of another Saturnian inert moon.In Sect. 3 we include a brief description of the hybrid code and the simulation set up.In Sect. 4 we present the simulation results, including the comparison with the magnetic field observations by Cassini.Finally, the summary and some additional implications and extensions of this study are presented in Sect. 5.

Saturn's moon Rhea
Rhea is the second largest satellite of Saturn.It has a radius of R Rh =764 km and orbits Saturn in a circular and equatorial orbit, with a semimajor axis of 8.75 Saturnian radii, R s (1 R s =60 268 km).Based on the value of the moon's mean volume density and the modeling of its gravitational field, it is believed that Rhea is composed by a mixture of ice and rock and therefore it should have the physical behavior of an insulator (Anderson and Schubert, 2007).The magnetic field signature from the first close Rhea flyby by Cassini to date (26 November 2005) is consistent with a plasma absorbing type of interaction with Saturn's magnetospheric plasma (Khurana et al., 2008).Equally, this means that plasma loss at Rhea's surface dominates any electromagnetic effect imposed by a sputtering-produced exosphere that probably surrounds the moon.This makes Rhea an appropriate choice for the purposes of our simulation attempts.
As mentioned earlier, Saturn's inert moons are numerous.However, selecting Rhea as the target body of the simulation instead of another inert moon is not an arbitrary choice.There are both physical and practical reasons behind this selection.
Rhea orbits Saturn at the edge the planet's E-ring (Baum et al., 1981), at location which can be seen as a "transition region" within the magnetosphere: there the magnetic field holds its dipolar configuration to a large extent but the magnetic field magnitude drops to values around 30 nT, while the plasma pressure is almost peaked (Sergis et al., 2007), leading to plasma beta values around unity.
Another important element for Rhea's interaction is the magnetospheric variability at its distance (see Table 1).Although the values of the plasma and fields parameters are quasi-stable for relatively large time scales (e.g.hours), the variance around the average values can be quite large.For instance, Sergis et al. (2007) have found that the plasma beta, β i , can vary about two orders of magnitude at the distance of Rhea, based on the abundance of energetic water group ions in the plasma.The plasma density values vary by a factor of 2-3, while the bulk plasma velocity can be between 80 and 90% of the rigid corotation.
On average, the plasma flow with respect to Rhea is subalfvénic (M A <1), supersonic (M S >1) and submagnetosonic (M MS <1).Such a combination of Mach numbers is similar to that found at the Galilean moons of Jupiter.
What is interesting is that as the Alfvénic Mach number, M A is close to 1, and given the range of values that the plasma parameters can take, extreme configurations of Rhea's space environment (e.g.high plasma density, low dipole magnetic field magnitude due to strong ring current, high plasma bulk velocity) could lead to M A and M MS ∼1.Then the conditions for a shock formation are satisfied and the presence of Rhea's exosphere cannot be neglected.We do not know, however, if such configurations are possible.
Magnetospheric variability is probably also significant for other inert moons, like for example Tethys and Mimas, that orbit Saturn at 4.89 R s and 3.11 R s , respectively.However, as these moons orbit Saturn in regions of strong magnetic field and low plasma velocity, plasma betas are always small and the interaction is always submagnetosonic.On that respect, Rhea's interaction with the magnetospheric plasma can be unique.
Adding to all the aforementioned unique characteristics of Rhea's interaction is the fact that the gyroradius of the water group ions is comparable to size of Rhea and to the width of the moon's wake, which means that ion kinetic effects can be important, especially if a weak exosphere is modelled around the moon.With a high resolution simulation set-up, even kinetic effects of protons can be resolved.Again, due to the strong magnetic field and the low plasma velocity at Tethys and Mimas, the ion gyration scales are much smaller and possible relevant effects (if any) would be very hard to illustrate.This makes Rhea a better choice for the application of a hybrid code.
Finally, as the plasma velocity relative to Rhea (∼60 km s −1 ) is higher than at Tethys and Mimas (∼34 and ∼16 km s −1 , respectively), the plasma convection time through a equal sized simulation box can be at least 2-3 times lower for a Rhea simulation than it is required for the two other moons.The same conclusion can be reached if we consider that in lunar type interactions, the sound speed is the characteristic velocity with which disturbances propagate (Samir et al., 1983).
One complication factor could come from the fact that Rhea might possess an extended dusty exosphere.This was hinted by data from the LEMMS instrument of Cassini which detected an extended energetic electron depletion region, with a scale equal to that of the moon's Hill sphere (∼15 R Rh diameter) -the region where Rhea's gravity dominates that of Saturn's (Jones et al., 2008).While we did not include any "dust component" in our simulations, we will use our results to show indeed that such a large interaction region cannot be reproduced by a hybrid code, if Rhea is simply modelled as an inert obstacle and with no mass deposited around it (in any form).

The simulation code
For the simulation of Rhea's interaction we make use of a hybrid code, where ions are treated as individual particles and electrons as a massless, charge neutralizing fluid.The code has been originally developed for the modeling the interaction of comets with the solar wind (Bagdonat and Motschmann, 2002), but has also been successfully applied to simulate the interaction of Mars and magnetized asteroids with the solar wind and of Titan with Saturn's magnetospheric plasma (Bößwetter et al., 2004;Simon et al., 2007a,b).
For the current simulation runs we use a slightly modified version of the code.The main modification is the inclusion of a resistivity term in the equations, which is necessary for a more physically correct description of the propagation of the magnetic field through the interior of Rhea.In the following subsections we include a summary of the code's basic features for the completeness of this paper.We refer the reader to these earlier studies for more detailed description of the code's main aspects.

Basic equations
The equations of the hybrid code are the following: -Equation of motion for the individual ions (single particle motion) with charge q and mass m: and where x and v correspond to the vector position and velocity of each individual ion, respectively.The E and B vectors are the electric and magnetic field vectors.
-Fluid, of massless, charge-neutralizing electrons (momentum equation) where u e describes the mean velocity of the electrons, n is the electron number density (equal to the ion number density, due to the assumption of charge neutrality) and P e is the electron pressure of the magnetospheric plasma.
-Electromagnetic fields equations: For the electric field equation we have: while the equation for the temporal evolution of the magnetic field is: The last three equations include an additional term compared to the previous applications of the code (in brackets).This has been introduced to describe the propagation of electromagnetic fields through a medium of finite resistivity, η.The inclusion of these terms is important for a self-consistent calculation of the electromagnetic field properties in the interior of Rhea.This is briefly described in Sect.3.2.
Equations ( 4) and ( 5) can be expressed in a simpler form, where the ion velocity, u i is not included and the electron mean velocity is instead included.The reason they are displayed in the given form is because in the hybrid code the ion velocity is the input parameter, while the electron velocity is derived under the assumption of charge neutrality and through the momentum equation.
-Electron population description, using an adiabatic scaling for the electron pressure: where k=2 (Bößwetter et al., 2004), n o is the initial number density of the plasma and β e is the electron beta.Using Eqs. ( 3) and ( 6), electron temperature effects can be described with this model.Because of the adiabatic description of the electron pressure, no relevant pressure term appears in Eq. ( 5) (Simon et al., 2007a).

Modeling the solid body of Rhea
Like in various geophysical applications studying induction effects within the Earth (e.g.Martinec, 1999), the solid body of Rhea is modeled as an Ohmic conductor.The last terms in Eqs. ( 4) and ( 5) emerge from the assumption of an Ohmic law for Rhea's interior, where the scalar resistivity profile η has to be specified for each simulation.As stated above, Rhea can approximately be described as an insulating solid medium.For this reason a high homogenous resistivity is applied which can be expressed by its reciprocal value, the conductivity σ Rhea =5×10 −7 S/m.This value has been chosen because a similar conductivity is estimated for the terrestrial Moon (Lipatov, 2002).Everywhere in the plasma the resistivity is assumed to be zero.In order to accommodate the solid nature of Rhea every particle impinging on its surface is removed from the simulation.
To self-consistently integrate the insulating body of Rhea into the plasma environment, Eqs.(4) and (5) are solved within the entire computational domain.Whether the moon's exterior or interior is modeled is exclusively controlled by means of the coefficients u i , n and η.This is implemented by using the resistivity profile defined above and setting the velocity u i and the reciprocal density 1/n to zero within the obstacle.

Simulation parameters
We carry out the simulation within a 3-D-box of 8×8×8 moon radii (R Rh ).We achieve a high spatial resolution using a cartesian grid with 120 cells along each box side.The coordinate system of the simulation is defined as follows: the positive x-axis is along the plasma velocity direction, the positive z-axis points is antiparallel to the magnetic field orientation and the positive y-axis points toward Saturn, opposite to the direction of the corotational (or motional) electric field, E=−u×B.
The center of the box is at [x, y, z] =[0, 0, 0], while Rhea is offset by about 1 R Rh towards the left side of the box (x<0).In this way the simulation area allocated for the wake is sufficiently large.Larger offsets can increase this area even more, but at the current stage we wanted to avoid possible interactions of the "conductive" Rhea with the left side of the simulation box.Compared to simulations of the Earth's Moon wake, our simulation box is sufficiently wider along the magnetic field direction (z-axis), as at Rhea the velocity parallel to the magnetic field (v ) is comparable to the perpendicular velocity (v ⊥ ) and this could result in a north-south extension of the wake.
In this paper the upstream plasma and fields parameters are chosen to have average values from those given in Table 1.The selected values are given in Table 2.Here we present the results from two, almost identical simulations.In the first ("ideal case") the magnetic field vector is perpendicular to the plasma bulk velocity vector.The second case is set up with non-zero B x and B y components, where This is mainly done for a better the comparison with the magnetic field data from the Cassini flyby of Rhea.
In summary, the upstream plasma velocity relative to Rhea is set to about 85% of the rigid corotation.For the ion composition description we use single ion species assuming that the plasma is dominated by water group ions (<m>=17 amu).The selected electron and ion temperatures correspond to β e of 0.09 and β i of 0.24.The magnetic field magnitude is set to ∼26 nT and the ion/electron density to 4 cm −3 .Therefore the interaction is subalfvénic, supersonic and submagnetosonic.
For the interior of Rhea the chosen resistivity value ensures a diffusion time of the magnetic field through the solid body which is similar to the plasma convection time across Rhea's diameter.The simulation time step (3×10 −3 m i /qB) satisfies the Courant condition.Normally, the simulation reaches steady state after ∼2500-3000 time steps.

Results
In this section we present the results from our simulation runs.We mainly focus on the "ideal case" scenario (B⊥u).The results from the run where B x , B y =0 are almost identical and we only refer to those for the comparison with the magnetic field measurements by Cassini (Sect.4.4).
The results from the simulation are summarized in Fig. 1.We will now discuss the structures observed in each individual parameter.

Plasma density
In all density plots (panels a, c, g) the most easily recognizable feature is Rhea's wake, identified as a dropout in the plasma density.The central region of the wake is formed due to the shadowing of the plasma flow from Rhea's volume.
Starting from the plot in panel (a) (yz plane, just behind Rhea), we see the wake of Rhea as an almost circular cross section with a diameter approximately equal to that of the moon.The central depletion region is slightly narrower than Rhea and densities are in most bins above 0.5 cm −3 .This reveals that the plasma expands relatively fast into the vacuum.
Moving on to panel (c) (xz plane), we see the structure of the wake along the magnetic field direction.As particles can easily be moved along the magnetic field lines, this plot can reveal the plasma temperature effects.The "umbral cone" of the wake becomes progressively narrower and disappears after less than 1 R Rh downstream of the moon.In the same time the plasma depletion region expands north and south of Rhea's geometrical wake.
This expansion can be explained in two different ways.The first explanation is purely geometrical and its introduced by Khurana et al. (2008).Within the radiation belts of Saturn, the ratio v ⊥ /v can be close to unity or much smaller, due to the rapid bounce motion of the charged particles along the magnetic field lines and the low plasma bulk velocity.This means that in the xz-plane the impact trajectory of the charged particles on Rhea can have a large angle with respect to the wake axis (x-axis) and this would result in a extended, north-south depletion region.For the same reason, particles within a certain energy range can even avoid absorption (e.g.sub-keV electrons), while for other populations (e.g.>keV electrons, v ⊥ /v 1) the flux tube convecting past Rhea can empty completely (Roussos et al., 2007).
The bounce motion results from the latitudinal variation of |B| along a dipole field line, which we cannot include in our simulation.Still, the effects of the bounce motion can be reproduced through the non-zero plasma temperature.This can be realized through a one-dimensional treatment of the plasma expansion into a vacuum, applied so far in many relevant studies.
The momentum and continuity equations, assuming charge neutrality and a Maxwellian distribution for the ions (as in our simulation), can be solved analytically with boundary conditions that describe a vacuum: t=0, z≤−1→n o =4 cm −3 and (−1<z≤0, t=0→n o =0).
These boundary and initial conditions define along the zaxis (perpendicular to the flow) a 100% plasma depletion inside z=−1 R Rh at t=0, where t is the convection time past Rhea's limb.The convection time can be approximated as t=(x−x limb )/u, where u is the plasma velocity, x limb is the location of the limb along the x-axis and x>x limb .These conditions are symmetrical for z>0, with the depletion occurring inside z=1 R Rh .The solution for the density is given by the following expressions (Samir et al., 1983): for (z+1)/t>−C s .For (z+1)/t≤−C s the solution is simply n(x, z)=n o .As mentioned earlier, the time t is expressed as: In Eqs. ( 7) and (8), C s is the ion sound speed, approximately 40 km s −1 .That solution tells us that for a given vertical cut through the wake, at a distance x−x limb downstream of Rhea's limb, the plasma density increases exponentially into the wake (z>−1).For for z≤−1 the density drops exponentially from its initial value (n o ).This density drop propagates outwards (z<−1) as a rarefaction wave with the ion sound speed.This exponential behavior is valid for (z−1)≥−C s t, up to the distance that the rarefaction wave has expanded.The solutions are illustrated in Fig. 2. We note that in this treatment, electromagnetic fields are not described in the equations.
As it can be seen in Eq. ( 7), the ratio of the plasma velocity to that of the sound speed, or equally the sonic Mach number, M s , will define how fast the wake will expand along the z-axis.For small M s this expansion will be rapid and will occur very close to the moon (as we see from our simulation results), while for large M s this expansion will occur at very large distances downstream of the absorbing body (e.g. as we see at the Earth's Moon).This description gives similar    shows the north-south configuration of the plasma density, velocity and electromagnetic fields, while the last panel of plots (z=0 cut/letters G to J) includes the same parameters, projected on Rhea's equatorial plane.In all panels, a small part of the simulation box is not shown, as it contains uninteresting regions (e.g.far upstream of Rhea).We also note that all vector parameters are plotted in the rest frame of Rhea.In several cases, perturbations in the vector quantities are small and the vectors seem unchanged compared to the background vector value and orientation (e.g.panel C).In this case, the perturbation is mainly ralized in the color coded vector magnitude.Additional plots in the plasma rest frame or with the background magnetic field subtracted, which reveal small perturbations better, are shown in Sect. 4.
results with the geometrical approach introduced by Khurana et al. (2008).
We can directly compare the simulation results with this simple model.A comparison is shown in Fig. 3.As it can be seen, the trend of the analytical solution is in good agreement   7).The position up to which the rarefaction wave has propagated is given by −(C s t+1).The dependence of the solution from the sound speed, tells us that the wake structure in the xz plane is determined by the plasma temperature.We note that in most previous studies the wake boundary is set at z=0 and the limb location at x limb =0, which results in a slightly different form of Eqs. ( 7).Since we want to in compliance with our simulation coordinate system, this slight modification in the equations has been introduced.
with our simulation results.The main difference is that the density depletion in the rarefaction region is slightly slower compared to what the analytical solution predicts.This probably results from the simplicity of the one-dimensional approach, and the lack of consideration of the electromagnetic fields in the equations.
What is well predicted is the distance up to which the rarefaction wave has propagated.This good agreement is found at almost all distances downstream of Rhea, up to x∼3 R Rh , where the the expansion extends beyond the limits of our simulation box.An alternative way to investigate the expansion width is by comparing the opening angle of the rarefaction region with that of a Mach cone, at the given sound speed C s in our simulation.The opening angle of a Mach cone, ϑ is given by the expression: The same equation can be expressed as a function of the sonic Mach number, M s , if we divide the nominator and the denominator by C 2 s .From Eq. ( 9) we get that ϑ should be ∼43 • .The opening angle of the rarefaction region in our simulation (panel c, Fig. 1) is about 45 • , almost in excellent agreement with the expected value.All these results show that the sonic Mach number is important for lunar-type interactions and that the one-dimensional approach of the plasma expansion into vacuum is fundamentally valid, at least for  7) solutions with the simulated data, for a distance of 1.1 R Rh downstream of Rhea's limb.The scattering of the simulation data points is due to numerical noise, deriving from the small number of particles per simulation cell (about 15 particles per cell on average).
the wake evolution in the plane that contains the bulk plasma velocity and the magnetic field (xz-plane in our simulation).
According to Eqs. ( 7), the wake expansion will continue until the density gradient across the z-direction will be zero.This theoretically happens at infinity, but with simple calculations we can find that very small gradients are achieved around 7-8 R Rh downstream of Rhea's limb.At these distances, the resulting plasma density is equal to ∼40% of the initial ambient plasma density, n o .This is of course unrealistic, otherwise charge particles would be depleted all over Saturn's magnetosphere.Practically, additional mechanisms will tend to bring the density values up to n o (e.g.magnetospheric diffusion, electromagnetic field perturbations, instabilities), which is actually what we see in our simulation.Given that, and considering also the variability of the space environment at Rhea's distance we would expect that the plasma wake can extend to a maximum distance not much more than 10 R Rh downstream of the moon.It is interesting to note here that the picture is different in energetic electrons as the wake refilling is only driven by magnetospheric diffusion, which is a relatively slow process, even at Rhea.Energetic electron wakes have been observed even 60 R Rh downstream of Rhea (Roussos et al., 2007).
Moving now to the density profile in the equatorial plane (xy-plane, Fig. 1, panel g), we see a much different wake configuration compared to that in the xz-plane.The plasma depletion does not expand much more beyond Rhea's geometrical shadow.It progressively becomes narrower and one can easily infer from the plot that it would not extend much more than 6 R Rh downstream.No similar temperature effects to those in the xz-plane can be seen in this projection.
This can be easily explained considering the single particle motion: from the maxwellian distribution of ions, we select an ion with energy greater or lower than the energy due to its bulk motion.The velocity vector for this ion is v i =v+v th , where v is the bulk velocity along the x-axis and v th an additional, random velocity vector with components v th,x , v th,y and v th,z .These additional components will give to the particle a non-zero velocity with respect to the convecting magnetic field lines.With the magnetic field oriented along the negative z-axis, the solution of the equation of motion tells us that the ion will move "freely" on the north-south direction, with a constant velocity v th,z .The perpendicular thermal velocity components v th,x and v th,y (or simply th,x +v 2 th,y ) will cause the ions to gyrate around the magnetic field lines.
This means that ion motion along the xy-plane is restricted by the magnetic field, while along the z-axis it is not.Therefore, following the geometrical approach by Khurana et al. (2008) we can understand why the plasma wake can only expand north-south and not east or west of Rhea.Particle cross-field motion is only possible in the presence of an electric field (as for example the corotating electric field which is responsible for the plasma bulk motion perpendicular to the magnetic field lines).The one-dimensional approach that we used previously is therefore not appropriate for the description of the wake structure in the xy-plane, as the constraints in the particle motion due to the magnetic field are not included in the equations.
The difference between the density wake structure along the magnetic field lines and perpendicular to it is better illustrated in a yz-cut through the simulation box, 2.4 R Rh downstream of Rhea's limb (Fig. 4).The wake has an almost rectangular shape and is at least four times wider along the z-axis.At the Earth's Moon this difference is much less pronounced in ions (due to the high sonic Mach number).It is more clear in electrons, which have a much higher thermal velocity than the ions (Bale, 1997).Jones et al. (2008) show that during the flyby of Rhea on November 2005, an extended energetic electron depletion was seen (∼15 R Rh wide) along a trajectory which is practically parallel to the y-axis.Figure 4 shows that such extended depletions cannot occur along the y-axis, unless additional absorbing material is present.If they can occur without the necessity of absorbing material around Rhea, that would require processes that cannot be described by our hybrid simulation approach.
Due to the gyration motion discussed previously, ions can impact Rhea even if their guiding center of motion does not intersect the moon.Given the ion temperature in our simulation, we can estimate a maximum gyration radius for the upstream plasma of about 300 km, or 0.4 R Rh (∼6 simulation cells).Still, upstream ions with such large gyroradii comprise only a small part of the total ion distribution and have small contribution to the total density.Therefore, any depletion region upstream of the moon, or any widening of the wake due to gyroradius effects is rather small and has insignificant effects in our simulation.Finally, the observation of the wake narrowing downstream of Rhea, in the xy-plane can be explained by two possible mechanisms: diffusion or plasma convection due to field disturbances.The first explanation is not very likely to be the answer, since diffusion alone tends to make the depletion region shallower and broader.The second possibility seems more appropriate, given the various structures we see in the electric field in the xy-projection.These are discussed in more detail in Sect.4.3.

Magnetic field
The magnetic field plots in the overview Fig. 1 are shown in panels (b), (d) and (h) (x=0, y=0 and z=0 cuts, respectively).As in the previous section, we begin first with the field configuration in the yz-plane, just behind Rhea.
Here we see two main features: a central region where the magnetic field magnitude increases by about 10% of the upstream value, and two almost similar field dropout regions, at the side of the magnetic field enhancement.The central depletion region enhancement region develops as a consequence of the requirement for total pressure balance.Since plasma pressure is lost on Rhea's surface, a field compression balances this loss in the wake.Because of ∇B=0 a field decompression (or rarefaction) occurs on each side of Rhea's cavity.
Fig. 5. Magnetic field perturbations in a the yz-plane, just behind Rhea.The direction of the plasma flow (and the wake) points out of the page.Here the background field has been subtracted to enhance visibility of the compression and decompression regions of the magnetic field.The blue circle indicates the location of Rhea's geometrical shadow.The length of the arrows (arbitrary units) is indicative of the perturbation amplitude.The peak amplitude of the residual By in this simulation is about 1.5 nT.The closed magnetic field loops and the inferred current flow vectors are also sketched.Similar instructive diagrams of the residual magnetic field have been constructed based on data from laboratory simulations of the solar wind interaction with unmagnetized or weakly magnetized obstacles (Podgorny et al., 1982;Dubinin et al., 1977).Fig. 5. Magnetic field perturbations in a the yz-plane, just behind Rhea.The direction of the plasma flow (and the wake) points out of the page.Here the background field has been subtracted to enhance visibility of the compression and decompression regions of the magnetic field.The blue circle indicates the location of Rhea's geometrical shadow.The length of the arrows (arbitrary units) is indicative of the perturbation amplitude.The peak amplitude of the residual B y in this simulation is about 1.5 nT.The closed magnetic field loops and the inferred current flow vectors are also sketched.Similar instructive diagrams of the residual magnetic field have been constructed based on data from laboratory simulations of the solar wind interaction with unmagnetized or weakly magnetized obstacles (Podgorny et al., 1982;Dubinin et al., 1977).
Several interesting features are visible: the central enhancement follows the trend of the density depletion and extends north and south of Rhea's cavity.Along the y-direction it is confined within the moon's diameter.The rarefaction regions are also extended north-south.In the y-direction their spatial extent reaches up to the edge of our simulation box.What is also interesting is that the magnetic field decrease is asymmetric with respect to the pointing of the corotational electric field, E.More specifically, it is stronger in the direction that E is pointing (y<0).In the same time, the peak of the magnetic field enhancement is slightly shifted in the opposite direction (y>0).
As magnetic field disturbances are only ∼10% of the upstream value, changes in the field orientation can be realized more easily if the background field is subtracted, as in Podgorny et al. (1982); Dubinin et al. (1977).This is shown in Fig. 5.The two main regions are also easily identified in this plot and the asymmetry along the y-direction is also clear.At the equator (z=0), opposite field perturbations in the ydirection cancel each other out, leaving only a perturbation in B z .Fig. 6.Magnetic field perturbations in a the xz-plane.The ground field has been subtracted to enhance visibility of weak disturbances.The peak amplitude of the residual B x is ab nT.Also here, positive and negative B x residuals cancel out o equator, leaving only a disturbed B z .Fig. 6.Magnetic field perturbations in a the xz-plane.The background field has been subtracted to enhance visibility of weak field disturbances.The peak amplitude of the residual B x is about 2 nT.Also here, positive and negative B x residuals cancel out on the equator, leaving only a disturbed B z .
The magnetic field residual vectors form closed loops, indicating the presence of a current towards the wake for y>0 (J in ) and a current that points away from the wake for y<0 (J out ).This is consistent with the idea that diamagnetic currents flowing along the wake boundaries and maintain the magnetic field configuration in and around the wake.
Switching to the xz-plane (Fig. 1, panel d), we see a magnetic field enhancement which correlates well to the spatial configuration of the plasma density depletion (Fig. 1,  panel c).This is a direct consequence of the fact that pressure balance has to be achieved everywhere.
Subtracting the background field, as before, we can see the perturbations in B x and B z clearly (Fig. 6).The magnetic field lines are drawn into the wake by the plasma expanding into the cavity.The magnetic field perturbations are usually expected in B y and B z , but not in B x (Khurana et al., 2008).However, in our case the perturbation of B x is rather significant.
At Rhea, plasma betas are close to unity and the flow is relatively slow, even a weak acceleration of the plasma e.g.along the −x direction can lead to notable changes in B x .For the Earth's moon for instance, where the sonic Mach number is greater than 10, the peak disturbance in the B x was only ∼2% of its upstream value, as it was found by Kallio (2005).
The last magnetic field plot shows |B| projected on the equatorial plane (z=0 cut, Fig. 1, panel h).As in the equatorial plane, B x and B y positive and negative perturbations cancel out, what we actually see is the change in the B z component.The B z reaches a peak value of about 28.5 nT, just behind Rhea, a value that drops further downstream.The rarefaction regions, and the asymmetry in the amplitude of the magnetic field dropout along the y-axis are also well illustrated in this plot.
A peculiar feature is the weak magnetic field pile-up upstream of Rhea (seen also in panel d, Fig. 1).This magnetic field enhancement is created due to an induction process, which is in detail discussed by Sonett and Colburn (1968).It can be summarized as follows.Within the interior of a resistive body which is surrounded by a constant plasma flow, an electric field arises which is equal to the motional electrical field E=−u×B (see panels f and j, Fig. 1).As a consequence of the Ohmic law, currents flow inside Rhea's solid body whose density is j=−σ Rhea u×B.The current path is closed at neighboring plasma regions and a magnetic field is induced, which is toroidal with respect to the motional electric field.At the dayside of Rhea this toroidal induction mode amplifies the external magnetic field as both fields are parallel and a pile-up in the total magnetic field results.To counteract this effect, higher resistivity values can be chosen to diminish the current and hence the induced toroidal magnetic field, but this translates into longer simulation runtime.As this weak field pile-up does not change the big picture derived from our results, we did not use higher resistivity.

Velocity and electric field
In the overview Fig. 1, we plot the plasma velocity and electric field only in the xz-and xy-planes (panels e, f, i and g).A plot of the velocity on a yz-projection is shown in Fig. 7.
We can see that just behind Rhea the plasma expands towards the center of the cavity.The spatial extent of the acceleration region is not much greater than Rhea's cross section.Further downstream the region becomes broader along the the z-axis.The non-zero v z can be easily explained by motion of plasma along the field lines, but the observed non-zero v y requires an electric field perturbation in the xy-plane.This perturbation should be associated to the presence of diamagnetic currents, inferred from our analysis in Sect.4.2 (see also Fig. 5).
The velocity plot in the xz-plane (Fig. 1, panel e) shows clearly a symmetric north-south expansion of the plasma into the cavity.Northward and southward expanding plasma streams should normally form counterstreaming populations  and set up a two stream instability, that could enhance the refilling of the plasma void.This has been observed in the lunar wake, in both data and simulations (e.g.Ogilvie et al., 1996;Farrell et al., 1997).
Transforming the same data in the plasma rest frame (Fig. 8), we can see that this expansion actually takes place along the disturbed field lines (compare this figure with Fig. 6).The expanding plasma does not seem to enter the wake in the form of accelerated beams (as has been identified at Earth's moon).Rather than that expansion takes place all over the flux tube connected to Rhea and Rhea's wake.
The plot for the E x and E z components is almost featureless but in panel (f) of Fig. 1 we can barely see a slight reduction of the electric field magnitude superimposed on the extended north-south plasma depletion region.This probably results from the slight reduction of v x (in the reference frame of Rhea) in this cutting plane, which introduces a perturbation in E y .
As along the y=0 cut electric field disturbances are insignificant and the wake refilling is driven by the expansion of the along the magnetic field lines, we can test the onedimensional, analytical approach that we applied in Sect.4.1 for the plasma density.Here, the resulting solution for v z is:   This equation practically tells us that v z should increase linearly into the wake, but as Samir et al. (1983) point out, the predicted extreme velocity values for regions of very low density are not physically valid.
Our results show the linear increase of the v z component into the wake, however the slope of this increase is much different than the slope defined from the simulated points (Fig. 10).Since the expansion occurs along the disturbed magnetic field lines, that have also a non-zero B x component, we also considered v = v 2 x +v 2 z as the characteristic expansion velocity (diamond points), but this introduces only minor corrections (v x is the x-velocity in the plasma rest frame).Again, what is well predicted is the location up to which acceleration into the void occurs, which shows undoubtly that disturbances propagate with the ion sound speed.
The differences observed between the one-dimensional approach and our simulation outline once more the importance of electromagetic fields in the plasma expansion and the complexities introduced by the three-dimensional approach.As mentioned earlier, electrostatic instabilities can also be important.
For instance, Halekas et al. (2005) found that if kappa functions are used to describe the distribution of electrons and ions, the analytical expressions for the density, the velocity and the potential drop deviate from the isothermal Fig. 9. Velocity and electric field in the plasma rest frame, projected on the equatorial plane (top and bottom panel respectively).10) solution (blue line) with the simulated data, for a distance of 1.2 R Rh downstream of Rhea's limb.Triangles show the variation of v z from our simulation results, while diamonds show the velocity variation when the residual v x is considered.
solutions of Samir et al. (1983).Several studies also discuss and emphasize the importance of the magnetic field (or the plasma betas) for the plasma expansion in a magnetized medium (Huba et al., 1992;Gisler and Onsager, 1992;Gisler and Lemons, 1989).
The most interesting structures in the velocity and the electric field are seen in the equatorial projections (Fig. 1, panels i, j).The plots show that on the one side of Rhea's limb, an increase of the electric field magnitude (panel j) and rather weak perturbations on the other side.This increase is much stronger for y>0.Further downstream the electric field increase is seen only for y>0, while for y<0 a dropout of |E| is found.
This hemispheric asymmetry is more clear in the velocity plot, where for y>0 plasma is accelerated and for y<0 the total velocity decreases.This probably explains the electric field enhancement for y>0 and the dropout for y<0.Several test runs were carried out with lower or even zero electron temperature to investigate whether this electric field enhancement is due to the electron pressure term in Eq. ( 4), but the overall picture did not change.This suggests that the observed effect results from the j×B force on the ions.This better explains why this electric field enhancement is seen only in the equatorial plane.A similar asymmetry has been reported in the three-dimensional lunar wake simulations by Kallio (2005).
Replotting these parameters in the plasma rest frame (Fig. 9), we find that close to Rhea, the electric field enhancement drives the plasma directly towards the cavity center.Further downstream, the asymmetric velocity structures seen in the rest frame of Rhea, appear as a plasma circulation pattern in the plasma rest frame, where two sectors of oppositely moving plasma are visible (Fig. 9, top panel).These structures are consistent with the diamagnetic current pattern inferred from Fig. 5.The associated electric field disturbances are not easily seen in this representation, except for regions close to Rhea (bottom panel).
There is an interesting analogy in our simulation results with experimental and simulation results presented by Hurtig et al. (2003), where they studied structures developing after injecting plasma beams in a curved magnetic field.The authors report the formation of a similar diamagnetic current system (e.g.see Fig. 6 from that paper), as well as asymmetric current structures in the "transition region" of their experiment (region where the field becomes curved).These results are partly attributed to the induced electric field that arises when the injected plasma (in our case plasma refilling the wake) experiences a time-changing magnetic field in the "transition region" (in our case the magnetic field change in the wake).The authors also highlight the importance of the j×B forces in this interaction, which we also believe that explains the effects seen in the equatorial plane in our simulation.
Finally, it is also interesting to note in Fig. 1 (panels f, g) that the electric field magnitude in the interior of Rhea is self-consistently calculated approximately equal to the external electric field, as expected (Sonett and Colburn, 1968).While predefining the electric field vector in the moon's interior is a reasonable approximation (e.g.Kallio, 2005), equaling it with the upstream corotational electric field might not be completely accurate, given the structures in the electric field that arise in the region of the obstacle.

Comparison with the magnetic field data
The Cassini spacecraft made its first close flyby of Rhea to date on November 2005, where a magnetic field signature, typical for a plasma absorbing body, was observed by the spacecraft's magnetometer (Dougherty et al., 2004).Cassini's flyby trajectory was almost perpendicular to the plasma wake of Rhea (parallel to the y-axis and about 0.65 R Rh downstream of Rhea) and slightly south of Rhea's equatorial plane (z∼−0.3R Rh ).A more detailed description of the flyby and the observations can be found in Khurana et al. (2008).
Before the encounter the magnetometer measured small values in B x and B y , and therefore for a more direct comparison with our results, a simulation was set up with non-zero B x and B y components, as we also describe in Sect.3.3.The overall results from that simulation do not differ from those presented for the "ideal case" scenario, since B x and B y are still much smaller than B z .
The comparison is shown in Fig. 11.As we can see there is a very good agreement in all three components of the magnetic field.Both the amplitude of the perturbations, as well as the the spatial extent of the perturbed regions are reproduced, which means the the main features of the interaction can be described by our code.The agreement in the asymmetric amplitude of the B y perturbation for y>0 and y<0 is especially interesting (panel b), meaning that asymmetries along the y-direction seen in our simulation are realistic.
Although the exact agreement of the data and the simulation is beyond the scope of this study, it is interesting to discuss the origin of several differences seen, such as in the exact increase of |B| in the wake (panel d).The measured peak increase in |B| is about ∼30% higher to what we see with in our simulation.There are two main sources of this discrepancy.
Probably the most important one is the choice of the plasma betas, which at Rhea can vary significantly.The plasma betas practically define how much plasma pressure is lost on Rhea's surface.For example, during several test runs with the plasma temperature set to 200 eV, the peak in |B| was above 29 nT in the wake, which is in much better agreement with the observations.The same can be achieved by varying the plasma density.
A less important source of the discrepancy is the choice of a finite conductivity for Rhea's interior.This leads to a small magnetic field pile-up upstream of Rhea which opposes to the energy loss due to plasma absorption.Therefore the required compression of |B| in the wake will not be as strong as when this pile-up is not present.
Another difference seen is that the measured In green color we also show the magnetic field signature on an equatorial trajectory in the wake along the y-axis, from the "ideal case" run.Fig. 11.Comparison of the simulation results with the magnetometer data from the Cassini flyby.In green color we also show the magnetic field signature on an equatorial trajectory in the wake along the y-axis, from the "ideal case" run.
the dipole field strengthens.In our simulation the dipole variation is not included in the background field.
As we see from the data, within the range of the simulation box, the measured |B| increases about 1-1.5 nT (or about 4% of the background value).We cannot define whether this variation induces any additional effect in the interaction.An additional relevant effect is that of the gradient drift, which in the Saturnian magnetosphere accelerates ions and decelerates electrons.Furthermore, the corotation velocity also varies slightly along the y-axis.All these could maybe enhance the observed asymmetries along the y-direction, but any such enhancement will probably be insignificant, given that the simulated area is rather small.In Fig. 11 we also plot the simulated magnetic field components from a trajectory through the wake, but along Rhea's equatorial plane (z=0) (green line).The data are retrieved from the "ideal case" run, where B x and B y are initially zero.What we want to point out is that with such a trajectory no variation in B x and B y is seen and the only indication about the nature of the interaction is in B z .While in most cases the increase in B z is sufficient to reveal that plasma absorption dominates, B z can also be enhanced within the current system region of an ionosphere.Therefore for very low altitude, equatorial flybys, the signature of B z alone could be misleading.For non-equatorial, downstream flybys, the shape of the B y variation can reveal the type of the interaction (Khurana et al., 2008).On the other hand, it is interesting to note that in the simulation run where B x and B y were non-zero, it was not possible to find a trajectory where B x and B y where undisturbed, as in the case of the "ideal-case" run.

Kinetic effects and phase-space diagrams
In this study we did not present simulations where a weak exosphere is modelled around Rhea, and therefore kinetic effects cannot be easily realized in our results.On the other hand, these results show that our code can give a valid description of the wake of an absorbing body, and its use can be extended to more complex scenarios.
For illustation purposes only, we show in Fig. 12 the escape of heavy ions from a test run where a weak exosphere was included.No differences were found in the resulting wake structure or the configuration of the electromagnetic fields, compared to what was presented in the previous sections.The magnetospheric plasma did not "notice" the presence of the low density exospheric ions.The corotational electric field penetrated down to the surface, accelerating the ions in spiral escape trajectories.The characteristics of these trajectories (spiral width and height) are in agreement with the analytical description given by Simon et al. (2007a).The initial acceleration of the exospheric ions along the electric field direction is also seen.
Apart from the exospheric ion escape, particle phase-space plots be constructed with detailed information about the particle distributions in each region.Here we show several examples from the "ideal case" simulation run, where a particle phase-space plot was converted to an ion spectrogram,  What stands out from this spectrogram is the peculiar shape of the ion wake.Given that the mean energy at each position defines the particle velocity, the shape of the ion wake is due to the fact that for positive and negative y, the ions are accelerated and decelerated, respectively.Similar wake structures in ion spectrograms have been observed at the moon Ogilvie et al. (1996), and have been reproduced at various lunar-wake simulations (Kallio, 2005;Birch and Chapman, 2001;Travnicek et al., 2005).
The width of the ion distribution represents the ion temperature at each point.In the near wake region, distributions become narrower, suggesting a drop in ion temperature.In  the central wake, the distribution is broad, consistent with an increase in temperature.At larger distances downstream, the distribution retains again its Maxwellian form and ion temperature drops.We note that the treatment by Samir et al. (1983) (Eqs. 7, 10) is based on an isothermal plasma expansion in the cavity.
Several additional interesting features can be found through the phase space diagrams: the top panel of Fig. 14 shows velocity phase-space plots of v z along a hypothetical trajectory that is parallel to the x-axis and crosses Rhea through its center.We can clearly detect that downstream, two distinct populations exist, with non-zero v z , while populations with v z ∼0 are filtered by Rhea.
Initially, one maxwellian with v z =0 could describe the distribution, but downstream, two maxwellians are required with v z =±30 km s −1 .As Rhea orbits in a dipolar region, these accelerated ions will become trapped populations mirroring at low magnetic latitudes and will continue drifting around Saturn, which means that plasma absorbing moons tend to make angular (or pitch angle) distributions in the magnetosphere more field-aligned.The presence of such populations is further enhanced by the electron and ion interaction with Saturn's E-ring dust and the associated neutral gas cloud.
An ion spectrometer with a 4 π angular coverage should therefore detect two ion populations along the magnetic field lines, during an equatorial wake crossing.Similar observations in the lunar wake have been reported by Ogilvie et al. (1996).
What is also interesting is that particles with high v z are the ones that first appear downstream of Rhea.This is consistent with the geometrical concept for explaining the phasespace distribution in the wake, introduced by (Khurana et al., 2008), as an alternative to the analytical treatment of the problem by Samir et al. (1983).A simple illustration of this   approach is shown in the bottom panel of Fig. 14.Particles with high v z that just avoid impacting on Rhea appear first in the wake due to the large angle of their trajectory with respect to the x-axis.

Summary and outlook
The present study focused on the magnetospheric interaction of Saturn's moon Rhea using a three-dimensional, hybrid plasma simulation code.The space environment of Rhea (low bulk velocity, high thermal velocity) is ideal to emphasize the relative importance of plasma parameters and to understand basic physical process that occur when plasma expands into a vacuum.
The main results of our case study are summarized below: -The plasma wake of Rhea probably does not extend much more than 10 R Rh downstream of the moon.The extension of the wake depends on the ion sonic Mach number.The length of this extension can be approximated by requiring that the density gradient along the x-axis (calculated from Eq. 7) is close to zero, given the good agreement of this analytic description by Samir et al. (1983) with our simulation results.
-Disturbances in Rhea's wake propagate with the ion sound speed, C s , as described by Samir et al. (1983).
-Due to the high thermal velocity with respect to the plasma bulk velocity, the depletion region behind Rhea extends along the magnetic field.The magnetic field magnitude increases in all this region to maintain the total pressure.The extension of the depletion region and the response by the magnetic field has been seen during the flyby of Saturn's moon Tethys by Cassini, below that moon's equatorial plane (Khurana et al., 2008).In the equatorial plane the density wake is restricted by the magnetic field within Rhea's geometrical shadow.
-Although several studies consider that B x in the wake should remain undisturbed, the magnetic field downstream of Rhea varies in all three components.This is also verified by the comparison with the Cassini magnetometer data.The B x perturbation in the wake probably reflects the difference between supersonic and subsonic lunar-type interactions.
-Positive and negative perturbations in B x and B y cancel out on Rhea's equatorial plane (z=0), leaving only a perturbation in B z .
-A decompression (rarefaction) of the magnetic field occurs on the side of Rhea's geometrical wake.The amplitude of the magnetic field dropout is asymmetric with respect to the direction of the electric field.This asymmetry is also seen in the magnetometer data (e.g. the B y component, Fig. 11, panel b).
-On the xz-plane the wake refilling occurs by plasma expansion along the field lines, while on the equatorial plane plasma is accelerated into the vacuum under the action of j×B forces.
-A diamagnetic current system is developed in the equatorial regions in the wake.The presence of this current system is revealed by the plasma residual velocity distribution on the equatorial plane (Fig. 9) and by the residual magnetic field projected on the yz-plane (Fig. 5).
-The asymmetric structures seen in the equatorial plane probably result from the action of the j×B forces on the ions.
-The interaction of a mass absorbing moon within a planetary magnetosphere tends to enhance the field aligned component of the particle angular (or pitch angle) distributions. - The assumption of a finite resistivity for the interior of Rhea can adequately and self-consistently describe the behavior of the electromagnetic fields in this region.
Besides these results, there are still many more issues that could be interesting to investigate in future studies.For instance, the structure of Rhea's wake under the different plasma environments, or any effects that result from the presence of a weak exosphere when the fast magnetosonic Mach number is greater than one, would be worth investigating.It would also be interesting to model how the system responds if there is a change in the upstream plasma configuration (e.g. during a plasma injection; see Hill et al., 2005).The latter scenario is theoretically very interesting: during a fresh injection event the flux of energetic ions (>50 keV) increases significantly, and these ions carry more than 70% of the total plasma pressure (Sergis et al., 2007).In the same time, most of these ions can escape absorption from Rhea due to their large gyroradii ( 1 R Rh ), their large bulk velocities and long bounce periods.As a consequence, the pressure drop in the wake will be small and therefore the low energy plasma wake could be sustained for large distances downstream, even greater than the 10 R Rh we predict with our simulations.
Through a fully kinetic approach, such as the one by Birch and Chapman (2001), high frequency or small scale phenomena (down to the Debye length, ∼20-30 m for Rhea) could also be studied.In particular, effects resulting from the fast expansion of electrons into the wake and the deviation from charge neutrality could play a significant role in the plasma dynamics downstream of Rhea.Such effects cannot be described in a hybrid code, due to the unavoidable assumption of charge neutrality.For a non-neutral plasma distribution in Rhea's wake, an additional term in the current will be added.In total, the current in this region will be given approximately by the following expression: In a hybrid code, the first term on the right hand side is zero, as n e =n i and cannot be described.If the second term, which can be calculated from our code, is more important then the assumption of charge neutrality will not affect significantly the overall results.
On the other hand, if this term is greater than the second term, then close to Rhea ion acceleration into the cavity from all sides will dominate and will form additional structures in the ion and electron velocity fields compared to what our simulation predicts.These additional structures will also be related to an ambipolar electric field resulting from a nonzero (n e −n i ) term.Kallio (2005) described this field analytically in his lunar-wake hybrid simulations and realized only small differences in his results compared to the case where this electric field was neglected.Still, Birch and Chapman (2001) have shown that a self-consistent calculation of the charge separation effects, results in much more complex structures in the wake that the analytical description of Kallio (2005) probably cannot reproduce.
Finally, increasing the resolution to the scales of the Debye length would also be interesting for the investigation of surface currents and charging processes (Farrell et al., 2007).Surface charging at Rhea could reach high values (kV or higher), as Rhea orbits in the ring current region of Saturn and it is therefore continuously exposed to energetic particles (see for instance Halekas et al., 2006).What makes such a region unique for investigation is also the fact that Rhea's dayside (where surface photoemission currents are generated) does not coincide with the hemisphere where Saturn's plasma is impinging, unless Rhea is situated on the dusk sector of the magnetosphere.For example, during the flyby of Rhea by Cassini, Rhea was at a local time of ∼12:00 and therefore photoemission currents were also created in the plasma wake, a case impossible for the Earth's moon.
Overall, the results presented in this study describe fundamental processes that occur downstream of plasma absorbing obstacles and are important for understanding the measurements that have been, or will be acquired by the plasma and fields instruments of Cassini during relevant moon flybys.Rhea is one of the many plasma absorbing moons of Saturn within the region of the planet's plasmasphere.Therefore our results should be relevant for most lunar-type interactions at Saturn.Given the fundamental nature of the "plasma vacuum behavior" problem, the Saturnian moon wakes represent an ideal physical laboratory which should be carefully investigated.

Fig. 1 .
Fig. 1.Overview of results for the "ideal case" run.Different parameters are plotted, at three different cuts through the simulation box.The top panel (x=0 cut/letters A, B) contains the plasma density and magnetic field just behind Rhea.The second panel (y=0 cut/letters C to F)shows the north-south configuration of the plasma density, velocity and electromagnetic fields, while the last panel of plots (z=0 cut/letters G to J) includes the same parameters, projected on Rhea's equatorial plane.In all panels, a small part of the simulation box is not shown, as it contains uninteresting regions (e.g.far upstream of Rhea).We also note that all vector parameters are plotted in the rest frame of Rhea.In several cases, perturbations in the vector quantities are small and the vectors seem unchanged compared to the background vector value and orientation (e.g.panel C).In this case, the perturbation is mainly ralized in the color coded vector magnitude.Additional plots in the plasma rest frame or with the background magnetic field subtracted, which reveal small perturbations better, are shown in Sect. 4.

Fig. 2 .
Fig. 2. Illustation of Eqs.(7).The position up to which the rarefaction wave has propagated is given by −(C s t+1).The dependence of the solution from the sound speed, tells us that the wake structure in the xz plane is determined by the plasma temperature.We note that in most previous studies the wake boundary is set at z=0 and the limb location at x limb =0, which results in a slightly different form of Eqs.(7).Since we want to in compliance with our simulation coordinate system, this slight modification in the equations has been introduced.

Fig. 3 .
Fig. 3. Comparison of Eq. (7) solutions with the simulated data, for a distance of 1.1 R Rh downstream of Rhea's limb.The scattering of the simulation data points is due to numerical noise, deriving from the small number of particles per simulation cell (about 15 particles per cell on average).

Fig. 4 .
Fig. 4. The asymmetric shape of Rhea's wake parallel and perpendicular to the magnetic field lines.The cut is at ∼2.4 R Rh downstream of Rhea (or at x=1.3 R Rh ).

Fig. 7 .
Fig. 7. Perturbations of v y and v z in an yz-cut just behind Rhea.The expansion of the plasma velocity inside the cavity is clearly seen.The peak amplitude of the v y or v z components is about 30 km s −1 .

Fig. 7 .
Fig. 7. Perturbations of v y and v z in an yz-cut just behind Rhea.The expansion of the plasma velocity inside the cavity is clearly seen.The peak amplitude of the v y or v z components is about 30 km s −1 .

Fig. 8 .
Fig.8.Perturbations of v x and v z in an yz-cut just behind Rhea, in the plasma rest frame.Plasma expands along the magnetic field lines, symmetrically with respect to the z=0 line.

Fig. 8 .
Fig.8.Perturbations of v x and v z in an yz-cut just behind Rhea, in the plasma rest frame.Plasma expands along the magnetic field lines, symmetrically with respect to the z=0 line.

Fig. 9 .
Fig. 9. Velocity and electric field in the plasma rest frame, projected on the equatorial plane (top and bottom panel respectively).Fig.9. Velocity and electric field in the plasma rest frame, projected on the equatorial plane (top and bottom panel, respectively).

Fig. 10 .Fig. 10 .
Fig. 10.Comparison of Equation 10 solution (blue line) with the simulated data, for a distance of 1.2 R Rh downstream of Rhea's limb.Triangles show the variation of vz from our simulation results, while diamonds show the velocity variation when the residual vx is considered.Fig. 10.Comparison of Eq. (10) solution (blue line) with the simulated data, for a distance of 1.2 R Rh downstream of Rhea's limb.Triangles show the variation of v z from our simulation results, while diamonds show the velocity variation when the residual v x is considered.

Fig. 11 .
Fig. 11.Comparison of the simulation results with the magnetometer data from the Cassini flyby.In green color we also show the magnetic field signature on an equatorial trajectory in the wake along the y-axis, from the "ideal case" run.Fig.11.Comparison of the simulation results with the magnetometer data from the Cassini flyby.In green color we also show the magnetic field signature on an equatorial trajectory in the wake along the y-axis, from the "ideal case" run.

Fig. 12 .
Fig.12.Illustration of the density (top) and velocity (bottom) of escaping heavy ions from a weak exosphere of Rhea.The plots shown are cuts through the equatorial plane of Rhea (xy-plane).The direction of the upstream plasma velocity and corotational electric field are shown on the top panel.Note that for this test run, Rhea was at the center of the simulation box.Here the heavy ion mass was chosen to be 17 amu and the exospheric scale height was set to 100 km.

Fig. 12 .
Fig.12.Illustration of the density (top) and velocity (bottom) of escaping heavy ions from a weak exosphere of Rhea.The plots shown are cuts through the equatorial plane of Rhea (xy-plane).The direction of the upstream plasma velocity and corotational electric field are shown on the top panel.Note that for this test run, Rhea was at the center of the simulation box.Here the heavy ion mass was chosen to be 17 amu and the exospheric scale height was set to 100 km.

Fig. 13 .
Fig. 13.Synthetic ion spectrogram along the Rhea flyby trajectory.The color represents the ion flux, normalized towards the maximum.

Fig. 13 .
Fig. 13.Synthetic ion spectrogram along the Rhea flyby trajectory.The color represents the ion flux, normalized towards the maximum.

Fig. 14 .
Fig. 14. (Top) Phase-space plot of Vz along a hypothetical trajectory which is parallel to the x-axis and crosses Rhea through its center.The white region in the plot denotes the space occupied by Rhea's volume.(Bottom) Sketch that explains the phase-space plot in the top panel.The dashed line shows the hypothetical trajectory along which the diagram of the top panel has been drawn.

Fig. 14 .
Fig. 14. (Top) Phase-space plot of Vz along a hypothetical trajectory which is parallel to the x-axis and crosses Rhea through its center.The white region in the plot denotes the space occupied by Rhea's volume.(Bottom) Sketch that explains the phase-space plot in the top panel.The dashed line shows the hypothetical trajectory along which the diagram of the top panel has been drawn.

Fig. 14 .
Fig. 14. (Top) Phase-space plot of V z along a hypothetical trajectory which is parallel to the x-axis and crosses Rhea through its center.The white region in the plot denotes the space occupied by Rhea's volume.(Bottom) Sketch that explains the phase-space plot in the top panel.The dashed line shows the hypothetical trajectory along which the diagram of the top panel has been drawn.

Table 1 .
Values and ranges of the various parameters describing Rhea's space environment.We note that the value range of plasma moments could partly be attributed to the moment calculation method.

Table 2 .
List of parameters used for the two simulation runs.