Local spherical harmonic power spectra from local magnetic or gravity data

SUMMARY

We present a method to calculate local spherical harmonic power spectra directly from spacecraft magnetic and gravity data with varying spacecraft altitude. Previously published applications of local spherical harmonic power spectra have been formulated for data that are available at a single collection altitude, such as data evaluated from a global spherical harmonic model. Our approach consists of first solving for local models from local data and then obtaining local multitaper spectra from the local models. We demonstrate with numerical tests that this approach can produce reliable results. Our method is particularly useful in situations where data coverage does not allow for calculating global magnetic or gravity field models, or where data quality or quantity varies regionally and where local models could yield superior resolution or quality over global models.

1 INTRODUCTION

Studies of planetary topography, gravity fields and magnetic fields often make use of spherical harmonic power spectra, or ratios of power spectra to estimate parameters such as elastic thickness or magnetic source depth. This estimation is carried out using analytical expressions for the power spectra that depend on the desired parameters and comparing these analytical spectra to the power spectra (or power spectral ratios) obtained from global topographic, gravity or magnetic measurements.

Planetary properties such as elastic lithospheric thickness or crustal magnetic source depths are expected to vary geographically, warranting localized studies of physical properties. To calculate local spherical harmonic power spectra, the method of Wieczorek & Simons (2005, 2007) and Dahlen & Simons (2008) multiplies gravity, topography or magnetic field data available at a single altitude with tapering functions of values close to zero outside of the region of interest. This approach is a spherical incarnation of the one-dimensional local multitaper spectrum of Thomson (1982). The local multitaper spherical harmonic power spectrum of Wieczorek & Simons (2005, 2007), and Dahlen & Simons (2008) is widely used to calculate regional variations of lithospheric elastic thickness from topography and gravity fields (e.g., Belleguic & Lognonné 2005; Grott & Wieczorek 2012; James et al. 2015; Genova et al. 2016; Goossens et al. 2022), see the overview by Wieczorek (2015) and to calculate magnetic source depths from magnetic fields (e.g., Lewis & Simons 2012; Wieczorek 2018; Gong & Wieczorek 2021). These applications used global spherical harmonic expressions (we call such expressions ‘models’ for the remainder of this paper) for the gravity field, the topography and the magnetic field. In principle, the method of Wieczorek & Simons (2005, 2007) and Dahlen & Simons (2008) could be applied directly to measurements of gravity, topography and magnetic fields, but the data would need to be at constant altitude and free from undesired signals such as external magnetic fields. Furthermore, the construction of high-fidelity global models depends on global data availability and ideally homogeneous quality.

Multiple recent and current spacecrafts have orbits that do not cover the surface of a planetary body evenly, for example, MESSENGER (Solomon et al. 2007), MAVEN (Jakosky et al. 2015) and Europa Clipper (Campagnola et al. 2019). The strong altitude variations and gaps in coverage of such missions require calculation of local power spectra from regional data with varying spacecraft altitude to extract spectral information. Here, we present a method for this purpose (Section 3). Our approach uses a combination of two previously published methods. In a first step, we calculate local models from local data using the method of Plattner & Simons (2017). This approach yields spherical harmonic coefficients that represent the data regionally only. In a second step, we calculate local power spectra from the local models using the approach of Wieczorek & Simons (2005, 2007) and Dahlen & Simons (2008). Our numerical demonstrations (Section 4) show that correct local multitaper spectra can be obtained from regional data with varying spacecraft altitude. A modified version of our method was used by Plattner & Johnson (2021) to extract source depth information of a core magnetic field feature on Mercury from regional data with varying spacecraft altitude. We discuss this application in Section 4.2.

Alternative methods for obtaining local spectra exist, such as projecting spherical models onto the plane and applying Cartesian tapering (e.g., McKenzie 1994), the spherical localization approach of Simons et al. (1997), spherical wavelet power spectra (e.g., Gerhards 2014; Audet 2014) and spectra obtained from Revised Spherical Cap Harmonic Analysis ‘R-SCHA’ (e.g., Thébault & Vervelidou 2015; Vervelidou & Thébault 2015; Thébault et al. 2018). The advantage of the multitaper method of Wieczorek & Simons (2005, 2007) is that it permits the direct comparison of data-derived local models with the myriad of analytical spherical harmonic power spectra, or ratios of spectra, that depend on parameters such as magnetic source depth (e.g., Jackson 1990; Voorhies et al. 2002; Langlais et al. 2014; Wieczorek 2018) or lithospheric parameters (e.g., Wieczorek & Phillips 1997; Wieczorek 2007), because Wieczorek & Simons (2005, 2007) provided a method to calculate expected multitaper local spectra from analytical global spectra. Here, we build on the success of the method of Wieczorek & Simons (2005, 2007) by extending it to local spacecraft data with varying altitude.

2 PREREQUISITES

