Surface current balance and thermoelectric whistler wings at airless astrophysical bodies: Cassini at Rhea.

Sharp magnetic perturbations found by the Cassini spacecraft at the edge of the Rhea flux tube are consistent with field-aligned flux tube currents. The current system results from the difference of ion and electron gyroradii and the requirement to balance currents on the sharp Rhea surface. Differential-type hybrid codes that solve for ion velocity and magnetic field have an intrinsic difficulty modeling the plasma absorber's sharp surface. We overcome this problem by instead using integral equations to solve for ion and electron currents and obtain agreement with the magnetic perturbations at Rhea's flux tube edge. An analysis of the plasma dispersion relations and Cassini data reveals that field-guided whistler waves initiated by (1) the electron velocity anisotropy in the flux tube and (2) interaction with surface sheath electrostatic waves on topographic scales may facilitate propagation of the current system to large distances from Rhea. Current systems like those at Rhea should occur generally, for plasma absorbers of any size such as spacecraft or planetary bodies, in a wide range of space plasma environments. Motion through the plasma is not essential since the current system is thermodynamic in origin, excited by heat flow into the object. The requirements are a difference of ion and electron gyroradii and a sharp surface, i.e., without a significant thick atmosphere. KEY POINTS: Surface current balance condition yields a current system at astronomical bodiesCurrent system possible for sharp (airless) objects of any sizeCurrent system is thermoelectric and motion through the plasma nonessential.

1. Introduction

Nonconductive objects in equilibrium with space plasmas are subject to the requirement that time-averaged positive and negative currents into (i.e., normal to) their surfaces, including ambient electrons and ions, and secondary and photoelectrons must balance at every point. A classic problem in spacecraft design and operation is determining the local conditions of this equilibrium, including local potential, sheath formation, and local surface charge [Garrett, 1981; Wipple, 1981]. Many studies have extended this physics to surfaces of planetary bodies [Roussos et al., 2010; Nordheim et al., 2014], e.g., to the problem of solar wind charging of the lunar surface [Farrell et al., 2010; Freeman and Ibrahim, 1975; Manka, 1973; Stubbs et al., 2014]. In this paper we consider the question of surface current balance in the presence of an ambient planetary magnetic field, as elucidated by Cassini spacecraft measurements of Saturn's satellite Rhea during two close encounters on 2 March 2010 and 11 January 2011.

As shown in Figure1, the flybys—designated R2 and R3, respectively—passed directly through the Rhea flux tube over the north and south poles, at 97 and 72 km altitude. As expected for a plasma-absorbing body, the Cassini fluxgate magnetometer (MAG) [Dougherty et al., 2004] measured (Figure2) an augmentation of the southward (−Z-directed) magnetic field due to the diamagnetic current circulating Rhea's flux tube and equator orthogonal to the field, which serves to compensate the absorption of plasma pressure by the body. Since Rhea orbits within Saturn's inner magnetosphere at 8.75 Saturn radii, it's orbital motion is overtaken by the corotating (collisionless) plasma, forming a leading side plasma wake which also carries a portion of the diamagnetic current (Figure1).

Figure 1

Rhea R2 and R3 north and south polar flyby configuration, showing the relation of the Moon, the flux tube, the ambient Saturn magnetic and corotation electric field directions, and the XYZ coordinate system. X and O denote fields and flows into and out of the page, respectively. Black dashed line: for reference, the R1 wake encounter. (a) North polar view showing (green dashed line) the expected diamagnetic current flow and the related magnetic field enhancement (pink X) which compensates pressure absorption by Rhea. -parallel ion diffusion causes rapid wake infilling downstream of Rhea, leading to diamagnetic current closure across the wake which may couple to Alfvénic field-aligned currents due to ion gyroradius effects (green X, O) [Simon et al., 2012]. (b) Equatorial view looking upstream, showing the observed -directed enhancement in the flux tube and orthogonal perturbations circling the flux tube edge (dotted lines).

Figure 2

Cassini MAG , , and measurements (black) versus spacecraft Y coordinate for the (top) R2 and (bottom) R3 encounters, with the differential (red, from Simon et al. [2012]) and integral (blue, this work) hybrid models also shown. The IHM assumes 4 (R2) and 7 (R3) cm−3 densities and a blurred Rhea resistivity profile to obtain approximate agreement with the DHM. The models reproduce the broad features of the data but not the sharp peaks aligned with the flux tube edge at Y ∼ ±1 (dashed lines).

However, MAG also measured intense perturbations of the magnetic fields in the X and Y directions (Figure2), i.e., orthogonal to the -directed ambient field, which is unusual. Such orthogonal perturbations are reminiscent of the field-draping pattern associated with the Pedersen current or pickup ions of a planetary atmosphere in relative motion to the ambient plasma, as seen, for example, at Enceladus [Simon et al., 2011], Io [Saur et al., 2002], and Europa [Volwerk et al., 2007]. In fact, Cassini's Ion Neutral Mass Spectrometer [Waite et al., 2004] found an O2 and CO2 atmosphere at Rhea during the R2 and R3 flybys [Teolis et al., 2010; Teolis and Waite, 2012], but the molecular abundance (∼2.5 × 1029 O2) yields conductivity values in the range 0.11–0.43 Ω−1, roughly 2 orders of magnitude less than necessary to explain the orthogonal field perturbations [Simon et al., 2012]. We have also analyzed (this work) Cassini Plasma Spectrometer (CAPS) [Young et al., 2004] data from the downstream R1 2005 and R1.5 2007 Rhea flybys and estimate Rhea's total pickup ion source to be ∼1022 ions/s (∼1 g/s) or ∼1.6 kA, still substantially less than the ∼14 kA (see below) required to explain the field perturbations. Progress was made by Simon et al. [2012], who showed that the short length (∼ 4 Rhea radii) of Rhea's plasma wake results in a downstream diamagnetic current, orthogonal to the flow and parallel to the ambient corotation electric field, which has sufficient magnitude to account for the strength of the orthogonal perturbations. They propose that the current and associated electromagnetic disturbance is propagated against the plasma flow by a standing current-carrying Alfvén wave, or “wing,” similar to the magnetospheric interaction of satellites possessing significant conductive atmospheres [Neubauer, 1998], and Rhea Alfvén wings have been confirmed by downstream MAG data from two distant (102 and 54 Rhea radii) Cassini flybys [Khurana et al., 2012]. Using a standard differential-type hybrid modeling (DHM) code (particle ions and magnetohydrodynamic (MHD) fluid electrons) that solves a differential equation for the electron fluid, Simon et al. [2012] have simulated the modification of Rhea's wake currents due to ion gyroradius effects and predict field perturbations generally resembling the measured signatures (Figure2). However, here we focus on the abrupt changes of , measured coincidently with the flux tube boundary on both flybys, which the DHM fit did not capture (Figure2). Each change occurred over short distances of roughly 100 km much shorter than the plasma's ∼500 km ion inertial length and therefore too short to be Alfvénic. The abrupt magnetic perturbation is an unequivocal indication of a field-aligned current at the flux tube boundary intersecting Cassini's trajectory. This flux tube current appears to be a consistent feature of Rhea's plasma interaction, as implied by its presence on two Cassini encounters in opposite hemispheres.

