Induced magnetic moments from a nearly spherical ocean.

The five largest planets all have strong intrinsic magnetic fields that interact with their satellites, many of which contain electrically conducting materials on global scales. Conducting bodies exposed to time-varying magnetic fields induce secondary magnetic fields from movement of eddy currents. In the case of spherically symmetric conducting bodies, matching magnetic solutions at the boundary results in relatively simple relations between the excitation field and the induced field. In this work, we determine the more complicated induced magnetic field from a near-spherical conductor, where the outer boundary is expanded in spherical harmonics. Under the approximations that the excitation field is uniform at a single frequency, the product of wavenumber and radius for the body is large, and the average radius of the body is large compared to the perturbation from spherical symmetry, we find that each spherical harmonic in the shape expansion induces discrete magnetic moments that are independent from the other harmonics in the expansion. That is, simple superposition applies to the magnetic moments induced by each perturbation harmonic. We present a table of the magnetic moments induced by each spherical harmonic up to degree 2 in the perturbed shape. We also present a simple formula by which the induced magnetic field may be evaluated for any arbitrary shape described by expanding the radius of the conducting body in spherical harmonics. Unlike the Earth, many moons in the Solar System are tidally locked to their parent bodies, and many also contain saline, subsurface oceans. Conductive material in these moons is therefore expected to be non-spherical. Accounting for the boundary shape of Europa's ocean will be critical for interpretation of Europa Clipper magnetic measurements near the moon, where the effects of quadrupole-and-higher magnetic moments will be most apparent. The results of this work permit magnetic studies considering non-spherical oceans of satellites for the first time.

Introduction

Induced magnetic fields result when a conducting body is subjected to a time-varying electromagnetic field. The time-varying magnetic field generates electric fields in accordance with Maxwell’s laws, which in turn drive currents in the conducting body, producing a secondary magnetic field. Because currents can only flow where there is conducting material, the shape of the conducting body affects the form of the resulting induced field.

The induced magnetic field for the configuration of a spherically symmetric conductor has been described by Lahiri and Price (1939) and Srivastava (1966), and summarized in application to planetary bodies by Parkinson (1983). Each of these authors examine cases of finite conductivity that vary only with radius. A significant body of work exists in the study of non-uniform conductivity structures within the Earth, but much of it is of limited application to other Solar System bodies. Some studies focus on inverting measurements to find induced magnetic moments for the Earth, e.g. Egbert and Kelbert (2012). Numerous authors consider forward modeling of a heterogeneously conducting Earth, especially in designing numerical solutions intended to aid in magnetic sounding (Fainberg and Zinger, 1980; Zhang and Schultz, 1992; Velímský and Martinec, 2005; Weiss, 2010). All geomagnetic studies benefit from the wealth of magnetic measurements taken over long periods of time in, on, and near the Earth. These extensive reservoirs of data help inform the variety of models that may be applied and make numerical simulations for this purpose far more practical.

In contrast, magnetic measurements of most Solar System bodies are very sparse. Earth’s moon is the only moon in the Solar System that, as yet, has had any spacecraft orbiting it. Past geocentric studies relying on a wealth of training data are of limited use in application to bodies with data-limited magnetic surveys—there is a need for techniques capable of constraining global-scale conductivity structures with sparse magnetic measurements. No previous work has evaluated the effect on induced magnetic moments from perturbing the outer surface of the conducting body. Crucially, as we will demonstrate in Section 3.3, deviations from spherical symmetry can alter the induced dipole moments, and result in induced magnetic moments of higher order than that of the excitation field, even when the excitation field is uniform. Although the signal from magnetic moments of greater than dipole order is expected to be small for existing spacecraft measurements (e.g., those from the Galileo mission), future missions are more likely to be able to resolve higher-order magnetic moments. For example, the Europa Clipper mission plans numerous flybys that will approach within 25 km of Europa’s surface (Campagnola et al., 2019; Lam et al., 2018), less than 2% of its radius in altitude.

Understanding the magnetic fields induced from nearly spherical conductors is of particular importance in the study of moon–magnetosphere interactions. The five largest planets all have global-scale magnetic fields, and all of them except Saturn have their dipole moments misaligned with their rotation axes. This generates oscillations in their magnetospheres at rotational, synodic, and orbital periods. All of the outer planets also have large moons whose self-gravity pulls them into nearly spherical shapes, and all large moons in the Solar System are tidally locked to their parent planets. Thus, they all experience asymmetric gravitational fields that perturb their shapes. Many of these moons also possess salty, global oceans of liquid water that are excellent conductors of electricity (Khurana et al., 2002) and promising candidate environments in the search for life elsewhere (Vance et al., 2018). There are, therefore, many near-spherical ocean worlds important to the study of astrobiology that are subjected to periodic magnetic fields. In fact, magnetic sounding has provided the strongest positive confirmation of an ocean of liquid water underneath Europa’s icy crust (Schilling et al., 2007; Hand and Chyba, 2007). Nevertheless, all past studies of magnetic sounding for planetary bodies have assumed spherical symmetry in modeling subsurface oceans, thereby modeling only induced dipole moments (Parkinson, 1983). Although this is a sensible first-order approximation, the Galileo flyby MAG data for Europa are not particularly well-fit by assuming a spherical ocean (cf. Schilling et al. (2007), Figs. 5–9). This leaves ample room for ambiguity that may be improved upon by considering non-spherical oceans, especially with the current scarcity of satellite magnetic measurements—only five flyby datasets have been widely studied for magnetic induction of Europa.