We first provide a summary of the methods upon which our approach builds. Specifically, we point out the difference between classical Slepian functions (Simons et al. 2006) and altitude-cognizant Slepian functions (Plattner & Simons 2017) in Section 2.3. We use the classical Slepian functions as tapers to calculate local power spectra from spherical harmonic models (Section 2.2) and we use the altitude-cognizant Slepian functions to calculate local models from regionally available data at varying spacecraft altitude (Section 2.3). Both the classical and the altitude-cognizant Slepian functions are constructed from spherical harmonics (Section 2.1).

2.1 Spherical harmonics and global power spectrum

Gravitational acceleration as well as magnetic fields in regions without electrical currents can be described as the spatial gradient of a scalar potential (e.g., Blakely 1995). If additionally, all mass or magnetic sources are within a planetary body of radius r_p, then outside of the planetary body, the scalar potential Φ can be described as a linear combination of spherical harmonics

where r, θ and φ are the radial position, colatitude and longitude, r_p is the planetary body’s radius, c_lm are the spherical harmonic coefficients describing the potential, L₀ is the lowest spherical harmonic degree (L₀ = 0 for the gravitational potential and L₀ = 1 for the magnetic potential) and Y_lm are the scalar spherical harmonics

where

for spherical harmonic degree l and order m. The spherical harmonic power spectrum for the coefficients c_lm is defined as

Gravity as well as magnetic potentials have a discipline-specific normalization for the spherical harmonics and the coefficients. For the gravitational potential, the convention is to use fully normalized spherical harmonics and normalize the coefficients to have a factor including the planetary mass and Newton’s gravitational constant before the summation in eq. (1), see, for example, Blakely (1995). The study of magnetic potentials typically uses the Schmidt quasi-normalized associated Legendre functions

where δ(m, 0) equals 1 for m = 0 and is 0 otherwise. For sources within the planetary body, a magnetic potential V(r, θ, φ) can be described outside of the planetary body as

where r_p is the planet’s radius and the coefficients g_lm and h_lm are the Gauss coefficients (e.g., Blakely 1995). The relationship between the Gauss coefficients g_lm and h_lm of eq. (6) and the spherical harmonic coefficients c_lm from eq. (1) is

Studies of the magnetic potential typically use the Mauersberger–Lowes power spectrum

which is related to the power spectrum S_l of eq. (4) through the renormalization

2.2 Local spectrum from global model

Spherical harmonic power spectra (eqs 4 and 8) express the global distribution of power per spatial wavelength. A local power spectrum for a chosen region can be calculated using the approach of Wieczorek & Simons (2005, 2007) and Dahlen & Simons (2008), which we summarize below.

An intuitive but inadequate approach consists of calculating a local power spectrum from the spherical harmonic coefficients of a potential field with values outside of the chosen region set to zero (‘boxcar tapering’). Multiplication in the spatial domain corresponds in the spectral domain to an operation reminiscent of convolution (Wieczorek & Simons 2007, their eq. 2.11). Boxcar functions for spherical caps (value 1 inside the cap and value 0 outside) have non-zero powers over a wide range of spherical harmonic degrees (Dahlen & Simons 2008) and this observation extends to boxcar functions with shapes other than a spherical cap. Thus, for local power spectra obtained from multiplication of potential fields with boxcar functions, each degree power will contain the spherical harmonic coefficients of all degrees of the global potential field.

Ideal tapering functions should have spectra that span a degree range as narrow as possible while spatially focusing on the region of interest (sometimes referred to as spatio-spectral localization). Simons et al. (2006) and Simons & Dahlen (2006) describe a construction of spatially concentrated yet spectrally limited functions for generic regions (‘classical Slepian functions’). For a chosen region R on the sphere Ω and a maximum spherical harmonic degree L_tap ('tapering bandwidth'), classical Slepian functions are constructed from linear combinations of spherical harmonics

by solving the optimization problem that concentrates the classical Slepian functions in the spatial region of interest

Solving this optimization problem is equivalent to obtaining the eigenvector with the largest eigenvalue for the eigenvalue problem

where the (L_tap + 1)² × (L_tap + 1)² matrix |$\bf K$| contains the entries

and the vector |$\bf q$| contains the coefficients q_lm. The eigenvalue λ indicates the relative fraction of energy of the function q(θ, φ) inside the region R. The construction from eq. (12) yields not one optimal tapering function, but (L_tap + 1)² functions q_j(θ, φ) with j = 1, 2, …, (L_tap + 1)², ordered by their eigenvalues |$\lambda _1 \ge \lambda _2 \ge \ldots \ge \lambda _{({L_\text{tap}}+1)^2}$| (Fig 1a–d).

A local Slepian-tapered spectrum |$S_l^{(1)}$| of a function f(θ, φ) can thus be calculated from the spherical harmonic coefficients of the product f(θ, φ)q₁(θ, φ), which spatially tapers the function f(θ, φ) using the best-concentrated classical Slepian function for the chosen region R and taper spectral bandwidth L_tap. The tapering bandwidth L_tap determines the bandwidth of spectral convolution, with the trade-off that reducing the width of the spectral convolution will increase the spatial footprint of the tapering function and thus of the local spectrum.