The concept of a field-aligned axially symmetric “wire” current was discussed by Santolik et al. [2011] and Simon et al. [2012], who noted that a current homogeneously distributed in the flux tube provides, at best, only a rough fit to the observed field perturbations. In Figure3 we have applied the axially symmetric form of Ampere's law in cylindrical coordinates to estimate the required (inhomogeneous) current distribution in the flux tube. The currents are of the order 10−9 A/m2, which is unfortunately below the ∼10−8 A/m2 discernible by CAPS, according to our analysis of the R2 and R3 ion and electron flux data acquired by the CAPS Ion Mass Spectrometer and Electron Spectrometer. The total flux tube current is ∼ 7 kA in each hemisphere, i.e., ∼14 kA in total assuming north/south symmetry. One can see that the field-aligned current, which flows out from Rhea, is concentrated at the flux tube edges as noted by Santolik et al. [2011], with a return “shielding” current apparently just outside the flux tube (Figure3).

Figure 3

(a) From R2 data (top): Ampere's law estimate of the -directed current versus radial distance from the flux tube (black line), compared with the IHM prediction for a sharp Rhea resistivity profile (red line), and the analytical expression (with for r >1) for the electron current needed to balance uniform irradiation of a sphere by a 270 eV, 106 m−3 ion plasma (blue dashed line). Note that the infinities of at r = ±1 are artifacts owing to the neglect of the finite electron gyroradius. MAG data (bottom) (black) and the expression (or for r >1) for the azimuthal wire field magnitude from (blue dashed line). (b) Same as Figure3a for R3. The broad hump in on R3 (not captured by the analytic expression) is a second-order effect due to the higher plasma (6 cm−3) than R2 (3 cm−3); see Figure8 and text. Both flybys show evidence for current peaks outflowing from Rhea near the edge of the flux tube (consistent with the analytical and IHM prediction) and an oppositely directed shielding current just outside the tube (not anticipated by the analytical expression but predicted by the IHM).

Figure 8

Same as Figure2 except that here we show the IHM results (blue) for a sharp Rhea resistivity profile. The R2 prediction closely resembles the wire model (Figure3a). The assumption of a higher density at R3 (6 cm−3) than R2 (3 cm−3), in accordance with RPWS data, reproduced the broad hump in on R3. The hump results from the second-order Alfvénic disturbance at scales larger than dsi (∼500 km), which is more significant at the higher R3 density. The IHM's rendering of ultrasharp features is limited by (i) binning resolution (currently 0.16R = 122 km) and possibly (ii) the steady state approximation (e.g., oblique whistler wavefronts near the flux tube edge are not modeled since they are temporally variable). Even with these constraints, the sharp Rhea assumption enables the IHM to successfully predict the magnitude and alignment of the peaks with the edge of the flux tube at Y ∼ ±1.

We consider here the concept sketched in Figure4, invoking (i) the difference of average ion and electron gyroradii rgi and rge (Table 1) and (ii) the requirement for current balance on the sharp Rhea surface to explain the observed current system. As shown the ions are absorbed into Rhea from (roughly) all directions, discharging across magnetic field lines with ease due to the large gyroradius, while the electrons are constrained to flow north/south along the field lines. Hence, the flux tube current is carried by the electrons, which flow most strongly along the tube rim, into Rhea's low-latitude surface to compensate ion deposition there. Simply considering the electron flow necessary to balance ion flux onto the sphere (assuming as an approximation uniform flux), we obtain to first-order for r ≤ 1 (r > 1), and a resultant -perturbation circular around the z axis of magnitude for r ≤ 1(r > 1), with radial coordinate r in units of Rhea radii R. As shown in Figure3, these simple expressions predict the basic shape and magnitude of the observed field perturbations and in particular the current and field perturbation maxima at the flux tube edges.

Figure 4

Schematic of the Rhea flux tube current system. Ions discharge into the body approximately uniformly (we neglect for simplicity the leading/trailing hemispherical difference due to corotation flow), while electrons are confined to flow along the field lines due to the smaller gyroradius. The requirement to balance ion and electron currents on the sharp surface yields a flux tube current flowing away from Rhea. The current is maximum on the flux tube edge due to the oblique angle of the magnetic field to the planetary equator.

Table 1

Field-Guided Wave Emissions for Flux Tube Thermoelectric Current Systems in Different Space Plasma Environmentsa

	rge	dse	rgi	dsi	References	Condition	Current-Carrying Wavesd,e
Solar wind (1 AU)	1.4	1.9	50	83	McComas et al. [2013] and Newbury et al. [1998]	dse < rgi ∼ dsi	whistler + partial Alfvén
Rhea R2 flyby	1.1	2.8	440	490	This work. Plasma parameters estimated from Cassini CAPS and MAG data.	dse < rgi ∼ dsi	whistler + partial Alfvén
Rhea R3 flyby	0.55	2.2	350	370	This work. Plasma parameters estimated from Cassini CAPS and MAG data.	dse < rgi ∼ dsi	whistler + partial Alfvén
Dione	0.084	1.2	55	200	Gustafsson and Wahlund [2010] and Wilson et al. [2008]	dse < rgi < dsi	whistler
Tethys	0.10	0.97	17	170	Simon et al. [2009]	dse < rgi < dsi	whistler
Jovian plasma (Io orbit)	0.0031	0.11	2.0	21	Kivelson et al. [2004]	dse < rgi < dsi	whistler
Jovian plasma (Callisto orbit)	2.7	13	160	2400	Kivelson et al. [2004]b	dse < rgi < dsi	whistler
Local interstellar space	4.6	19	200	810	Burlaga et al. [2013], Gurnett et al. [2013], and Slavin and Frisch [2002]	dse < rgi < dsi	whistler
Earth GS orbit (6.6 Re)	0.024	1.68	0.79	72	Keyser et al. [2009]	dse ∼ rgi < dsi	partial whistler
Earth low orbit (400 km)	0.00003	0.017	0.0034	2.9	Roble et al. [1978]c	δ < rgi < dse	none, currents carried entirely by flux tube anisotropy

There are four distinct regimes: δ < rgi < dse (current carried by anisotropy), dse < rgi < dsi (whistler), rge < dsi < rgi (whistler + Alfvén), r > dsi (Alfvén), though as shown some environments are borderline between regimes. Whistlers are the predominant wave mode in a wide range of environments. Distances in km.

Conditions highly variable along orbit. Using average plasma properties.

Satisfying δ < rgi, δ = 1 cm. Ion and electron mean free paths (∼0.1 and 10 km at 400 km) >> gyroradii [Kelley and Heelis, 1989], but whistler excitation nevertheless inhibited by the strong terrestrial field (average ∼35 μT) and cold ions (∼0.1 eV).

“Whistler + partial Alfvén” refers to rgi ∼ dsi cases where spreading whistler wings may weakly excite Alfvén wings at sufficient distances from the body.

Results are applicable to the flux tube current and do not preclude the possibility of secondary Alfvénic disturbances from the wake or flux tube.

Crucially, this current system only exists if the plasma absorber's surface is sharply defined, unlike, e.g., an absorber surrounded by a significant thick atmosphere. The DHM has numerical difficulty in the case of the sharp absorber, which can account for the problems in resolving the flux tube current. We discuss this issue below, and using a solution algorithm adapted to the case of a sharp surface, we demonstrate the agreement of self-consistent simulations with the current system described above.