In most cases, little is known about planetary interiors. Interpretation of spacecraft magnetic measurements requires a careful accounting of the plasma environment and ionospheric effects in addition to planetary fields (Schilling et al., 2007). However, there remains a need for techniques capable of constraining the interior structure. As we will show in Section 3.3, asymmetry in subsurface oceans yields induced magnetic moments of greater than dipole order; signals from these will be much more pronounced during the close flybys planned for the Europa Clipper mission. Interpretation of measurements from these close flybys will require models capable of predicting higher-order magnetic moments resulting from interior structure (Soderlund, 2019; Vance et al., 2020), from plasma sources (Rubin et al., 2015) and from putative plumes (Arnold et al., 2019; Jia et al., 2018; Roth et al., 2014). Failing to account for degeneracies in the possible sources of higher-order magnetic moments will invalidate the results of such investigations.

There are strong indications that inhomogeneity in the ice shell of Europa will contribute to large-scale variance in thermal transport through the crust. In particular, cryovolcanism (Quick and Marsh, 2016) and upwelling of diapirs (Tobie et al., 2003) are expected to sustain latitudinal variation in the thickness of Europa’s ice shell. A compensatory variation in the outer radius of Europa’s ocean will influence the flow of electric currents generated by the excitation field applied by Jupiter. Careful examination of magnetic measurements near Europa may be able to yield information about the shape of the ice–ocean boundary, but this would be impossible without an appropriate model for the induced magnetic field.

It is important to note that the expectation of non-spherical conducting bodies within Solar System moons is in contrast to the Earth. Smaller size, lower bulk densities, and proximity to parent planets apply proportionally stronger gravitational asymmetries. Patterns of tidal heating and heat escape may further distort the conductor outer boundary away from a spherical shape. In addition, strong oscillatory magnetic fields from the rapidly rotating gas giants combine with these effects to result in larger signals from non-spherical conducting bodies in satellites compared to induced fields from the Earth.

In this work, we evaluate the magnetic moments induced from a spherical conducting body whose outer surface is perturbed in terms of spherical harmonics. The induced moments are expressed in terms of the harmonic coefficients of the excitation field. Most smoothly varying, near-spherical shapes may be approximated by retaining only low-order harmonics. In Section 3.3, we prove that expanding to first order in our perturbation parameter leads to simple superposition from shape harmonics to resultant magnetic moments—each harmonic appearing in the perturbed conductor’s shape generates proportional magnetic moments that are independent of the other harmonics. This means that any shape whose departure from spherical symmetry may be approximated by a linear combination of spherical harmonics of degree n = 1, 2 has induced magnetic moments that are trivial to calculate with the results of our method, using Table 1. This treatment enables, for the first time, consideration of subsurface oceans that are not spherically symmetric. Applying this work to the particular case of Europa, examined through the lens of magnetic measurements from the Galileo mission, is done in a companion work (Styczinski and Harnett, 2020, in prep.).

Table 1

Results of multiplying derivatives of normalized spherical harmonics by shape harmonics S, calculated in terms of derivatives of other normalized spherical harmonics. The functions Y(θ, ϕ) and S(θ, ϕ) are tabulated in Appendix A. Values for Γ may be read directly from the table, and are color-coded based on n′ of the resulting harmonic (online only). n, m index the column, p, q match the subscripts of the shape functions S, and n′, m′ identify the term within the corresponding table cell. Entries are calculated by multiplying together the functions found in Appendix A, then solving for the proportionality constant of the highest order resulting harmonic, then using that factor to calculate the next highest order resulting harmonic, and so on. Although the selection rules are nearly identical, the values for Γ are not the same for and S. However, because the radial expansion terms vanish, we only tabulate the values for the transverse terms here. In this work, we tabulate only p = 1, 2, as magnetic moments of higher than octupole order are expected to be undetectable at spacecraft orbital distances for Solar System moons.

Shape	∂Y_1−1∂θ	∂Y_10∂θ	∂Y_11∂θ
S 1–1	1245∂Y_2−2∂θ	1225∂Y_2−1∂θ−32Y_00e^−iϕ	12215∂Y_20∂θ
S 10	1215∂Y_2−1∂θ+1232Y_00e^−iϕ	12415∂Y_20∂θ	1215∂Y_21∂θ−1232Y_00e^iϕ
S 11	12215∂Y_20∂θ	1225∂Y_21∂θ+32Y_00e^iϕ	1245∂Y_22∂θ
S 2–2	132435∂Y_3−3∂θ	13835∂Y_3−2∂θ−83Y_1−1e^−iϕ	138175∂Y_3−1∂θ−415∂Y_1−1∂θ
S 2–1	131635∂Y_3−2∂θ+23Y_1−1e^−iϕ	1364175∂Y_3−1∂θ−43825∂Y_1−1∂θ	1324175∂Y_30∂θ+12825∂Y_10∂θ
S 20	1318175∂Y_3−1∂θ+35∂Y_1−1∂θ	1327175∂Y_30∂θ−15∂Y_10∂θ	1318175∂Y_31∂θ+35∂Y_11∂θ
S 21	1324175∂Y_30∂θ+12825∂Y_10∂θ	1364175∂Y_31∂θ−43825∂Y_11∂θ	131635∂Y_32∂θ−23Y_11e^iϕ
S 22	138175∂Y_31∂θ−415∂Y_11∂θ	13835∂Y_32∂θ+83Y_11e^iϕ	132435∂Y_33∂θ

Methods

We are concerned with the induced magnetic field for a configuration analogous to that of global-scale oceans of large moons in our solar system. The primary magnetic field applied to these bodies is dipolar at typical orbital distances, and neglecting contributions from magnetospheric plasma oscillations, is effectively uniform across the body (Zimmer et al., 2000). For simplicity, we therefore restrict our focus to an excitation field that is uniform across the conducting body. Similar methods may be applied to excitation fields with greater spatial variation, but the relatively simple case of a uniform field is sufficient for application to moons. We apply our method to Europa in detail in a companion article (Styczinski and Harnett, 2020, in prep.) to constrain the asymmetry that may be present in its subsurface ocean.