Because of the spatial pattern of the classical Slepian function q₁(θ, φ) (see e.g., Figs 1a and c), the ‘single-taper’ local spectrum |$S_l^{(1)}$| samples the data in the region R unevenly (Wieczorek & Simons 2005). This can be remedied by calculating single-taper local spectra |$S_l^{(j)}$| from the product of a function f(θ, ϕ) with each taper q_j(θ, ϕ) and then averaging these to create a multitaper local spectrum

The weights w_j for the single-taper local spectra S^(j) need to sum to 1 (Wieczorek & Simons 2005). We use

A Mauersberger–Lowes multitaper local power spectrum |$R^\text{mt}_l$| can be calculated from the multitaper power spectrum |$S^\text{mt}_l$| using eq. (9).

To compare the performance of multitaper and boxcar local spectra, we create an example function with a spectrum containing most of its energy at four spectral peaks at degrees l = 10, 20, 30 and 40 (Fig. 2). This function is created as a sum of four functions F₁, F₂ and F₃. Function F₁ is constructed from spherical harmonic coefficients that are zero for all spherical harmonic degrees except degree l = 10 and has normally distributed random coefficients for all orders m for degree l = 10. To localize this function to the Southern Hemisphere, we multiply the spatial pattern of F₁ by the spatial pattern of the single best concentrated classical Slepian function for R = the Southern Hemisphere and L_tap = 1. The function F₂ is constructed using the same approach, but for the Northern Hemisphere, where the spatial expansion of spherical harmonic coefficients that are zero except for all orders m for l = 40 is multiplied by the single best classical Slepian function for R = the Northern Hemisphere and L_tap = 1. The nonzero coefficients are random normally distributed. The function F₃ is created from a spherical harmonic expansion with all coefficients set to zero except those for all orders of spherical harmonic degree l = 20 as well as degree l = 30. The nonzero coefficients are random normally distributed. The spatial expression of these coefficients is localized to R = North America by multiplication with the single best concentrated classical Slepian function for L_tap = 2. The spatial expression of F₁ + F₂ + F₃ is shown in Fig. 2a and its spectrum in Fig. 2d. The localization of each function F₁, F₂ and F₃ leads to a spectral spreading of each single-degree peak over the L_tap + 1 degrees greater and less than each peak. For functions F₁ and F₂, the spectral peak is spread over degrees l = 8 to 12, and degrees l = 38 to 42, respectively. For the function F₃, the spectral peaks are spread over degrees l = 17 to 23 and l = 27 to 33 (Fig. 2d).

For this function F₁ + F₂ + F₃, we calculate local spectra for North America using the boxcar approach (Fig. 2b) as well as the multitaper local spectrum approach (Fig. 2c) for tapering bandwidth L_tap = 4. We expect that a spectrum of this function focusing on North America should show peaks for l = 20, 30 and 40, but not the peak for l = 10, because that spectral peak is focused on the Southern Hemisphere. We also expect that the peaks of spherical harmonic degrees l = 20 and 30 have correct relative strengths, because both fields are equally concentrated over North America. Finally, we expect that the local spectrum contains a peak for l = 40 that is weaker than the corresponding peak in the global spectrum, because the signal with peak at l = 40 is distributed over the entire Northern Hemisphere. The boxcar local spectrum (Fig. 2d) shows a small but non-zero peak for l = 10. This is an artifact caused by the spectral aliasing of the broadband boxcar tapering. The boxcar local spectral peak at degree l = 20 is greater than its peak at l = 30 even though the original spectral peak at l = 20 is smaller than its peak at l = 30. For the boxcar local spectrum, parameter estimations that rely on relative spectral power per degree, such as magnetic field local source depths (e.g., Jackson 1990; Voorhies et al. 2002; Langlais et al. 2014; Wieczorek 2018) may yield spurious results because of wrong relative powers. On the other hand, the multitaper local spectrum (Fig. 2d) does not show substantial power for l = 10 and correctly represents the relative powers of the peaks at l = 20 and 30.

2.3 Local models from local data

To calculate single-taper and ultimately multitaper local spectra of a potential field Φ in eq. (1) or V in eq. (6), we need to obtain spherical harmonic coefficients of the products of the potential field with the single taper function. If measurements of the potential field or its radial derivative are available at a constant altitude, then obtaining the spherical harmonic coefficients of this product may be possible directly from the data. In many situations, however, the data are distributed over a range of altitudes and may contain contributions other than the targeted internal field. It may thus be desirable to first obtain a model of the (internal) potential field and calculate the multitaper local power spectrum from this model. If global data coverage is available and data quality is globally (sufficiently) homogeneous, then global spherical-harmonic models can be calculated from the data and used to evaluate the potential field on a grid for tapering. If, however global data coverage is not available or data quality is spatially varying, then we need to resort to a local approach.