2. Differential Hybrid Models: The Sharp Surface Problem

Differential-type hybrid models are discussed at length in literature [Lipat="jgra0119-8881-math-000ov, 2002; Müller et al., 2011], and thus, we shall only summarize the basic approach here. The DHM bins the rectangular simulation space into a hierarchical Cartesian grid and randomly initializes a statistically large sample of ions with (typically) Maxwellian-distributed speeds as individual particles upstream of the plasma-absorbing obstacle, advancing the ion positions and velocities in finite time steps Δt according to the equation of motion (neglecting gravity):

where , q, and m are the ion velocity, charge, and mass and and are the electric and magnetic fields at the position of ion i. Ions striking the body are removed from the simulation, and the ion/electron density n is estimated from the number of simulated ions in each grid cell, with quasi-neutrality n = n = n assumed. The steady state electric field satisfies the force balance condition

where ρ denotes the resistivity due to ion-electron collisions, the flow velocity of the massless electron fluid, the electron pressure tensor divergence, and the definition of current is used to express the condition in terms of . The DHM solves the magnetic diffusion equation

obtained from equation (2) using Faraday's law and Ampere's Law without the (assumed negligible) displacement current term, which ensures continuity . Typical boundary conditions are to fix to (Dirichlet condition) at the upstream face of the simulation box (with the ambient field measured far from the body), its spatial derivatives to zero (Neumann condition) across the downstream face, and periodic conditions at the side faces [Bagdonat, 2005]. is often approximated by the scalar pressure gradient [Roussos et al., 2008], but for anisotropic pressures P⊥ ≠ P∥ this generalizes (assuming a gyrotropic electron velocity distribution) to

containing magnetic mirror and curvature forces in the third and fourth terms [Hu and Denton, 2009]. The DHM discretizes the spatial derivatives with finite cell widths Δx and iterates equation (3) over sufficient finite time steps to converge the field to steady state, i.e., , setting and the (diverging) n− 1 term to zero within the body. The model is suitable for the study of the ion gyroradius effects at bodies, e.g., comets [Bagdonat, 2005; Müller et al., 2011], Enceladus [Kriegel et al., 2009], and Titan [Sillanpää, 2008], with conductive atmospheres possessing smoothly varying resistivity.

Airless objects like Rhea present a dilemma for the DHM since, as evident on expansion to , the resistivity term in equation (3) explodes due to (i) large ρ inside the body and (ii) the convergence of ∇ρ to a Dirac delta function on the body surface, as Δx → 0. The ∇ρ catastrophe at the body surface will, in fact, give rise to numerical artifacts including false current systems. The standard preventative measure is to simply remove the surface resistivity step by smoothing, e.g., through suppression of high spatial frequency components [Holmstrom, 2013], or imposition of an arbitrary function [Müller et al., 2011; Vernisse et al., 2013] that smears some resistivity into the space outside the body, with the required smoothing width a increasing with internal resistivity ρ0. The implicit assumption is that the smoothed surface is sufficiently sharp to accurately render the physics of the plasma interaction. Realistically, the smoothing imposes a possible limitation on the DHM's applicability [Vernisse et al., 2013, p. 39], since insufficient internal resistance allows too much current to permeate the body, while excessive ρ0 may inhibit exterior currents by requiring too much surface smoothing. The Δt and Δx used (which are limited by computing speed) may in some cases not allow for any physically realistic combination of ρ0 and a to be simulated. Other variations of the DHM encounter similar difficulties; e.g., Holmstrom et al. [2012] neglect the resistivity terms in their study of the lunar solar wind interaction but must still smooth n at the object surface to avoid numerical differentiation of the surface step. The issue is not solvable, e.g., by imposing a surface boundary condition , since the condition is already inherently “satisfied” by the DHM, in a rough sense, within the smeared “surface” [Simon et al., 2012], simply by way of the fact that in equation (3). Rather, the difficulty is numerical and lies in imposing a resistivity and/or density boundary onto the DHM with sufficient sharpness to model the relevant physics.

3. Our Solution: Integral Hybrid Model

Though differential expressions like equation (3) are often preferred due to computational efficiency, integral equations tend to be better behaved numerically and may sometimes be necessary in special cases, e.g., solutions involving large derivatives [Greengard and Rokhlin, 1991; Ledvina et al., 2008] as at Rhea. We have therefore implemented the physics of the electron fluid in the form of integral equations to treat Rhea's sharp surface—an approach we designate here as the integral hybrid model (IHM). We treat the ions as particles exactly as in the DHM but diverge from the DHM's derivation by using to include ion and electron current explicitly in the force balance condition:

where and . From equation (5) we solve for the -perpendicular electron current in terms of

where , A = enρ, , and the P and H subscripts signify the alignment of the components with the Pedersen and Hall directions in the case. Inside the body where , equations 6a, 6b, 6c yield , with K and n (in the context of the solid) the Hall coefficient and charge carrier density (i.e., n = − 1/eK), respectively. However, the Hall terms are insignificant if ρ is very large inside the body, i.e., , with in the infinite resistivity limit. In the other extreme case of A = enρ → 0 (i.e., low resistivity and/or density), equations 6a, 6b, 6c yield , i.e., , far from the body (where ion gyroradius effects on are negligible). With the parallel current components included (below), these conditions generalize to in the limits of (i) low density (few charge carriers), (ii) large resistivity (high ion/electron collision rates), and (iii), with some exceptions, low resistivity. Therefore, equations (6a), (6b), and (6c) anticipate, for example, that a thick planetary ionosphere with a gradient of ρ versus altitude would (as expected) conduct most current in the altitudinal “shell” of intermediate resistivity. The exceptions to (iii) which allow for , i.e., current flow even in the low-resistivity case, are (1) ion gyroradius effects acting on , (2) pressure gradients resulting in diamagnetic currents, and (3) sources or sinks along a field line which draw parallel electron current. Rhea's flux tube current system is an example of exception (3).

Whereas the DHM enforces continuity via Ampere's law, here we directly impose the condition , by requiring the parallel part of to satisfy , i.e., . We integrate this expression along the field lines using the identity , valid for integration on paths along field lines connected to , to solve for :

where , and and are the field line entry/exit points from the simulation space, assuming no closed field lines. The constant of integration M varies across field lines but is constant along a field line. This constant provides the unique for which the total parallel resistive force on the field line particles (i.e., the integral of from to ) is zero, as required in steady state. For the special case of zero field line resistance , M defaults to , resulting in equal and opposite parallel currents at each end of the field line. The role of M is to constrain the field line electron current “drawn” into or out of the simulation box to the (unique) solution which balances both currents and (parallel) forces. Here we assume equal field line resistance beyond and , i.e., in both directions outside the simulation box, such that the exterior regions make no net contribution to the force balance.

We compute the steady state electric field as the negative gradient of the potential − ∇φ:

where is the ambient corotation field and is the potential due to on the simulation box wall. Here is the field line integral of the parallel (Ohmic) electric field —given (from equation (5)) by . Finally, we evaluate from via the Biot-Savart law:

where is the ambient (Saturnian) magnetic field. Shielding effects prevent the simulated field in the simulation box interior from “seeing” the exterior currents, and therefore, equation (9) is well approximated by evaluating the volume integral only over the simulation domain (rather than all space). We minimize any errors near the boundaries by appending several external grid cell layers (typically eight layers, 1000 km total thickness) to the box faces. We continue (along the field lines) the currents and , as evaluated at the boundaries of the original box, into the new grid cells, and evaluate equation (9) in this extended box. These “pseudo” grid cells are only used for purposes of equation (9) and have no role in other steps of the IHM.

Writing equations 6a, 6b, 6c, 7a, 7b, 7c–9 in shorthand (excluding here fixed model inputs, e.g., ρ): , , , , with , we can construct the self-sufficient expression encapsulating all the physics in terms of the composite function ℋ:

which we implement computationally in the form of nested subroutines. Unlike the DHM which satisfies a condition (equation (3)) on , one can see in equation (10) that the IHM implements the same physics (Figure5) by satisfying a condition on , i.e., .

Figure 5

Computational flow of the electron fluid portion of the integral hybrid model. Asterisks: steps where relaxation parameters are applied; see text.

We bin the simulation space into a Cartesian grid as in the DHM and approximate the field line and volume integrals in ℋ by the midpoint Riemann summation method. The line integration algorithm applied to equations 7a, 7b, 7c and (8) iteratively retraces the field lines intersecting every grid cell as they evolve throughout the simulation. We converge the IHM to the self-consistent steady state by iterating a two-step sequence: (1) evaluate equation (10) repeatedly as an iterated function sequence until relaxes to an attractive fixed point then (2) update n and by executing one or more time steps of ion motion (equation (1)). Since ℋ is nonlinear, we use damping (relaxation) parameters (which weight the average of the current and prior iterations) of 0.5 and 0.2 on and (Figure5) to establish convergence in step (1). There is flexibility in the number of ion time steps in (2): we find that propagating the ions for a gyroperiod (t = 2πm/eB0 = 48 s), i.e., sufficient for the ions to adjust to the updated fields from (1), tends to minimize the IHM convergence computation time. Hence, with a time step of 0.75 s, sufficiently small to resolve the gyromotion, we run each iteration of (2) for 64 time steps. After achieving ±10% precision on the magnitudes of , , , , and n, we begin gradually increasing the number of time steps to 100, while averaging results with previous steps to reduce statistical uncertainty in and n. The simulation is finally terminated upon achieving ±1% precision. In accordance with equation (10), , n, and fully describe the system state in the model, with , , and obtainable at any time via equations 7a, 7b, 7c–9.

Contrary to the DHM, the electron fluid portion of the IHM makes no attempt to evaluate the temporal evolution of as in equation (3). Accordingly, by solving directly for the steady state, the IHM bypasses another inherent DHM limitation which occurs in two of the cases mentioned above: (i) high resistivity (e.g., an insulator such as Rhea) and (ii) low density (e.g., vacuum). Both cases lead to a divergence of terms (i.e., diverging magnetic diffusivity) on the right-hand side of equation (3). The difficulty has its origins in the DHM's neglect of the displacement current, which (by Maxwell's equations) determines the field propagation versus time in regions of low current, i.e., . Numerical tricks are required by the DHM in these limits, such as setting n−1 to zero at low density or assuming an arbitrary high (but not too high) resistivity in vacuum regions or within the insulating body [Holmstrom, 2013]. By contrast, the IHM has no difficulty in the low-current (i.e., high diffusivity) regions, such as Rhea's interior, since the steady state magnetic fields are readily given at any point in space by the “action at a distance” principal embodied in equation (9). We note also that the steady state displacement current is by definition zero and hence of no relevance to the time-independent IHM electron fluid algorithm.