The excitation field applied to the conducting body takes the form with static and dynamic components that are complex in general. The measurable magnetic field is found by taking the real part of any complex expressions. Coefficients in any expression may be complex, except where noted. Although Eq. (1) accounts for the possibility of multiple oscillation frequencies applied to the conducting body, we will restrict our approach to a single frequency. Superposition permits an independent handling of each excitation frequency, so our method may be applied repeatedly for a combination of frequencies.

In the absence of electric currents, the magnetic field satisfies Laplace’s equation and therefore may be described by the gradient of a scalar potential ψ:

These equations are valid outside the conducting body if we neglect currents in the plasma. Although the currents present in the plasma environment around the body are not generally negligible, the principle of superposition permits us to consider each contribution to the net electromagnetic response independently—the net magnetic field is then the sum from each individual contribution. In this work, we consider only the induced magnetic moments generated by the interaction of the primary excitation field with the conducting body, as this is the dominant interaction that induces magnetic fields from within solar system moons.

The near-spherical shape of the conducting body encourages our use of spherical polar coordinates r = (r, θ, ϕ). Solutions to Eqs. (2) and (3) are easily represented in these coordinates by the following general form (Jackson, 1999): where B and B are complex coefficients for the external and induced magnetic fields, Y are the spherical harmonics, and a is a unit of radial distance. Orthonormal expressions for spherical harmonics up to n = 3 are listed in Appendix A. Magnetic potentials proportional to positive powers of r can only be generated from outside of the conducting body under examination; because we are assuming a uniform excitation field, we keep only the n = 1 terms for B. The B terms in Eq. (4) are those of the multipole expansion, so each B is proportional to, and thus represents, an induced multipole moment. We assume the magnetic potential for the excitation field oscillates sinusoidally, so time dependence is added to Eq. (4) by multiplication of e− as in Eq. (1).

Within the conducting body, the dynamic excitation field induces electric fields that drive currents, so B cannot be represented by Eq. (3) in this region. Instead, we must use a diffusion equation for B, derived from combining Maxwell’s laws with Ohm’s law:

For simplicity in deriving our model, we neglect movement of conducting material within the body, which can itself induce secondary fields (Saur et al., 2010; Vance et al., 2020). Considering only the oscillatory magnetic field, taking the time derivative of B is equivalent to multiplication by −iω. We can thus rewrite Eq. (5) in terms of a diffusion constant k, and arrive at a vector Helmholtz equation:

In general, μ and σ are functions of position and will vary throughout the body. On planetary scales, μ is well approximated by μ, even for bodies containing large amounts of ferromagnetic materials (Saur et al., 2010). We further assume that σ is uniform within the conducting body. Although this may not be valid on global scales, and the conductivity will increase with depth in a well-mixed ocean (Vance et al., 2018), we will later assume that the average outer radius of the conducting body a is large. When the imaginary part of ka is large, the radial dependence of the magnetic field within the conducting body becomes an exponential function whose strength decreases with depth. When the thickness of the conducting layer τ is significantly larger than the skin depth , deeper conductivity structure will be well-screened from the oscillating excitation field, and ceases to affect the induced field (Neubauer, 1999). When this condition is satisfied, any material below the uppermost conducting layer may be safely ignored. For simplicity and clarity in our derivations, we also make the approximation that any material outside the primary conducting body has negligible conductivity. Validity of these approximations is discussed in Sections 2.3 and 4, and a diagram of the interior structure model we apply is depicted in Fig. 1.

Fig. 1.

Near-spherical 3-layer model applied to calculate the induced magnetic field for a near-spherical conducting ocean. Conductivities σ and radii r are labeled for each region, as are the body radius R, average ocean depth τ, and the conducting body outer radius r(θ, ϕ). For the example case of Europa, an expected iron core and rocky mantle are depicted. Layers are not to scale; r(θ, ϕ) is discussed in detail in Section 3.2.

General solutions to Eq. (6) for the configuration at hand must be consistent with a poloidal field; since they are induced by an external field, there will be no toroidal field component (Moffatt, 1978). Poloidal fields take the following form: where the poloidal potential P is a scalar function of position. The quantity in square brackets in Eq. (9) is the angular momentum operator , of which the spherical harmonics Y are eigenfunctions:

If we suppose P is separable, we can expand it in spherical harmonics: where A are constant coefficients determined by the boundary conditions and are functions we must determine from other relations. As we later satisfy the boundary conditions with this functional form of P, the uniqueness theorem confirms that this is the physically correct representation, validating the supposition that P is separable.

Inserting Eq. (13) into Eqs. (9)-(11) and utilizing Eq. (12) yields expressions for the components of the magnetic field within the conducting body in terms of :

We can now make use of these expressions along with Eq. (6) to find a differential equation for . Linearity of the ∇2 operator allows us to consider only a single n, m term, as the same equations will apply to all terms. The component of Eq. (6) reads as

Inserting Eqs. (14)-(16) and again exploiting the angular momentum operator, we arrive at a Bessel equation for :

Solutions to this equation are linear combinations of the particular solutions, which are the spherical Bessel functions of the first and second kind, j and y:

C and D in Eq. (19) are constants determined from the boundary conditions. The solutions y diverge at the origin, so D must always be zero for the innermost region. As our solution for now contains arbitrary coefficients C and D, we may absorb the coefficients A and rename C and D to C and D.

For later convenience, we also require expressions for and . Eqs. (20) and (21) can be manipulated to obtain which we now also define as and , respectively.

We can now write general expressions for the magnetic field in all regions. From Eqs. (3) and (4), the external magnetic field follows where δ1, is the Kronecker delta function, which reflects our consideration of a uniform excitation field. Inserting Eq. (19) into Eqs. (14)-(16), and absorbing A into C and D, the internal magnetic field follows

Y is an eigenfunction of with an eigenvalue of im, so we later replace these derivatives. The operator is not conveniently replaced; we evaluate it as needed when we examine a near-spherical boundary in Section 3.3.