One could consider using the classical Slepian functions for high maximum spherical harmonic degrees as alternatives to the spherical harmonics to solve for local potential fields. Indeed, if only the few classical Slepian functions that are well concentrated within the data region are used, then the resulting potential field will be local. However, to optimally concentrate the classical Slepian functions within their region R, substantial energy is allocated to their high spherical harmonic degrees. For data at spacecraft altitude, such focus on high spherical harmonic degrees leads to noise amplification under downward continuation because of the factor (r_p/r)^{l + 1} in eqs (1) and (6).

To avoid this issue, Plattner & Simons (2017) developed spatially concentrated functions that also minimize noise amplification caused by the factor (r_p/r_s)^{l + 1} for a mean data radial position r_s. For a maximum spherical harmonic degree L_max, a target region R, and a mean data radial position r_s, the ‘altitude-cognizant Slepian functions’ of Plattner & Simons (2017) can be constructed by solving the eigenvalue problem

where the (L_max + 1)² × (L_max + 1)² matrix |${\bf \tilde{K}}$| contains the inner products

of the functions

Eigenvectors t_j of |${\bf \tilde{K}}$| for 1 ≤ j ≤ L_max contain the coefficients t_{j, lm} which can be expanded to the potential field expression of the altitude-cognizant Slepian functions

and to their corresponding spatial-derivative expression

As for the classical Slepian functions in Section 2.2, the construction in eqs (16)–(20) yields (L_max + 1)² altitude-cognizant Slepian functions, ordered in decreasing order by their eigenvalue μ_j with j = 1, …, (L_max + 1)². The altitude-cognizant Slepian functions with the largest eigenvalues are best suited to solve for a potential field on the planet’s surface from spacecraft data within the region R while minimizing noise amplification under downward continuation. Subsequent (those with smaller eigenvalues) altitude-cognizant Slepian functions incorporate an increasing amount of energy in the higher spherical harmonic degrees. This allows for resolving finer structures but is also more prone to noise amplification. Altitude-cognizant Slepian functions for the lowest eigenvalues resolve the region complementing R and focus on the highest spherical harmonic degrees. Examples of altitude-cognizant Slepian functions for R = Africa, L_max = 40 and r_s = 6671 km with r_p = 6371 km are visibly different from the classical Slepian functions with the same spherical harmonic bandwidth L_tap = 40 constructed for the same region (Figs 1c–f).

To calculate a potential field Φ on a planet’s surface within a region R, given spacecraft data of the potential’s spatial derivative |$-\boldsymbol{\nabla } \Phi$| within R only, we can use the first J altitude-cognizant Slepian functions for eigenvalues μ₁ ≥ … ≥ μ_J to solve the linear system of equations

for the coefficients b_j for 1 ≤ j ≤ J, where (r_i, θ_i, ϕ_i) are the data locations, with 1 ≤ i ≤ n_data. The resulting potential field Φ for the coefficients b_j for 1 ≤ j ≤ J,

is representative of the true potential field within the region R only and tends towards zero outside of R.

The choice for the number J of altitude-cognizant Slepian functions to use in eq. (22) is data specific and problem dependent and can be considered a regularization parameter of the inverse problem. Choosing J too small will produce a potential field that is not sensitive to short-scale noise, but which may not reproduce short-scale features of the potential field. On the other hand, choosing J too large may create artifacts in the solution. Section 4 presents a strategy to choose a suitable value for J.

The functions of eqs (16)–(20) originally presented by Plattner & Simons (2017) were constructed for the problem of solving for a potential field on the planet’s surface from spatial derivative data (magnetic |$\boldsymbol{B}$|-field) at spacecraft altitude. For alternative problems, such as calculating a local potential field model on the planet’s surface from potential field data at spacecraft altitude, the derivation can be adapted by omitting the spatial derivative in eq. (18).

3 IMPLEMENTATION: LOCAL SPECTRUM FROM LOCAL DATA

Given data only within a region R, we can calculate a local multitaper power spectrum for the region R following the two steps

• calculate a local potential field Φ using the approach of Section 2.3,

• calculate a multitaper local power spectrum |$S^\text{mt}_l$| from Φ using the approach of Section 2.2.

The resulting local multitaper power spectrum depends on several choices. The region R is typically determined by data availability or quality, or geographic choice. Only data within the region R should be used when solving for Φ in step (i) to avoid fitting altitude-cognizant Slepian functions that have small values outside of the region R to data outside of R. The maximum spherical harmonic degree L_max in step (i) should be chosen to be sufficiently large to resolve the desired spherical harmonic degrees but should not exceed limitations imposed by the spatial density of the data. An optimal value for the number J of altitude-cognizant Slepian functions for step (i) can be found using the misfit minimization approach shown in Section 4. The value for r_s in step (i) is typically chosen as the mean data radial position. Step (ii) additionally requires choosing the tapering bandwidth L_tap. Each degree power |$S^\text{mt}_l$| in the multitaper spectrum is correlated with the powers for degrees l − L_tap to l + L_tap. The value for L_tap thus needs to be chosen to be sufficiently small to yield a sufficiently large number of uncorrelated degree powers in the multitaper spectrum. On the other hand, L_tap needs to be sufficiently large to provide at least a few well-concentrated classical Slepian functions for the region. The region R used for step (ii) does not need to be the same as the region used for step (i), but can be a subregion, as shown in Section 4. Appendix  A details how analytical spectra such as those presented by Jackson (1990), Voorhies et al. (2002), Langlais et al. (2014) and Wieczorek (2018) can be used to extract spectral parameters (e.g., source depth) from multitaper local spectra.