Another benefit of integral equations is that the known constants of integration (M, , and in the IHM) are sufficient to establish a unique solution, thereby taking over the role of boundary conditions in a differential scheme. In this respect the IHM is more realistic, circumventing the requirement of the differential approach for assumptions, inherently somewhat arbitrary, about, e.g., , , and (in the DHM's case) at the simulation box walls. Rather, the values everywhere in the domain are already given uniquely by the integral equations 7a, 7b, 7c–9 as constrained by M, , and , respectively. Considering as an example the upstream boundary, the IHM gives the unique realistic solution for , i.e., very close to but spatially variable, rather than the exact equality typically enforced by the DHM.

While not required for Rhea, we note that the IHM is readily adaptable to bodies possessing a more significant neutral exosphere. Together with an assumed exospheric density profile, the required modifications are (i) an extra electron-neutral collisional resistivity term in equation (5) and (ii) inclusion of ion-neutral collisions and a pickup ion source.

The IHM maintains numerical stability in the sharp surface limit for several reasons. First, the integration over ρ in equation (8) mitigates the differentiation of the surface step by ∇. Similar reasoning applies to the discontinuity of and due to ion and electron deposition on the body surface, since the divergent derivatives, i.e., (for r = 1) as Δ x → 0, are safely reintegrated by equations 7a, 7b, 7c. Finally, while both the DHM and IHM exhibit a surface spike if , the issue is trivial in the IHM since the parallel component of the spike is integrated by equation 8. The perpendicular component contributes to a surface spike of , but the effect is integrated out by equation (9). These considerations are demonstrative of the advantageous numerical properties of integral equations in systems containing sharp discontinuities.

4. IHM Results at Rhea: Flux Tube Current System

We applied the IHM to Rhea by setting the average ion mass to 16 amu [Wilson et al., 2010], the ambient (upstream) ion velocity distribution to a bidirectional Maxwellian with T = 65 eV and T = 225 eV [Wilson et al., 2010], the Rhea-referenced corotation speed v0, and to 50 km/s and 22 nT southward [Khurana et al., 2008; Thomsen et al., 2010; Wilson et al., 2008], respectively, and to anti-Saturnward. We simulated a less dense plasma at R2 (n = 3 cm−3) than R3 (n = 6 cm−3), in accordance with the Cassini radio and plasma wave instrument (RPWS) [Gurnett et al., 2004] density measurements during these flybys [Roussos et al., 2012]. We model Rhea as a sharply defined sphere of infinite ρ, which forces in the body. We typically add a small amount (104Ωm) of artificial ambient resistivity to ensure numerical stability as is also commonly done in differential-type codes [Ledvina et al., 2008]. The simulation space was a 9R-sided cube, binned (uniformly) into 56 × 56 × 56 cells of width Δx = (9/56)R. We calculated from equation (4) using P∥ = nT and P⊥ = nT and estimated the average ambient R2(R3) electron temperatures T ∼ 30(8)eV and T ∼ 60(16) eV from the CAPS ELS data, setting P∥ = P⊥ = 0 inside Rhea. Additionally, we tracked B field lines intersecting Rhea throughout the simulation and reduced T and T by 70% on these field lines to include the effect of high-energy electron absorption seen by ELS [Jones et al., 2011]. As expected, we find that the resulting electron pressure loss enhances in the flux tube but does not contribute to the and components (which exhibit a similar response even if constant is everywhere assumed).

In Figure6 we show the IHM plasma density n and and streamlines on the Y = 0 (XZ) and equatorial Z = 0 (XY) planes. Ion absorption yields a thin density cavity surrounding Rhea and a downstream wake which expands more rapidly in Z than Y due to unrestrained -parallel ion diffusion. Contrary to the equatorial electron streamlines which flow smoothly around Rhea, Figure6 reveals that the equatorial ion flow is more responsive to the obstacle, exhibiting more cross-field curvature into Rhea and the wake than the electrons due to the larger ion gyroradius. This behavior is reversed in XZ to satisfy on Rhea's surface, with the -parallel electron flow more responsive to the obstacle than the ions (Figure6), converging more sharply into Rhea and the wake to compensate the cross-field ion convergence.

Figure 6

IHM Rhea results showing plasma density XY and XZ cross sections (blue-to-green: increasing density, with low-density wake on right) and the ion (black solid streamlines) and electron currents (red dashed streamlines, shown). Note that the orthogonal cross-section planes bisect Rhea, and therefore, only the northern anti-Saturn quadrant of the sphere is visible in the figure. The simulation assumes a sharp Rhea resistivity profile. Ions curve more strongly into Rhea on the equatorial plane due to the large gyroradius. In response, electrons curve more strongly into Rhea in the vertical XZ plane along the magnetic field to maintain surface current balance.

The sum of and yields the expected flux tube current system, with a current converging toward Rhea's equator and a commensurate outgoing current concentrated on the tube rim. The visualization of Figure7 shows how the field-aligned current flows north/south away from Rhea, roughly compensating the convergence of perpendicular current in the equatorial plane (see “ring” of negative around Rhea's equator in Figure7). The cylindrical shell wire generates a field circulating the flux tube, with the IHM correctly predicting the position (and the magnitude) of the peaks in on the flux tube edges during R2 and R3 (Figure8). Accordingly, the total field ( plus perturbation) twists about the flux tube, and the field-aligned current is therefore somewhat helical. A shielding current of opposite helicity, directed toward Rhea, also flows just outside the flux tube, not unlike that implied by Ampere's law from the MAG data in Figure3. Such flux rope topologies are characteristic of electron MHD disturbances, which exist on scales sharper than the ion inertial length dsi = c/ωpi (∼500 km or 0.65R at Rhea), and tend to propagate in tandem with whistler mode disturbances [Stenzel et al., 1999] as we discuss in the next section. The wire field also introduces a large-scale second-order perturbation to the ions, on MHD size scales (i.e., larger than dsi). The effect is insignificant at R2 but more substantial at R3 due to the higher plasma density, yielding a broad “hump” in the R3 component as also seen in the data (Figure8). As can be seen in Figure8 the DHM, as expected, has little difficulty capturing the large-scale (Alfvénic) [Simon et al., 2012] hump. However, the IHM goes a step farther, by also succeeding to render the sharp small-scale features at the flux tube edge (Figure8) which are indicative of current balance on the sharp Rhea surface.