Solving for the B, C, and D coefficients in these equations is accomplished by applying Maxwell’s laws at the common boundaries between each region. On each boundary surface, Maxwell’s laws dictate that the normal component of B must be continuous, and the tangential components of B/μ must be continuous whenever there are no surface currents confined to the boundary itself. As we assume μ = μ within the body, B is continuous everywhere, and the components of the magnetic field for adjacent regions are equal on the boundary surface.

Internal boundary conditions

Boundaries interior to the outer surface must match component values according to Eqs. (27)-(29) at each boundary surface. In this work, we assume that interior boundary shapes are spherically symmetric. At a radius r describing the outer boundary of a lower layer l with wavenumber k and an upper boundary u with k, the internal boundary conditions read

Because the interior boundary surfaces are spherically symmetric, the C and D coefficients are independent of m.

As first introduced by Srivastava (1966), a recursion relation may be obtained from the internal boundary conditions. Factoring out the C coefficients and expressing the ratio D as Λ, we can solve for in terms of by dividing Eq. (30) by Eq. (31). For a given interior structure model, j and y are functions of known quantities k and r, so is completely determined by the structure model and . In the innermost layer, which contains the origin, Λ are all zero so that the resultant magnetic field remains finite. The equations above may therefore be repeatedly applied to find Λ for each successive layer, all the way to the outer boundary. This is a necessary step to reducing the number of variables in the exterior boundary conditions to a solvable set—the outermost Λ is completely determined from the interior structure of the conducting body.

External boundary conditions

Combining Eqs. (24)-(29), we obtain the boundary conditions that apply at the outermost boundary for a conductor of any shape: where we have omitted the i/ sin θ factor common to all B terms. The shape of the outer boundary is defined by a function r(θ, ϕ), which we insert in place of r in the above boundary condition equations. The transverse boundary conditions (Eqs. (33) and (34)) provide redundant information regarding the B and C coefficients.

Every term in the boundary condition equations contains a spherical harmonic; these functions are all mutually independent. From the boundary condition equations, we can therefore obtain a set of equations that are linear in the B and C coefficients. These linear equations are sufficient to determine all coefficients analytically if we make a pair of approximations:

∣kr∣ → ∞, which is reasonable for very high conductivity, for very large conductors of moderate to high conductivity, or for rapidly varying magnetic fields (Parkinson, 1983). As a corollary, the skin depth s becomes small in this regime, suggesting the interior becomes screened from oscillating fields.

The boundary surface may be adequately represented with a few low-order spherical harmonics.

The first approximation above is discussed in detail in Section 2.3. The second approximation is necessary because, as we will see in Section 3.3, all harmonics in the boundary shape adjust either the dipole moments or quadrupole moments relative to the spherically symmetric case. However, these adjustments are proportional to the amount of representation a given harmonic has in the boundary shape, so harmonics that have comparatively small coefficients result in smaller adjustments to the induced magnetic moments. We must therefore assume the series expansion of the boundary may be truncated, i.e. the boundary shape is dominated by low-order harmonics.

Finally, the induced magnetic field outside the conducting body is found by taking the real part of

B will always be expressed as linear combinations of B. The complex values B are determined, for planetary bodies, by measuring the excitation field from the parent planet. Although the values measured for the excitation field are purely real, the complex coefficients B are determined by expressing the oscillatory components of the field in terms of complex spherical harmonics Y.

High-conductivity approximation

Europa has an outer radius of 1560 km and Jupiter has a synodic rotation period of 11.1 h with Europa. Supposing Europa has a nominal ice shell thickness of 30 km, this approximation requires that Europa’s ocean has a conductivity of order 1 S/m or greater. For a conductivity comparable to that of Earth’s ocean at 2.75 S/m (Hand and Chyba, 2007), ∣kr∣ is approximately 35; next-to-leading-order terms are about a factor of 20 smaller than the leading-order term. Measurements from the Galileo mission confirm that the orientation and strength of the induced field is close to that of a perfect conductor (Zimmer et al., 2000), which justifies the use of this approximation. We therefore retain only leading-order terms in 1/kr when applying this approximation.

The radial functions take the following asymptotic forms (Heald and Marion, 2012):

With these expressions, at the bottom of the ocean Eqs. (30) and (31) become independent of the lower layers and yield a result for Λ:

The second term in Eq. (38) results in a negligible change to the imaginary part of Λ, but contributes a real part of order . Although this is small for large ∣ka∣, it multiplies an exponential that grows as in both and , so it is important to retain as part of applying the high-conductivity approximation. We can also see from Eq. (38) that the real part of Λ scales as , causing deeper oceans to converge more closely to an exponential form exp {−ikr}.

For the approximate values above for Europa, the ocean skin depth s is about 60 km. Gravity measurements from the Galileo mission indicate that the ocean depth is likely of order 100 km (Anderson et al., 1998), somewhat larger than . As a consequence, the Bessel functions will be dominated by the growing exponentials in Eqs. (36) and (37). We can then replace the Bessel function expressions found in the external boundary conditions with simpler forms: where we have absorbed constant factors into the final coefficient . We will also later need which may be obtained using the l’Hôpital rule and recursion relations for Bessel functions.

Results

First, we examine the case of spherical symmetry to prove consistency with this well-established solution. Next, we apply our method to an arbitrary near-spherical boundary shape. Finally, we summarize our results with a simple formula for calculating induced moments for an arbitrary near-spherical boundary.

Induced magnetic field for a spherical conductor

A boundary surface with spherical symmetry is defined by r(θ, ϕ) = a. Inserting this into the boundary condition equations yields an exact solution with our approximations—only n = 1 terms have non-zero B and C coefficients. Collecting terms proportional to each Y (or ) yields equations of the form where indicates that this is the zeroth-order solution for B. The B and B equations offer redundant information. With the expressions from Section 2.3 we apply our kr → ∞ limit. Keeping only first order terms in 1/kr, we get:

Solving for in Eq. (44) and inserting into Eq. (43), we gain a term that is reduced by a factor 1/ka relative to the other terms. Discarding this term, we can conclude

Because the excitation field is uniform, containing only n = 1 harmonics, the induced magnetic field is a pure dipole for this boundary shape.

Near-spherical boundary shape

We must now define the surface r(θ, ϕ) for a near-spherical boundary that we will insert into the boundary conditions. Because we wish to represent physical surfaces, r(θ, ϕ) must be purely real. Furthermore, we endeavor to represent subsurface oceans of Solar System moons; to do so, we must limit the deviation from spherical symmetry. This ensures that the surface described by r(θ, ϕ) remains at or beneath the planetary surface, when a is the average outer radius of the conducting ocean. With these considerations in mind, we choose a surface of the form where ε, which is purely real, is the maximum deviation from spherical symmetry, χ are purely real or purely imaginary constants, and S are proportional to spherical harmonics of the same indices.

The harmonic functions S are chosen such that they represent a functional dependence in accordance with the matching normalized harmonics, but have proportionality constants chosen such that they have a maximum range ∣S∣ ≤ 1. Shape harmonics in our chosen normalization are listed in Appendix A.

The χ coefficients must be purely real or purely imaginary because we require that linear combinations of S be purely real to describe a physical surface. χ represent the relative amount of each harmonic present in the boundary shape. Typically ∑∣x∣ = 1; this is not strictly true because the range of S20 is not symmetric about zero. Coefficients χ must be chosen so that the maximum value of r(θ, ϕ) is a + ε, or the minimum value is a – ε, in accordance with the physical shape to be approximated. We will later show in Section 3.3 that, under our approximations, the magnetic moments for any arbitrary choice of χ may be readily calculated from Table 1.

Surfaces described by Eq. (46) are near-spherical in that we make the approximation ε ≪ a. Equivalently, ε/a ≪ 1, and we retain terms up to first order in ε/a only. This enables us to use a Taylor expansion in the boundary conditions that truncates quickly, adding only one term containing a product of two spherical harmonics. Products of spherical harmonics may be expressed as a linear combination of different harmonics (Condon and Shortley, 1951). This results in “mixing” of harmonics in the excitation field from n = 1 into other n, so a uniform excitation field induces magnetic moments of quadrupole order or higher for this shape, in addition to the original dipole moments. Table 1 contains the results of multiplying each n = 1 spherical harmonic with each shape harmonic, up to p = 2. These are the products we will need in order to find linear equations for the resultant spherical harmonics that appear in the boundary conditions.

Induced magnetic moments for a near-spherical boundary shape

Let us now insert our near-spherical r(θ, ϕ) into the boundary conditions. To first order, a Taylor expansion of a function f(r) about r(θ, ϕ) has terms

The r power series that multiply B and B then take the form a(1 + nε/a ∑ χ). The interior field terms become

Applying these expansions to the exterior boundary condition equations (Eqs. (32)-(34)) results in new terms added to the n, m series:

The terms are identical to the terms with replaced by , so for brevity we focus only on a single transverse equation. Because we assume the boundary surface is near-spherical, there will always be at least one induced moment of degree n = 1 that is non-zero, with m matching the excitation field. Higher-order magnetic moments will be proportional to ε/a, so the only terms that will be non-negligible are those where B in Eqs. (50) and (51) have n = 1; all other terms will be of second order or higher in ε/a and we discard them.

We can further simplify the terms with the help of our zeroth-order solution. Using Eqs. (41) and (42) to replace C, and using Eq. (40) to make the relevant substitution for the high-∣kr∣ limit we get

The products will be sums of of degree and order where N = 0, 1, 2, … The resulting products are given by where Γ are constant coefficients, presented in Table 1 for n = 1, up to p = 2. The subscripts n, m represent the degree and order of the ‘input’ harmonic Y; p, q represent a shape ‘operator’ mixing input harmonics to outputs according to the particular shape of the boundary; n′, m′ represent the degree and order of the ‘output’ harmonics Y.

As terms for all n, m are added to the boundary condition equations, we can manipulate the subscripts on to more easily reduce the equations later. It is convenient to now swap the n, m and n′, m′ subscripts in Eq. (55) so that the output harmonics match the indices for the other terms in the n, m series. The subscripts on can refer to input or output harmonics, so we keep them n, m as before. With these adjustments, Eq. (53) results in

This operation is equivalent to collecting all series terms proportional to . We now define a quantity Ξ that captures the effect on the induced magnetic moment B based on the boundary shape for a given excitation field:

We can now evaluate how the near-spherical boundary shape will impact the induced magnetic moments. The orthogonality of the spherical harmonics may be exploited to obtain a set of 2 × 2 systems of equations that contain only the unknowns C and B for particular values of n and m. These orthogonality conditions are detailed in Appendix A. For the B conditions, we now need to express in terms of other harmonics:

We can now reduce the boundary condition equations. Multiplying both sides of the B boundary condition equations (Eq. (32)) by and integrating over a unit sphere is equivalent to replacing Y by δ. Similarly, multiplying both sides of the B boundary condition equations (Eq. (33)) by and integrating over a unit sphere is equivalent to replacing by . Because every term that survives after applying the delta functions is proportional to this quantity, which can never be zero, it divides away. In addition, all Y± terms in Table 1 vanish from this operation due to the orthogonality conditions. After these reductions, the external boundary conditions finally become the following set of linear equations:

Solving Eq. (60) for C and inserting into Eq. (59), we can solve directly for each induced moment. There is ultimately only a small change from the spherically symmetric solution in Eq. (45):