4 PROOF OF CONCEPT AND APPLICATION EXAMPLE

The examples in this section demonstrate that our approach can successfully calculate a locally representative potential field as well as a local power spectrum from spatially limited data with varying radial position (Section 4.1). We also summarize an application of a modified version of this approach (Section 4.2).

4.1 Proof of concept

To be able to show that the resulting fields and spectra are correct, we need to know the ‘true’ potential field that gives rise to the data. To this end, we use the magnetic field model of Langlais et al. (2019) as the ‘true’ potential field and create ‘synthetic’ data by evaluating the spatial derivative of the potential field at data locations that simulate spacecraft data available within a region only.

The radial position of our data has a uniform random variation between 300 and 350 km above Mars’ mean radius of r_p = 3393.5 km. We limit the data locations to a spherical cap of semi-opening angle 30° centered at longitude 160°W, latitude 69°S. We set the data locations to be along tracks with one degree longitudinal spacing and sample each track every 0.2 degrees. To simulate noise, we added zero mean Gaussian random values of standard deviation 10% of the standard deviation of the corresponding magnetic field component to the data. The added random noise standard deviations are 4.3 nT for the radial component, 3.7 nT for the colatitudinal component and 2.0 nT for the longitudinal component (Fig. 3a).

From these data, we aim to solve for a regional representation of the magnetic field within the region of data availability, as well as for a regional spherical harmonic power spectrum for a smaller cap of semi-opening angle 20° centered on the same location as the spherical cap of data availability. We select the maximum spherical harmonic degree L_max = 134 to solve for the regional magnetic field from the synthetic data, because this is the maximum spherical harmonic degree of the model that gave rise to the synthetic data.

Similar to typical spacecraft data, our synthetic data have a higher spatial resolution along the tracks than between the tracks and the spatial data density is higher closer to the poles than at lower latitudes. To ensure that the uneven spatial sampling does not affect our results, and to create a way of assessing our results using a bootstrapping approach, we subsample the data using an equal area random selection approach (Figs 3b and c). The choice for the subsampling retention percentage should counteract the spatially uneven sampling caused by the spacecraft tracks. Here, we present two choices for retention percentages, 30% and 10%.

To find an optimal number J of altitude-cognizant Slepian functions, we calculate models for a range of values for J using the retained data from the subsampling and obtain the model root mean squared error (rmse) using the data points excluded from the model calculation. We repeat this 100 times and use the minimum of the average rmse curve to find an optimal value for J (Figs 3e and f). We observe an increase in rmse with high values for J because structure that is compatible with individual subsampled data sets and which is retained in the models with large numbers of Slepian functions is not required by the excluded data. Hence, selecting J corresponding to the minimal rmse yields a model with structure that can be resolved by the retained data and which is required by all data. This effect is more pronounced when a lower percentage of data is retained, leading to less structure that can be resolved by the retained data. For our example, for equal area random subsampling retaining 30% of the data we obtain J = 549 (Fig. 3d). For 10% random subsampling, J = 458 (Fig. 3e). Higher data altitudes also lead to lower values for J because the signal is upward continued and thus less small-scale structure is recoverable from the data.

We demonstrate our method for an equal area random data subselection with 30% data retention (Fig. 4). For J = 549, we calculate models for 100 equal area randomly subselected data (each 30%) to obtain a mean model (Fig. 4a) and standard deviation (Fig. 4b). The spatial structure of the mean regionally representative spherical harmonic model (Fig. 4a) resembles the global model that was used to create the synthetic data (Fig. 4c) within the region R, but it is missing shorter spatial scales as well as the signal outside of the region, as is visible in the difference between the mean model and the original global model (Fig. 4d).

For each of the regionally representative spherical harmonic models, we calculate a multitaper local power spectrum for a smaller region R′, which is concentric with the region R, but has a semi-opening angle of 20°. We use a tapering bandwidth L_tap = 8, which allows for two classical Slepian functions (out of 81) with more than 50% of their power within the spherical cap of semi-opening angle 20°. This choice for the spherical cap and tapering bandwidth is identical to the choice made by Lewis & Simons (2012). The local power spectra (Fig. 4e) obtained from the data (Fig. 3d) match the local power spectrum obtained from the true model up to spherical harmonic degree 70 (Figs 4e and f). Between spherical harmonic degree 70 and 85, the local spectra obtained from the synthetic data show an increasing divergence among individual solutions (Fig. 4f). For degrees greater than 85, the local spectra obtained from the synthetic data have substantially less power than the local spectrum obtained from the model of Langlais et al. (2019).