Figure 7

IHM prediction of the Rhea current system showing the flux tube current flowing out from Rhea and oppositely directed shielding currents outside the flux tube. Red lines: northward (positive Z) flowing currents. Blue lines: southward (negative Z) currents. Equatorial cross section shows the diamagnetic current spiraling into Rhea (black streamlines) and the divergence of the perpendicular current component (color). One can see that the outgoing flux tube current compensates the ion convergence at Rhea's surface (red ring around Rhea). The simulation assumes a sharp Rhea resistivity profile.

The IHM gives a good approximation of the DHM fit (Figure2) after smoothing Rhea's surface resistivity step using (as an example) a ρ(r) = ρ0/(exp[(1 − r)/w] + 1) profile with ρ0 = 6 × 105Ωm and a smoothing width w = 0.2 Rhea radii. In this scenario (Figure9) the equatorial diamagnetic current is more circular around Rhea, with the convergent component mostly cut off due to the extension of resistivity outside the body, enabling ions to frictionally drag electrons across field lines into the surface. By breaking the circuit, the external resistivity layer automatically balances the surface current and suppresses the flux tube current system. A dichotomy in equatorial develops on the Saturn-facing and anti-Saturn-facing flanks of Rhea as seen in Figure9, yielding ingoing and outgoing field-aligned currents on these flanks, respectively. Unlike the sharp resistivity case (see next section), this DHM current system requires relative motion between the object and the plasma and depends on object size (due to the ion gyroradius effect) and plasma β as noted by Simon et al. [2012]. For example, the current system vanishes in simulation runs where we assume no bulk plasma flow past Rhea and diminishes when we reduce β, either by initializing a colder ion population or increasing the magnetic field. We concur in part with the explanation of Simon et al. [2012] in terms of cross-flow diamagnetic currents in the infilling downstream wake but would add that the variation of ion flow speed around the obstacle (a gyroradius effect) [see Roussos et al., 2008, Figure1-I] also appears to play a role in the simulation. The fast ions concentrated on Rhea's Saturn-facing flank couple to the electron fluid via the resistivity smeared outside the body, thereby forcing electrons anti-Saturnward into Rhea through induction and generating current divergence (positive ) and a corresponding ingoing field-aligned current (outgoing electrons) on this flank (Figure9). By contrast, the anti-Saturn convergence (negative ) results from deceleration of ions entering the wake across the field lines [Roussos et al., 2008, Figure1-I], which produces an outgoing field-aligned current (ingoing electrons) from this flank.

Figure 9

Same as Figure7 assuming a blurred Rhea resistivity profile in the IHM. The divergence (blue region) and convergence (red region) of the equatorial currents on the Saturn and anti-Saturn flanks of Rhea are compensated by a dichotomous current flowing vertically to/from Rhea. This current system (unlike Figure7) increases in magnitude with ion gyroradius and plasma β and requires relative motion between the object and the plasma. Moreover, the diamagnetic current is more circular around the body, i.e., prevented from spiraling directly into Rhea's surface by the extension of resistivity outside the body. As expected, the overall current system resembles that obtained in the DHM due to the smooth resistivity assumption.

Outside the flux tube the IHM and DHM yield similar results irrespective of the assumed sharpness of the body. For example, as shown in Figure10 for the 26 November 2005 R1 Cassini wake flyby (Figure1), both the IHM and DHM succeed to fit the MAG data, showing the expected intensification of in the wake. Roussos et al. [2008] give many detailed plots of the estimated plasma moments and fields around Rhea; their simulations agree with our IHM results outside the flux tube.

Figure 10

Fit to the R1 wake flyby (Figure1), with the DHM (red, from Roussos et al. [2008]) and IHM results (blue) both showing the expected increase of in the wake (see Roussos et al. [2008] for a detailed analysis of the shapes in all components). The IHM assumes a sharp Rhea, with = 25 nT and n = 6 cm−3. Dotted lines: Y ∼ ±1 radii. This example demonstrates the general agreement of the DHM and IHM on the plasma properties and fields outside the flux tube.

5. Current System Wave Modes and Physical Requirements

In the simplest terms Rhea's flux tube current can be understood as a direct consequence of the surface potential, which charges negative to equalize the ion and electron fluxes everywhere on the surface [Roussos et al., 2010]. Accordingly, the electron current distribution in the flux tube (Figure4) results from precipitation of electrons with sufficient energy to overcome the surface potential. However, we can gain additional perspective into the physics of Rhea's plasma interaction by also considering the wave modes excited by the body, their contribution to the currents and magnetic topology around Rhea, and the implications of the sharp surface. A basic dilemma is that (due to the surface step) the flux tube perturbations are too short to be MHD, implying that the classic Alfvén wing model traditionally applied to planetary satellites possessing significant exospheres/ionospheres [Neubauer, 1998] cannot give a complete description of Rhea's flux tube current system.