It is important to note that some B will contain terms with excitation moments of differing m′—for example, B may contain a term proportional to B, as is the case for non-zero χ21. This is true for all boundary surface harmonics of non-zero q. It is immediately apparent that taking ε → 0 in Eq. (61), equivalent to reverting to a spherically symmetric boundary surface, recovers the solution for that configuration (Eq. (45)).

Discussion

Here we discuss some key notes regarding our results, elaborate on how they may be applied, and detail the limits of their application.

Superposition of shape harmonics

Eq. (61) is linear in all χ and Γ, which is a powerful result. Values for χ, which remain attached to terms stemming from shapes S of the same p, q, are all mutually independent. This tells us that the changes to the induced magnetic moments driven by the non-spherical boundary are unique to each boundary harmonic. In other words, the moments induced by each harmonic in the near-spherical boundary shape obey simple superposition. Provided our χ yield real values for r(θ, ϕ) with a maximum value of a + ε, they are entirely arbitrary. We can therefore choose any boundary shape as a sum of harmonic functions as in Eq. (46) and immediately read values from Table 1 to write the magnetic moments induced by this shape.

Application of results

Although Eq. (61) appears complicated, its use is made rather simple by considering the entries in Table 1. Eq. (61) essentially states that induced moments of n, m differ from the spherically symmetric case by 3ε/2a times each value from the table that multiplies , times B, with m′ determined by the column. For example, consider the B magnetic moment. The spherically symmetric case has B = B/2, so we obtain our result by adding to that. There are three entries in Table 1 that multiply a Y11 harmonic: , , and . The second subscript for each is m′, so we multiply the first of these by B, the second by B, and the third by B. Therefore, the B induced magnetic moment for the near-spherical case is and the induced magnetic field for this moment is found by taking the real part of Eq. (35) with this B moment inserted.

We can obtain the same result by examining Eqs. (61) and (54):

We are seeking the contributions to the n = 1, m = 1 term.

n′ = 1, from our consideration of a uniform field, so the selection rules also tell us that p can be 2, 4, 6, … We assume that p = 3 and higher harmonics may be neglected in this work, so p = 2.

m must equal q + m′, from the selection rules (recall that Eq. (61) has swapped the primed quantities from Eq. (54)). Therefore, there are contributions to this magnetic moment only for shapes S with q = 0, 1, 2.

Finally, we multiply B by the coefficients for the table entries where p = 2, q + m′ = 1 and insert this product in place of ΓB in Eq. (61). As before, the table cells are those for , , and .

Our method can readily be applied to near-spherical shapes described by a function r(θ, ϕ) = a + ε ∑ with p continued to arbitrarily large degree. The values for Γ in Table 1 are simply derived from algorithmic algebra. Entries are calculated by multiplying together the functions in Appendix A, then solving for the proportionality constant of the highest-order resulting harmonic, then using that factor to calculate the next highest-order resulting harmonic, and so on. This procedure can be used to calculate Γ, which are essentially renormalized Clebsch–Gordan coefficients, for any combination. Our choice of normalization in the shape functions S, as well as the need for combinations of large angular momenta indices, obfuscates application of standard Clebsch–Gordan tables.

We limit our table of Γ values to p = 2 because the intended application does not require higher-order terms—in our companion article (Styczinski and Harnett, 2020, in prep.), we study the induced magnetic moments for Europa, examining a variety of shapes for the ice–ocean boundary. For planetary bodies, global deviations from spherical symmetry in terms of shape harmonics are expected to be significant only for low-order harmonics. Although more local-scale variation is expected as well, the relative amount χ of each shape harmonic represented in the boundary surface r(θ, ϕ) drops quickly with increasing degree p, and local variation is necessarily represented by large values of p. This is essentially a statement that the spectrum of shape harmonics is very red for planetary bodies—weighted toward low frequencies. That is, the power represented by high-order harmonics is negligible, provided small-scale variations are not periodic on a global scale. Induced magnetic moments resulting from each shape harmonic are proportional to χ, so the low-order shape harmonics will dominate the contributions to the induced moments, too.

Recent work by Vance et al. (2020) explores the magnetic effects of oceanic flows in the Galilean moons, including Europa. These authors find that, under the approximation of a steady field applied by Jupiter, expected oceanic flows induce measurable signals at orbital distances; they do not attempt to solve the more complicated problem of magnetic induction from a varying Jovian field applied to a circulating ocean. The application of a steady-field approximation pairs the results of Vance et al. (2020) well with our results—under the high-conductivity approximation we apply, the time-varying field is entirely excluded from the interior of the ocean, which then experiences a steady field from Jupiter. Combining these results with our model permits calculation of a global induced magnetic field for Europa that is not degenerate in conductivity and boundary shape, even for a single inducing frequency. For highly conducting oceans, no spherically symmetric induction model is capable of breaking this degeneracy without either signals measured at multiple frequencies (Saur et al., 2010; Seufert et al., 2011) or complicated modeling of pressure- and temperature-dependence of subsurface conductivity (Vance et al., 2020). Evaluating a global induction model that accounts for induction by oceanic flows may be important to precisely determine signal provenance for future missions, and failure to include these effects increases the uncertainty of measured ocean properties. However, such a synthesis is beyond the scope of this work.

Application to a specific example

Here we consider a specific example and apply our result to planetary magnetic sounding. Computer code written in Python for the calculations in this section is provided in the Supplementary Material. For calculations relating to the Galileo E26 flyby, the data set GO-J-MAG-3-RDR-HIGHRES-V1.0 was obtained from the NASA Planetary Data System (PDS).

Tobie et al. (2003) considered tidal heating and heat transport through Europa’s ice shell. From numerical models of thermodynamics, these authors concluded that the shell is roughly 22.5 km thick and may be as much as 5 km thicker at the sub- and anti-Jovian points than at midlatitudes. The maxima in thickness can be approximated by choosing in Eq. (46). The minima in thickness at midlatitudes require other shape harmonics, but for simplicity we consider only one real harmonic—this linear combination of shape harmonics is equivalent to choosing r = a + ε sin2 θ cos 2ϕ.