We next demonstrate that the discrepancy for high spherical harmonic degrees between the local spectra based on the data and the local spectrum from the model of Langlais et al. (2019) (Figs 4e and f) is caused by the spacecraft altitude, not by the method itself. For this purpose, we create a synthetic data set on the planet’s surface r = r_p = 3393.5 km without altitude variation, but with added Gaussian noise of standard deviation 10% of the standard deviation of the data for each component (115.4 nT for the radial component, 91.4 nT for the colatitudinal component and 69.0 nT for the longitudinal component). We use the same approach as for the demonstration in Fig. 4 for an equal area random data subselection percentage of 30%. This subselection percentage yielded an optimal value of J = 728 for the number of altitude-cognizant Slepian functions to be used. The resulting mean local model matches the model of Langlais et al. (2019) well within the region of data availability and is close to zero outside (Figs 5a and b). The local spectra calculated from the synthetic data fit the local spectrum of the model of Langlais et al. (2019) well (Figs 5c and d).

In Figs 4 and 5, the region of data availability is larger than the region for which we calculate the spectrum. This is chosen to account for the footprint of the local spectral calculation in which the field outside of the region influences the result. This influence diminishes with distance from the region, depending on the tapering bandwidth L_tap. Our third demonstration (Fig. 6) investigates the spectrum obtained when the region of data availability equals the region for which we calculate the spectrum. For that purpose, we use the same data as for our second demonstration (Fig. 5) but remove data outside of the spherical cap of semi-opening angle 20°, centered at longitude 160°W, latitude 69°S. This is the spherical cap we had used to calculate the spectrum in demonstrations 1 and 2 (Figs 4 and 5).

To calculate a regionally representative model and local power spectrum from these synthetic data, we use the same approach as for Figs 4 and 5. We use equal area random subselection retaining 30% of the data within the spherical cap of semi-opening angle 20° to calculate local models using a range of values for J, the number of altitude-cognizant Slepian functions constructed for this region and average spacecraft radial position r = r_p, with L_max = 134. We repeat this process 100 times to find the optimal number J = 824 of altitude-cognizant Slepian functions as described in Section 3 and used in Figs 4d and e. For J = 824, we then calculate 100 local models from 100 equal area randomly subselected data with 30% retention rate to obtain an average regionally representative model (Fig. 6a). This average model fits the model of Langlais et al. (2019) well within the region of data availability and is close to zero outside this region (Fig. 6b). The local power spectra obtained from the local data are different from the local power spectrum obtained from the model of Langlais et al. (2019), particularly at spherical harmonic degrees less than 25, as well as between degrees 40 and 75 and for degrees greater than 110 (Figs 6c and d). These discrepancies are a consequence of the power spectral localization effectively including contributions from outside the region of interest that are not represented in the local model.

Finally, to study the effects of the choice of the retention percentage for the equal area random data subselection on the local models and spectra, we repeat the demonstration of Fig. 4 for 10% subsampling (Fig. 7) instead of 30% subsampling. The resulting local spectrum faithfully follows the power spectrum of Langlais et al. (2019), but deviates at a lower spherical harmonic degree of 62 instead of 70 (cf. Figs 7c and d to Figs 4e and f) as expected by the reduced complexity captured in the local model.

4.2 Application example

Plattner & Johnson (2021) used a modified version of the method of Section 3 to determine the source depth of a magnetic field anomaly on the planet Mercury from MESSENGER spacecraft data. We discuss this implementation of the method of Section 3 further here to demonstrate its utility in a case with a challenging data set. After subtraction of the magnetospheric model of Korth et al. (2017) from the data, as well as subtraction of large-scale magnetic field patterns in the local time frame, Plattner & Johnson (2021) used non-zonal altitude-cognizant Slepian functions to calculate local magnetic field models. These non-zonal altitude-cognizant Slepian functions were constructed from spherical harmonics with order m ≠ 0, to account for having subtracted axisymmetric magnetic signals from the data. In contrast, the altitude-cognizant Slepian functions in Section 2.3 are based on spherical harmonics for all orders m.

The challenging analysis of large-scale magnetic fields in the MESSENGER data required a different approach than shown in Section 4.1 for the determination of J, the optimal number of altitude-cognizant Slepian functions to use in the local model calculations. For each subsampled data set, Plattner & Johnson (2021) selected those values for J which yielded the steepest local spectrum for their range of resolved spherical harmonic degrees (see detailed discussion in Plattner & Johnson 2021). From these local spectra, Plattner & Johnson (2021) calculated source depths by solving for the source radius using a localized non-zonal analytical spectra. The localized non-zonal spectra were obtained from the global analytical spectrum of Langlais et al. (2014) using the localization method of Wieczorek & Simons (2005, 2007), which we summarize in Appendix  A. Synthetic tests comparable to those of Section 4.1 allowed Plattner & Johnson (2021) to demonstrate that their approach yielded correct results.