We start by comparing, in Figures1 and 2, the field perturbation produced by the sharp object and that of the smooth resistivity approximation. One can see in the smooth resistivity case (Figure11) a typical field-draping pattern, as would be expected for the Figure9 dichotomous current system discussed above. We expect a smoothly resistive object to excite modes with wavelengths of similar size to the object [Neubauer, 1998], which for bodies like Rhea larger than the ion inertial length yields magnetosonic and shear Alfvén waves characteristic of an MHD plasma. Only the Alfvén wing persists to large distances due to (ideally) perfect -parallel constructive interference of the Alfvén waves [Gurnett and Bhattacharjee, 2005]. One can see that the draped field pattern of Figure11 tilts at roughly the Alfvén angle θ = atan(v0/v) = 23°, with v0 = 40 km/s and v = 89 km/s the plasma flow and Alfvén speeds, respectively.

Figure 11

Blue streamlines: the IHM prediction of the magnetic field perturbation (i.e., total field minus constant background) in the smooth resistivity (R2, n = 3 cm−3) case. Equatorial cross section shows the plasma density (blue-to-red: increasing density) and the diamagnetic current circling Rhea and the wake (black streamlines). The perturbation is entirely southward directed as expected to compensate pressure loss around Rhea and the wake. One can see the field draping in the flow direction, as expected for an Alfvénic current system.

This result stands in stark contrast to the sharp surface case shown in Figure12, which exhibits (in addition to some field draping in the wake) the prominent flux tube helical perturbation discussed above.

Figure 12

Same as Figure11 for the sharp resistivity case, with southward and northward perturbations colored blue and red. The helical perturbation is due to the flux tube electron current. Contrary to the Alfvén disturbance from the wake, the flux tube feature has negligible draping angle to the ambient magnetic field, consistent with the high whistler group velocities of order 103 km/s.

Unlike the smoothly resistive object, we anticipate that the sharp object may also excite (in addition to MHD waves) short wavelength modes as the plasma responds to accommodate the surface step, with the near-parallel field alignment of the flux tube perturbation consistent with propagation of most of the excited wave energy along much faster than the Alfvén speed. On plotting the different branches of the plasma dispersion relation out to large wave numbers (Figure13), we find that the average wave numbers of the sharp flux tube perturbations (∼100 km, i.e., centered at ∼ 6 × 10−5 m−1) are well above the range of MHD waves. The spatial scale of the surface plasma absorption interaction is given by the topography (∼ 1 km) [Nimmo et al., 2010] on large scale and the Debye length δ on small scale (∼0.05 km) and therefore tends to excite modes in the ∼ 6–300 × 10−4 m−1 wave number range. The effect of plasma absorption at the sharp surface is twofold: (1) short wavelength (i.e., ∼0.05–1 km) electrostatic wave excitation from the surface sheath and (2) anisotropization of the flux tube electron velocity distribution onto which the surface step is “imprinted” as a sharp flux tube boundary (e.g., as exhibited in the CAPS ELS spectra [Jones et al., 2011; Santolik et al., 2011]). Here we briefly address the question of sheath waves, before discussing the generation of waves further up the flux tube by anisotropy-driven instabilities.

Figure 13

(top) Dispersion relations with Rhea plasma parameters approximated by cold plasma equations [Gurnett and Bhattacharjee, 2005; Hollweg, 1999; Stringer, 1963]. (bottom) Group velocity versus wave number averaged over all wave vector angles θ. Vertical dashed lines denote wave numbers for relevant length scales. Blue dashed line: limiting for simulation box size. Black dashed: average for the perturbations on the flux tube edge (Figure8). Olive dashed: surface topography and electron gyroradius scale, which roughly define the lower limiting for surface electrostatic wave excitation. Note that the surface excites wave numbers well above the MHD Alfvén range. Negatively sloped straight lines (Figure13, top): inverse ω versus relationship of the limiting positive (lower red) and negative (upper black) surface sheath charging timescales. Cross-hatched zone: approximate spatiotemporal scale of possible surface wave excitation, lying near or between the charging timescales and topographical and Debye spatial scale, and near several slow-propagating resonances. Nonlinear coupling via anisotropy-driven instabilities (i) at the surface and (ii) farther up the flux tube causes energy flow (curved arrows) into faster-propagating modes (Figure13, bottom) which travel large distances from Rhea. The asterisk denotes the field-aligned whistler beam (Figure14).

Figure 14

(right) Group velocity (GV) hodographs with direction (arrows) and magnitude (color) versus wave number (logarithmic radial coordinate) and propagation angle to , for the three low-frequency branches of the Rhea ambient plasma corresponding to the fast (bottom), slow (middle), and shear Alfvén (top) modes in the MHD limit. Blue arrows: GV directions in the plasma rest frame. Red hash marks: GV directions in the Rhea frame. Color: GV magnitude in the Rhea frame, shown in linear scale to demonstrate the negligibility of the GV in all modes except kinetic Alfvén (top, oblique angles only) and whistler (bottom, colored stripe spanning all angles). (left) GV angular distributions (white-to-black: low-to-high intensity, angle relative to in the Rhea frame) assuming isotropic wave vector excitation. Surface waves excited at topographical and Debye scales are far beyond the range of the Alfvén and slow shocks (top, middle) but can nonlinearly couple to the fast-propagating whistlers (bottom). Note for whistlers that the GV arrows align nearly to at all propagation angles (bottom; see preponderance of vertical arrows in the whistler band), forming a fast-propagating (mean GV of ∼6400 km/s) field-aligned beam near 10−3.56 m−1 with a small ∼1° spread. The measured perturbation (black dashed line) is centered at ∼10−4.2 m−1 in the whistler regime, i.e., well beyond MHD spatial scales. Inset: RPWS 17:41:26.39 UTC R2 flyby measurement (solid line) of intense of whistler frequency excitation in the Rhea flux tube from Santolik et al. [2011], compared to predicted spectra for the 10−3.56 m−1 whistler beam (black dashed line: assumes cosine distribution of wave vector propagation angles; red dotted line: isotropic wave vectors).

Sheath electrostatic waves may (i) acquire free energy from the electron/ion velocity distributions, which are highly anisotropic at the surface, and (ii) exchange thermal energy with the random component of the particle velocities [Gould, 2002]. The interaction frequencies are roughly speaking bounded by the limiting charging timescales, i.e., ion(electron) flux ratioed to the unit charge column density σ = ε0T/e2λ (i.e., surface charge plus sheath space charge) to repel the particles (where the σ expression treats the plasma penetration distance of the surface electric field as equivalent to the wavelength λ = 2π/k of the surface-coupled electrostatic wave). As shown in Figure13 with limiting frequencies ω = 2π/τcharging = 2πF/σ = C/k (where ), the relevant spatiotemporal scales correspond closely with (i) electron plasma (Langmuir) and upper hybrid waves near the plasma frequency (fep∼15 kHz), (ii) ion plasma waves (fip∼80 Hz), and (iii) electron cyclotron (fce∼650 Hz) waves (and, in a hot plasma, Bernstein modes), i.e., all slowly propagating resonance modes. The RPWS indeed detected intensification of all of these resonances in Rhea's flux tube [Santolik et al., 2011], though we note that the electron plasma resonance undergoes additional anisotropy-driven amplification farther up the flux tube [Santolik et al., 2011]. As the slowly propagating resonant waves accumulate energy near the sheath, some wave energy flows (Figure13) into fast-propagating off-resonance frequencies through nonlinear interactions. The principal nonlinear interaction pathway is by way of electrostatic wave-particle interactions in the sheath, which contribute to the electron velocity anisotropy by scattering thermal electrons into the surface [Iess et al., 1998]. Additionally, a major contribution to the anisotropy also comes from fast ambient electrons arriving with sufficient energy to overcome the sheath potential without significantly interacting with the sheath waves. The anisotropy and associated electron current is propagated up the flux tube [Santolik et al., 2011], where thermodynamic free energy is continuously exchanged between the electron velocity distribution and different wave modes (and, likewise, between waves modes via the distribution).