The uniform excitation field Jupiter applies to Europa can be approximated as where the directions are in Europacentric EϕΩ coordinates— is along Europa’s orbital velocity vector and is approximately toward Jupiter’s center. As in the work of Schilling et al. (2007), we take the amplitudes of oscillation to be B = 84 nT, B = −210 nT. From these we obtain B ≈ −425 nT, B = 0, B ≈ −182 nT. We will obtain induced magnetic moments from four cells of Table 1, where the and columns of the excitation field meet the S2–2 and S22 rows of the boundary shape. From Eq. (61), the non-zero induced magnetic moments for this configuration will be

Of the flybys the Galileo spacecraft made past Europa, the E26 flyby had perhaps the most significant magnetic measurements confirming the subsurface ocean there. Data from the E26 flyby ruled out an intrinsic magnetic moment (Schilling et al., 2004; Kivelson et al., 2000), and the strength of the signal encountered has encouraged the interpretation that the subsurface ocean is quite saline (Hand and Chyba, 2007). The E26 flyby had a closest approach altitude of 373 km over Europa (Kivelson et al., 2000). With the induced moments we have calculated for the shape just described, at this distance the induced magnetic field accounts for about 54.25 nT of the observed signal. If we instead assume spherical symmetry in the ice–ocean boundary, at this distance the induced magnetic field accounts for a negligible difference to the expected signal, approximately 54.27 nT. The observed signal therefore differs by only about 0.02 nT in the interpretation of Galileo data based on this boundary shape, and Europa resides in a region where plasma fluctuations routinely force uncertainties on the order of 10 nT (Schilling et al., 2004). It is therefore apparent that the Galileo data are insufficient to resolve the shape suggested by Tobie et al. (2003). However, this does not imply that the Galileo data are incapable of resolving a non-spherical ice–ocean boundary; our companion work (Styczinski and Harnett, 2020, in prep.) explores what limits may be placed from the existing data. The degree of departure from spherical symmetry in the example shape we applied may be quantified by the ratio ε/a, which for this shape was 2.5 km/1537.5 km ≈ 0.0016. A boundary shape deviating more from spherical symmetry will have a proportionally greater effect on the induced magnetic field.

The upcoming Europa Clipper mission will have several flybys within 25 km of Europa’s surface. At this low altitude, but with the same flyby orientation as Galileo E26, the difference in field magnitude from assuming the ice–ocean boundary shape described by Tobie et al. (2003) is instead about 0.06 nT. In order to produce a difference in the magnitude of the induced field of greater than 1.0 nT, assuming the same flyby geometry and shape harmonics in the boundary, E must be at least 24.5 km, with a correspondingly thicker average ice shell. This unlikely scenario would result in an ice shell of vanishing thickness at the leading and trailing points that is 49 km thick at the sub- and anti-Jovian points. Examining the effect on individual components of the magnetic field is more likely to produce a measurable signal—we take a naive approach here only to demonstrate how our model may be applied.

Even for the shape described in Eq. (63), a lander at the sub-Jovian point would likely experience such a measurable signal. Our model predicts a difference of about 1/2 nT in the vertical component () of the magnetic field compared to the spherically symmetric case at peak times during Europa’s synodic period with Jupiter, despite a negligible change in field magnitude. Fig. 2 shows the approximate vertical component of the net magnetic field (excitation field + induced field) experienced by a lander at the sub-Jovian point on Europa’s surface through a synodic period with Jupiter. The vertical component experiences a larger signal because under the high-conductivity approximation, the radial component of the time-varying field is entirely canceled at the outer boundary of the conducting surface (Parkinson, 1983). At Europa’s surface, the conducting boundary is 22.5 km below the lander for the spherically symmetric case, and 25 km below the lander for the non-spherical example case. The additional difference results in induced currents generating the measured magnetic fields that are farther from the measurement point, so the time-varying radial field is not as effectively canceled. This highlights the importance of considering the overall structure of the magnetic field measured by visiting spacecraft, and other harmonics in the boundary shape may have a more pronounced response. In addition, future magnetic measurements from Europa’s surface must account for possible asymmetry in the ice–ocean boundary, lest these signals be incorrectly attributed to another source.

Fig. 2.

The vertical component of the net magnetic field at the sub-Jovian point on Europa’s surface for two cases: an asymmetric ice–ocean boundary with a shape described by Eq. (63), based on Tobie et al. (2003), and a spherically symmetric ice–ocean boundary with the same average ice shell thickness, 22.5 km. The maximum difference, which occurs when Jupiter’s dipole moment nods directly toward or away from Europa, is approximately 1/2 nT for this component of the magnetic field. Our results demonstrate that even modest variation in ice shell thickness must be accounted for in understanding magnetic measurements by a lander on Europa’s surface.

Although our calculated values are small compared to systematic uncertainties for Galileo MAG due to plasma effects, we do not mean to suggest that the magnetic effects of an asymmetric ice–ocean boundary are negligible for orbiting spacecraft. Europa Clipper will have numerous close flybys at varying combinations of latitude, longitude, and System III (Jovian) longitude that will present wider variation in magnetic conditions. Plasma measurements by the PIMS instrument (Grey et al., 2018) will also aid in far better characterization of Europa’s plasma environment than was possible with the compromised Galileo data return, which will support better removal of plasma signals from magnetic data. Furthermore, the overall shape of the induced magnetic field varies with differing induced moments, and changes in global field will be better captured by least-squares comparison to flyby measurements than by individual locations. In this work, we endeavor only to demonstrate how our results may be applied in the study of magnetic sounding. For a more thorough application to Europa in constraining its interior structure, see our companion article (Styczinski and Harnett, 2020, in prep.).