Plattner & Johnson (2021) found that the source depth of their magnetic anomaly placed it near the core–mantle boundary of Mercury and further estimated an uncertainty in this source depth. This result presented the first study of a non-axisymmetric anomaly in Mercury’s core magnetic field. The magnetic field anomaly spatially overlaps with a crustal topographic bulge, which was linked to a deep seated mantle anomaly (James et al. 2015). This spatial overlap could be indicative of a connection between core, mantle and crustal dynamics made possible by Mercury’s thin mantle.

The same data of Mercury’s magnetic field that was used by Plattner & Johnson (2021) was also used by Thébault et al. (2018), who used a different method to build a local field model. Thébault et al. (2018) used R-SCHA (Thébault et al. 2006; Vervelidou & Thébault 2015), which solves for a regional potential field inside a spherical cone using boundary conditions that satisfy null flux conditions across the region boundaries (Thébault et al. 2006). R-SCHA model coefficients can be related to spherical harmonic coefficients via infinite sums (e.g., Thébault et al. 2018). However, Thébault et al. (2018) did not extract a source depth from their R-SCHA power spectrum and did not present analytical parameterized power spectra that would allow for source depth calculations directly from R-SCHA spectra. Instead Thébault et al. (2018) noted that the resulting R-SCHA spectrum evaluated at a depth of 400 km appeared to decrease slightly as a function of a spatial scale parameter, which was interpreted to be in general agreement with Mercury’s core–mantle boundary as determined from gravity measurements.

5 DISCUSSION AND CONCLUSIONS

We present a method to calculate local power spectra directly from local data arising from potential fields such as magnetic or gravity fields. Our demonstrations (Section 4.1) show that this method can calculate local multitaper power spectra from local data with varying altitude. We note that these local spectra do not represent the degree powers globally, but within the region of data availability only. Because of this, the resulting local spectra need to be compared to localized spectra, such as expected multitaper spectra, not to global power spectra. A modified approach of our method allowed for the study of a core magnetic field feature on Mercury (Section 4.2, Plattner & Johnson 2021), by comparing the data-derived multitaper local spectra to expected local spectra calculated from analytical parameterized spectra (Appendix  A). These results elucidate the usefulness of this approach, which is advantageous given the increasing number of spacecraft orbit or flyby designs that collect data with uneven regional and altitude coverage.

Our demonstrations (Section 4.1) show that local spectra from global spherical harmonic models tend to incorporate variations outside the target region. Our local modeling dampens the field outside the region and thus reduces the effect of the field outside the region on the local power spectrum. This is apparent when comparing the spectra obtained from local models for a larger region (Fig. 5c) to the spectra obtained from local models for the same region as the spectral calculation (Fig. 6c). However, the manner in which the locally modeled field transitions to near zero outside of the target region depends on the spatial pattern of the altitude-cognizant Slepian functions. Thus, if the region for the spectral calculation is chosen to be the same as the region for the local model, then the spatial pattern of the altitude-cognizant Slepian functions may affect the local power spectrum. This effect can be reduced if the region for the spectral calculation is chosen to be a subregion of the region of local model calculation, as we did in the first two demonstrations (Figs 4 and 5).

Equal area random data subsampling plays a crucial role in our demonstrations. It is used to transform clustered data (higher sampling along tracks) into an equal area data distribution. It is used also to estimate an optimal regularization parameter (the number J of altitude-cognizant Slepian functions) in the model calculation, to create a range of local models, and ultimately, to allow for an assessment of uncertainty in parameters calculated from local spectra. The choice of subsampling density affects the short spatial scales of the resulting local models and thus the local power spectrum portion for the higher spherical harmonic degrees (cf. Fig. 4 to Fig. 7). The choice for the retention percentage of the random equal area subselection should depend on the spatial distribution of the available data. We recommend repeating the calculation of local spectra and their parameters for different choices of subsampling percentage to assess the effect of this choice.

When applying this method to real data, the local model approach may need to be adapted to suit the special needs of each data set. For magnetic field data, we may need to solve for internal and external magnetic fields simultaneously (Plattner & Simons 2017) and calculate the power spectrum from the internal field. When solving for a core magnetic field model from a sufficiently extensive data coverage or low maximum spherical harmonic degree, the best choice for J may be (L_max + 1)², the full set of altitude-cognizant Slepian functions. This would render step (i) in Section 3 equivalent to solving for spherical harmonics instead of altitude-cognizant Slepian functions.

In principle, the presented method can be modified to use other methods such as R-SCHA (Thébault et al. 2006; Vervelidou & Thébault 2015; Thébault et al. 2018) instead of altitude-cognizant Slepian functions to obtain local models, and then calculate local spherical harmonic power spectra based on those local models. Section 4.1 demonstrates that local spherical harmonic power spectra from models based on altitude-cognizant Slepian functions are equivalent to local spherical harmonic power spectra from global spherical harmonic models. An open question is whether local spherical harmonic power spectra from R-SCHA or other local models show the same equivalency.