The current system dissipates (i.e., closes) gradually up the flux tube as free energy escapes into fast-propagating wave modes not guided along the field lines, including high-frequency (>15 kHz) electron acoustic waves (∼4400 km/s) shorter than the ∼0.05 km Debye length and low-frequency (<80 Hz) ion acoustic waves (∼40 km/s) above kilometer wavelengths. RPWS, in fact, detected enhanced broadband excitation at Rhea both above fpe (presumably electron acoustic and/or doppler-shifted plasma waves) in association with “bursty” high-frequency emissions and at frequencies below ∼40 Hz in the ion acoustic range [Santolik et al., 2011]. By contrast, free energy exchange with field-guided wave modes may facilitate the survival of the flux tube current system to large distances from Rhea. Current-carrying kinetic Alfvén waves do have high field-aligned group velocities of order 103 km/s, but energy flow into this mode is likely minor since the waves only occur in a small sliver of phase space at oblique wave vector angles above ∼87° (Figure14).

However, whistlers at wavelengths above the electron inertial length dse = c/ωpe ∼3 km are an obvious candidate to perpetuate the currents since they (i) are well focused along and (ii) have the highest average group speed (peaking at ∼8000 km/s) of any mode on the dispersion plot (Figure13, bottom) spanning all wave vector angles (Figure14). Physically, the whistlers are the continuation of the MHD fast waves to wavelengths shorter than dsi for which electron motion dominates the plasma wave response, which is consistent with our findings that the flux tube current is carried by the electrons. Intense whistler activity inside the Rhea flux tube was confirmed by RPWS (Figure14, inset), likely deriving energy from the electron cyclotron-resonant instability [Santolik et al., 2011]. The whistlers may be initiated in part at the surface itself and farther up the field lines due to cyclotron-resonant (parallel energy ∼230 eV) electrons reflected from negative surfaces [Jones et al., 2011]. Similar processes are known to operate at the Moon, where whistler emission is stimulated by solar wind electron interactions with crustal magnetic fields [Halekas et al., 2006] and electron reflection from negatively charged surfaces [Halekas et al., 2012] and lunar wake ambipolar electric fields [Farrell et al., 1996; Nakagawa et al., 2003]. Comparing the group velocity field of the three low-frequency dispersion surfaces, versus wave number and vector angle θ to (Figure14), we can see that only the whistler group velocities exhibit alignment parallel to at all θ in Rhea's reference frame. The narrow group velocity cone angle of ∼1° for near ∼ 10− 3.56m− 1 (Figure14, bottom left), or ∼23 km wavelength, is within the ∼100 km thickness of the measured field perturbations (Figure2). The near-parallel whistler beam alignment is due to the high average group speed, ∼6400 km/s, i.e., much faster than the Alfvén and slow shocks (89 and 37 km/s) which exist at much lower wave numbers below ∼10−6 m−1 (i.e., wavelengths > dsi) and at substantial angles to (23 and 44°) in Rhea's reference frame. For this reason the whistler Poynting flux is field aligned [Santolik et al., 2011].

Slow-propagating oblique wavefronts at the flux tube edge may contribute to the sharpness of the magnetic perturbations sampled by MAG, by analogy, e.g., with oblique whistler wave fields near perpendicular shocks [Hellinger et al., 2007]. For instance, the flux tube edge may contain wave components with temporal variability on small scale not considered by the (steady state) IHM. The lack of these wave components in the IHM may explain the model's underestimation of the steepness of the perturbations (Figure8). The phase speed approaches ∼ v = 89 km/s [Stringer, 1963] (in the plasma rest frame) for perpendicular whistler propagation at wavelengths the size of the flux tube perturbation (∼100 km), and the wave normals can pass through 90° to (no resonance cone) since the frequencies (f ≃ v/λ = 0.9 Hz) are below the lower hybrid frequency (fLH = 4.2 Hz) [Gurnett and Bhattacharjee, 2005]. At shorter wavelengths (f → fLH) the perpendicular phase velocity slows further as the resonance cone regime is approached, with the waves becoming electrostatic lower hybrid oscillations [Bell and Ngo, 1990]. However, wavefronts not oriented precisely perpendicular to also have significant transverse electron currents and magnetic degrees of freedom [Bellan, 2013], and such waves directed near radially to the flux tube may contribute to the shielding current and fields near the tube rim (Figure3). Wavefronts drifting out of phase with the sharp flux tube boundary (in the electron current and velocity distribution [Santolik et al., 2011]) may be rapidly Landau damped [Zhang et al., 1993], thereby transferring energy back to the electrons [Cattell et al., 2012], analogous, for instance, to the damping of oblique precursor whistler waves observed drifting ahead of sharp planetary bow shocks [Gary and Mellott, 1985; Orlowski et al., 1995]. Current-carrying “whistler wings” consisting of oblique wavefronts [Stenzel, 1999] have also been demonstrated in laboratory experiments in which the motion of a conducting electrode through a magnetized plasma was simulated by current pulse superposition [Stenzel and Urrutia, 1990, 1993], and Cerenkov whistler emissions are known to form from the motion of small (< dsi) magnetized asteroids through the solar wind [Baumgartel et al., 1997; Gurnett, 1995; Omidi et al., 2002; Simon et al., 2006; Wang and Kivelson, 1996].