Limits of applicability

The approximations we have made in order to obtain solutions for B limit the applicability of our results. We now discuss the implications for each approximation we make and necessary considerations for the application of our results to physical situations.

Chiefly, we have assumed that kr is effectively infinite; in real applications, it will be large for planetary bodies and good conductors, but not infinite. We are therefore assuming that the induced moments are sensitive only to the shape of the outermost boundary surface. For intended future applications of these results, examining the effects of non-uniform conductivity structures σ(r) expected to be present in ocean worlds, this approximation cannot be applied. Instead, numerical values for kr, j(kr), and y(kr) for layers of approximately uniform conductivity must be inserted into the boundary conditions between each layer. An approach to do this for the spherically symmetric case was first published by Srivastava (1966), and summarized by Parkinson (1983). Concentric near-spherical shapes may be possible to combine with the technique pioneered by Srivastava, further extending the utility of our results. Whether simple superposition of boundary shapes would continue to hold for finite kr is a topic for future study.

The assumption that μ = μ for the conducting body is completely overshadowed by the assumption that kr → ∞, because k is proportional to μ, which can only be μ or greater. In any future applications where kr is finite, non-unity μ/μ must be considered and incorporated, but for planetary applications the difference will still be insignificant.

Also overshadowed by assuming infinite kr is any contribution to the induced magnetic moments by buried conducting layers of much higher or lower conductivity. For example, all known ocean worlds contain rocky mantles beneath their ice and ocean layers. Except when very well hydrated, rocks are typically very poor conductors (Khurana et al., 2002). The finite depth of conducting oceans on these bodies will limit the induction response, especially for very low frequencies, where ∣kr∣ takes smaller values. In addition, deep layers of high conductivity, as in the case of an iron core, which may be present within Europa, could contribute to the induced magnetic field. As described above, numerical calculations would be required to accurately represent finite-depth layers of differing conductivity, whether it grows or shrinks with depth.

We also assumed that the excitation field applied to the conducting body may be treated as uniform. This is a standard assumption, as the complexity introduced by non-uniform excitation fields is great: Magnetic moments matching the order of the spherical harmonics in the non-uniform field are induced, even for the case of a spherically symmetric boundary. For a near-spherical boundary, powers of r in the new terms for the excitation field introduce more mixing between harmonics, because they must be replaced by r(θ, ϕ) in the boundary condition equations. However, deviations from uniformity encountered in physical applications are typically small, or other approximations, like considering only far-field solutions, render the higher-order terms insignificant. The boundary condition equations could be updated to reflect a non-uniform field without difficulty, but that task is left to the applications that may require it.

The assumption that we may safely neglect terms of second order or higher in ε/a is dictated by the sensitivity of measurement equipment for the induced magnetic field. For Europa’s relatively thin ice shell, which could average 30 km thick (Hand and Chyba, 2007), the maximum possible value for ε/a = 30/(1560 – 30) ≈ 0.02 for a nominal 1560 km planetary radius. Differences in the induced magnetic moments compared to the spherically symmetric case will be of this order of magnitude. In practice, the space around planetary bodies is far from a vacuum, and this reduces the need to consider higher-order terms in ε/a. The planets’ extensive magnetospheres contain plasma flows that contribute short-period fluctuations to the measured magnetic fields. These extraneous signals are difficult to model, and even harder to extract from measurements, diminishing the precision that may be afforded in measuring the induced magnetic fields of target bodies. Incorporating higher-order terms in ε/a would be best accomplished by numerically solving the boundary conditions for B, and would likely result in a violation of the superposition of boundary shape harmonics.

We chose to tabulate values only for shapes of harmonic order p = 1, 2 due to rapid decline of the magnitude of χ with increasing p. As explained in Section 4.2, our methods can readily be applied to shape harmonics of higher order, but high-order shape harmonics are not expected for planetary bodies. For applications where higher-order shape harmonics are needed, the necessary Γ values can easily be calculated.

Conclusions

This work represents the first appearance in the literature of a solution to the boundary value problem of the induced magnetic moments for a near-spherical conducting body. Although the approximations we make limit the applicability of our results, the conditions under which these approximations are valid are critical areas of study. Our approximations readily apply to the many ocean worlds scattered throughout the Solar System. Any planetary body with globally distributed conducting material that is subjected to a uniform oscillatory magnetic field may be described by our results. Synergy with other recent work (e.g., Vance et al. (2020)) may even permit independent determination of conductivity and ocean shape, which are degenerate under our approximations and the standard spherically symmetric models of Zimmer et al. (2000) and Schilling et al. (2007). In addition, near-spherical conducting bodies on a laboratory scale may be described by our results, provided they too are subjected to a uniform excitation field and are highly conducting.

Calculations using a simplified model with the results we present here demonstrate that future missions, especially those that land on Europa’s surface, must account for asymmetries in the ice–ocean boundary. Our results suggest that at the surface, signals on the order of nT may be observed for even modest variations in ice shell thickness. Furthermore, these signals are likely to be resolvable at orbital distances, including for Europa Clipper, especially if measurement precision reaches 0.1 nT. A failure to include subsurface asymmetry will result in incorrect attribution of these signals to other processes, and increase uncertainties in valuation of subsurface properties.

This work was motivated by the lack of resources in the literature for this very problem; we apply our results to magnetic sounding of Europa in a companion article (Styczinski and Harnett, 2020, in prep.). Therein, we constrain the departure from spherical symmetry that may be present in the shape of the ice–ocean boundary beneath Europa’s crust, to the extent possible with Galileo data. We accomplish this by comparing the induced magnetic moments for boundary shapes representing characteristic extremes that may be expected from geological, orbital, and thermodynamic constraints to magnetic measurements from the Galileo mission. Calculating the induced magnetic moments using the results of this work made that study possible.