For gravity fields from satellite-to-satellite tracking, such as GRACE (Tapley et al. 2004), GRACE-FO (Flechtner et al. 2014) and GRAIL (Zuber et al. 2013), an energy balance approach (Jekeli 1999; Guo et al. 2015; Shang et al. 2015; Zeng et al. 2015) can be used to turn inter-satellite range variations into potential differences. To study such data, the local model approach in Section 2.3 can be adapted as described in that section to solve from potential field measurements instead of spatial derivatives of the potential field.

ACKNOWLEDGEMENTS

The material is based upon work supported by NASA under award no. 80NSSC19K1426. AMP thanks Frederik Simons for discussions that led to this research.

DATA AVAILABILITY

Software developed for this paper is available from the website https://www.doi.org/10.5281/zenodo.7838081 and is built on the software https://www.doi.org/10.5281/zenodo.598177.

References

APPENDIX A: ESTIMATION OF SPECTRAL PARAMETERS

Comparing a power spectrum (eqs 4 and 8) to an analytical spectral expression (e.g., Jackson 1990; Voorhies et al. 2002; Langlais et al. 2014; Wieczorek 2018) allows estimating parameters such as the magnetic source depth. To extend such comparisons to the realm of local power spectra, Wieczorek & Simons (2005, 2007) and Dahlen & Simons (2008) presented an approach to calculate expected local multitaper power spectra |$\tilde{S}_l$| from global power spectra S_l by multiplication with a coupling matrix. The following formulation is valid for our choice of multitaper weights w_j (eq. 15) and for the spectral normalization in eq. (4),

where the coupling matrix

In eq. (A2), the Wigner 3-j symbols |$\binom{{l}\,\,{p}\,\,{l'}}{{0}\,\,{0}\,\,{0}}$| are zero, unless l + p + l′ is even and |l − p| ≤ l′ ≤ l + p (‘triangle condition’). If both of these conditions are met, then

where q = (l + p + l′)/2 and

see Messiah (1962, their eqs C.23a and C.23b). We note that due to the triangle condition, all elements of |$M_{l,l^{\prime }}$| with |l − l′| > L_tap are zero and thus in the expected multitaper spectrum in eq. (A1), spherical harmonic degrees l and l′ with |l − l′| > L_tap are not mixed.

To apply this formulation to the spectrum R_l used by the magnetic community (eq. 8), we first renormalize the analytical expression R_l to obtain the global spectrum S_l using eq. (9). Then, we calculate the expected multitaper spectrum |$\tilde{S}^\text{mt}_l$| from S_l using eq. (A1). Finally, we transform |$\tilde{S}^\text{mt}_l$| to the Mauersberger–Lowes normalization |$\tilde{R}^\text{mt}_l$| using eq. (9). Magnetic fields do not have monopoles, thus for global magnetic fields R₀ = 0. However, when calculating local multitaper spectra R^mt and expected multitaper power spectra |$\tilde{R}^\text{mt}$| for L_tap > 0, the power at degree l depends on degree powers |$R_{l-{L_\text{tap}}}$| to |$R_{l+{L_\text{tap}}}$|⁠. For local spectral analysis, we thus include |$R^\text{mt}_0$|⁠.

Calculating such expected multitaper spectra from global spectra requires two assumptions. First, the (unknown) global function pertaining to the coefficients that yield the power spectrum is generated by a globally uniform and isotropic process. Second, the coefficients that yield the power spectrum are random variables with zero mean and a variance depending on the degree l (Wieczorek & Simons 2005). This may not be the case for fields arising in nature. For Earth’s core magnetic field (Finlay et al. 2016), the expected Mauersberger–Lowes multitaper spectrum |$\tilde{R}^\text{mt}_l$|⁠, calculated from the global spectrum R_l for tapering bandwidth L_tap = 3 closely follows the multitaper local spectrum |$R^\text{mt}_l$| for a spherical cap of semi-opening angle 20° centred on the prime meridian at latitude 60°N (Fig. A1a). On the other hand, for Mars’ crustal field and the same tapering bandwidth L_tap = 3, the expected Mauersberger–Lowes multitaper spectrum |$\tilde{R}^\text{mt}_l$| differs substantially from the multitaper local spectrum |$R^\text{mt}_l$| for a spherical cap of semi-opening angle 20° centred on longitude 160°W and latitude 69°S (Fig. A1b). This is a consequence of the variation of Mars’ crustal field over the planet’s surface. Mars’ Southern Hemisphere has a stronger magnetic field than its Northern Hemisphere. As a result, the multitaper local power spectrum focusing on the spherical cap in the Southern Hemisphere has more power than the expected multitaper power spectrum averaging the global magnetic field (Fig. A1b).