A noteworthy property of the flux tube current system is that motion of the body through the plasma is not required as for a (smooth) conductor, with the currents retained even in IHM runs with the plasma at rest relative to Rhea. The reason is that the current system is thermodynamically, rather than inductively driven, being excited solely through plasma absorption, i.e., by heat flow from the plasma into the body. The currents are therefore proportional to the ion implantation flux into the surface, i.e., to the density times the ion thermal speed, as verified by our IHM runs. A related example of “thermal” field-aligned plasma currents may be the “short-circuit” effect first discussed by Simon [1955] in magnetically confined laboratory plasmas, in which cross-field thermal ion diffusion, e.g., to the axial walls of a conductive containment vessel [Zhilinskii and Tsendin, 1980], can induce field-aligned electron eddy currents which short across field lines at the end walls [Drentje et al., 2002]. The competition of these currents with electron ambipolar diffusion has been the focus of recent theoretical discussion [Fruchtman, 2009] and simulation studies [Lafleur and Boswell, 2012], and thermally driven electron shorting currents turn out to be a significant design consideration for electron cyclotron ion sources [Schachter et al., 2008]. Analogous phenomena are common in other areas of physics as well. An interesting example is a thermoelectric circuit, with the difference of Seebeck coefficients between two materials enabling temperature gradient-driven currents (via the material's different densities of electronic states), playing the same role as the difference of rgi and rge for a plasma-absorbing body (via the different phase space densities of ion/electron trajectories intersecting the sharp absorber).

The limiting lateral spatial scales of the current system about the flux tube are defined by rgi and rge. The current system may propagate in tandem with whistler waves as long as the spatial extent of the whistler characteristics are within the rgi-rge size range. Current closure begins at the distance up the field lines for which the whistler beam spreads to become larger than rgi (assuming rgi > rge), provided that rgi < dsi, with dsi the minimum spatial scale for MHD waves. Taking 1° group velocity cone angle as a lower limit for Rhea's whistlers (Figure13), we find that the Rhea whistler wing may extend at most to 30,000 km up the field lines before current closure begins, i.e., outside our simulation box, but well short of the ∼1.6 × 106 km distance to Saturn. However, at Rhea's L shell rgi and dsi are similar (Table 1), and therefore, the current system and associated magnetic anomaly may spread sufficiently to significantly perturb the ion trajectories prior to current closure, thereby weakly exciting Alfvén wings which can propagate currents much further. Evidence for branching of currents between fast field-aligned electron “beams” and slower-propagating Alfvén waves has also been found at Io, with the observations of multiple Jovian auroral spots in association with Io's Alfvén wings [Bonfond et al., 2008]. Whistler wing spreading north/south of Rhea might also scatter energetic electrons away from Rhea, possibly explaining [Santolik et al., 2011] the puzzling energetic electron depletion measured near Rhea on several of the Cassini flybys [Roussos et al., 2012]. The depletion was at first interpreted as possible evidence of a Rhea ring system [Jones et al., 2007], but nondetection of rings by Cassini imaging appeared to rule out this explanation [Tiscareno et al., 2010].

Strong magnetic fields and/or cold or tenuous plasma (for which rgi is shorter than the electron inertial length dse = c/ωpe) are at the limit of the whistler wing description, since whistlers transition into electron cyclotron oscillations for wavelengths below dse due to electron inertia. In this limit the current system closes at a distance up the flux tube sufficient for the dissipation of the electron velocity anisotropy, which may occur as free energy flows from the velocity distribution into nonfield-guided (short wavelength) waves, such as Langmuir waves via Landau damping. At more extreme (and unlikely) high field intensities and low densities for which rgi < δ, the current system is completely inhibited since the surface is, from the plasma perspective, only as sharp as its Debye sheath. In the opposite (also unlikely) limit, i.e., very weak fields and/or hot or dense plasma for which rge > dsi, the spatial scales are entirely in the MHD range, allowing for an entirely Alfvénic current system. These criteria assume collisionless plasma: with collisions the current system propagation distance is also limited by the ion and electron mean free paths. Note that the conditions for the existence of the current system depend only on plasma parameters, and not object size (Figure15), e.g., a spacecraft or a planet may generate a current, so long as the surface is sharp (e.g., a airless body, or one with an atmosphere much thinner than rgi). Simply by absorbing plasma from its environment, even the Cassini spacecraft should excite a weak whistler-generating current system with a total estimated current of ∼50 nA, and minimum spatial scale (given by rge) about a kilometer at Rhea, i.e., much larger than the spacecraft. As shown in Table 1 the conditions for thermoelectric whistler wing formation are satisfied in a wide range of space plasma environments.

Figure 15

Schematic of the flux tube current system for a small (e.g., spacecraft), medium, and large plasma absorber (relative to the ion gyroradius). The current system exists for absorbers of any size, provided there is a difference of electron and ion gyroradii. We note that even large bodies on MHD size scales can excite the currents due to the surface step, which contains frequency components above the MHD range. (top) Ion (green) and electron (black) flows. (bottom) Total current (red).

6. Conclusions

Cassini's findings at Rhea yield clues on the fundamental plasma process of surface current balance on a sharp-absorbing body and raise the fascinating question of how this balance is achieved in a magnetic field for which ions and electrons approach the body on different trajectories. In this paper we have demonstrated from first principles, and with self-consistent modeling, that the difference of ion and electron average gyroradii yields a current system in the flux tube with unique properties: (1) motion through the plasma is not required since the current is produced thermoelectrically through work done by heat flow into the object, (2) the current system can form with objects of any size (e.g., the relation of the object size to the ion gyroradius is not critical), and (3) an object with a sharp surface (much sharper than rgi) is necessary, i.e., without a significant thick atmosphere. We demonstrate that standard modeling approaches employing differential equations encounter numerical problems with the sharp surface, and we therefore implement a so-called integral hybrid model to solve the physics with an integral equation approach better suited to cope with the surface step. We suggest that this type of approach can serve as a starting point for future treatments of a new class of modeling problems involving sharp-surface plasma absorbers, such as airless planetary bodies and spacecraft, immersed in magnetized space plasmas. Finally, we show on the basis of the plasma dispersion relations and Cassini RPWS results at Rhea that whistler waves stimulated by anisotropy-driven instabilities in the flux tube and interaction with surface sheath electrostatic waves on topographic scales may facilitate transmission of the current system up the flux tube.

The discovery of a Rhea flux tube current system addressed in this paper, and that of “secondary” Alfvén wings [Khurana et al., 2012] from the wake diamagnetic [Simon et al., 2012] and flux tube currents (Figure8, at R3), paint a new and surprisingly complex picture of the general interaction of objects with space plasmas beyond the scope of cold plasma theory. Four fundamental current systems can now be distinguished (the first two familiar and the last two new): (1) Alfvén wings from a conductive atmosphere in relative motion to the plasma, (2) the diamagnetic current (Figure1a), (3) the secondary Alfvén wings, and (4) the flux tube current. Only current system (1)—that generated by a conductive atmosphere—exists in the cold plasma limit but requires motion through the plasma (the canonical pure induction case). The remaining three current systems require a warm plasma, and (2) and (4)—the diamagnetic and flux tube currents, respectively—can form in a stationary plasma (the pure thermoelectric case). In a moving plasma (2) and (4) can give rise to secondary Alfvén wings, i.e., current system (3). Hence, (3) has a hybrid induction-thermoelectric character, with kinetic effects enabling pathways (specifically, current systems (2) and (4)) for momentum exchange (via Alfvén wings) between the object and the moving plasma. A significant difference between (2) and (4) is the inverse B0 dependence of the diamagnetic current (2), implying, e.g., a negligible diamagnetic contribution to Alfvén wings in low β plasmas, as is common (see, e.g., the case of Dione and Tethys discussed by Simon et al. [2012]). By contrast, the flux tube current (4) is proportional only to ion density and thermal speed as discussed above and is not quenched by strong magnetic fields alone (except in the extreme rgi < δ limit mentioned above). Flux tube current systems are therefore expected to form under quite general conditions, to satisfy surface current balance, e.g., at the other Saturnian (or other planetary) satellites, or objects in the solar wind such as asteroids and Earth's Moon. This potential broad implication of Cassini's results at Rhea starkly illustrates the necessity for new modeling approaches as discussed here and for future in situ measurements at solar system bodies to further elucidate the fundamentals of the plasma interaction and, in particular, the role of thermoelectric current systems and surface current balance.