# On the accuracy of palaeopole estimations from magnetic field measurements

F. Vervelidou, V. Lesur, A. Morschhauser, M. Grott, P. Thomas

► To cite this version:

F. Vervelidou, V. Lesur, A. Morschhauser, M. Grott, P. Thomas. On the accuracy of palaeopole estimations from magnetic field measurements. Geophysical Journal International, 2017, 211, pp.1669-1678. 10.1093/gji/ggx400. insu-03596869

HAL Id: insu-03596869

https://insu.hal.science/insu-03596869

Submitted on 4 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the accuracy of palaeopole estimations from magnetic field measurements

F. Vervelidou,<fnref n="1" /> V. Lesur,<fnref n="2" /> A. Morschhauser,<fnref n="1" /> M. Grott<fnref n="3" /> and P. Thomas<fnref n="3" />

Accepted 2017 September 22. Received 2017 September 14; in original form 2017 January 11

## SUMMARY

Various techniques have been proposed for palaeopole position estimation based on magnetic field measurements. Such estimates can offer insights into the rotational dynamics and the dynamo history of moons and terrestrial planets carrying a crustal magnetic field. Motivated by discrepancies in the estimated palaeopole positions among various studies regarding the Moon and Mars, we examine the limitations of magnetic field measurements as source of information for palaeopole position studies. It is already known that magnetic field measurements cannot constrain the null space of the magnetization nor its full spectral content. However, the extent to which these limitations affect palaeopole estimates has not been previously investigated in a systematic way. In this study, by means of the vector Spherical Harmonics formalism, we show that inferring palaeopole positions from magnetic field measurements necessarily introduces, explicitly or implicitly, assumptions about both the null space and the full spectral content of the magnetization. Moreover, we demonstrate through synthetic tests that if these assumptions are inaccurate, then the resulting palaeopole position estimates are wrong. Based on this finding, we make suggestions that can allow future palaeopole studies to be conducted in a more constructive way.

**Key words:** Geomagnetic excursions; Magnetic anomalies: modelling and interpretation; Palaeomagnetism; Reversals: process, time scale, magnetostratigraphy; Satellite magnetics.

## 1 INTRODUCTION

A magnetization distribution that has been induced by a centred dipole carries the information of the location of the poles of the inducing dipole, which is usually known as the palaeopole position. This fact is frequently used in terrestrial palaeomagnetic studies, where the direction of a rock sample's magnetization is estimated through laboratory measurements. The same idea has also been applied to estimate palaeopole positions for seamounts on Earth as well as palaeopole positions on Moon and Mars. In these cases, only a limited number of studies deal directly with rock samples (see, e.g. Cournède et al. 2012), while most studies estimate magnetization directions from magnetic field measurements. For this purpose, a variety of techniques have been suggested and used. The most common technique consists in selecting isolated magnetic field anomalies and then fitting signals generated by sources of simple form (see, e.g. Hood & Zakharian 2001; Arkani-Hamed 2001; Frawley & Taylor 2004; Boutin & Arkani-Hamed 2006; Hood et al. 2007; Quesnel et al. 2007; Berguig et al. 2008; Takahashi et al. 2014; Plattner & Simons 2015; Tsunakawa et al. 2015). A different technique consists in fitting the crustal magnetic field to a grid of equivalent dipolar sources (see, e.g. Langlais et al. 2004; Langlais & Purucker 2007; Milbury et al. 2012). Another method followed by

Sprenke & Baker (2000) and Milbury & Schubert (2010) considers a crust which has been magnetized by one single, centred inducing dipole. Alternative approaches include the method of Parker (1991) (see also Oliveira & Wieczorek 2017; Thomas et al. 2017), which is based on the assumption that the magnetization is of uniform direction, the statistical method of Sprenke (2005), and the method of Ditty & Ravat (2014, 2015), who determine palaeopoles from edge effect magnetic field anomalies.

The various studies have led to a large variety of palaeopole estimations, both for the Moon and Mars, and no consensus concerning the results has been reached [for Mars see, e.g., table A1 of Milbury et al. (2012), and for the Moon compare, e.g., fig. 6 of Arkani-Hamed & Boutin (2014) to fig. 3 of Takahashi et al. (2014) and fig. 6 of Oliveira & Wieczorek (2017)]. These discrepancies have been explained in terms of significant polar wander, but also in terms of noisy data, lack of high-resolution measurements (e.g. Quesnel et al. 2007), and differences in data processing and modeling techniques (e.g. Boutin & Arkani-Hamed 2006). Additionally, a good fit to the magnetic field data is not sufficient to ensure the accuracy of the obtained palaeopole locations because inverting magnetic field measurements for a magnetization distribution suffers from inherent non-uniqueness (e.g. Runcorn 1975; Maus & Haak 2003). The studies of Biswas & Ravat (2005a,b) address this problem. They

<footnote n="1">¹Helmholtz Centre Potsdam—GFZ German Research Centre for Geosciences, D-14473 Potsdam, Germany. E-mail: foteini@gfz-potsdam.de</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, F-75005 Paris, France</footnote>

<footnote n="3">³Institute of Planetary Research, German Aerospace Center, D-12489 Berlin, Germany</footnote>
show that different sources, resulting in different palaeopole locations, can explain equally well a given magnetic field pattern. In order to tackle the problem of non-uniqueness, additional assumptions about the sources are needed (see, e.g. Milbury et al. 2012, who use gravity data to constrain the location of the magnetized sources). Overall, many studies point out that the accuracy of the inferred palaeopoles strongly depends on the assumptions made about the underlying sources (e.g. Parker 1991; Sprenke & Baker 2000; Arkani-Hamed 2001; Sprenke 2005; Boutin & Arkani-Hamed 2006; Hood et al. 2007; Milbury et al. 2012; Tsunakawa et al. 2015), but the effect of these assumptions on palaeopole estimates has not been systematically investigated in the current literature.

In this study, we are interested in clarifying in which way and to what extent the assumptions made about the sources affect the inferred palaeopoles. We do so by conducting two synthetic tests. In the first test, we demonstrate what is the effect of introducing wrong assumptions about the null space of the magnetization. In the second test, we demonstrate what is the effect of introducing wrong assumptions about the spatial wavelengths of the magnetization that are not constrained by magnetic field measurements. The paper is divided into four parts: the first part presents the methodology employed to produce synthetic magnetic field data from a given magnetization distribution and then to infer palaeopole locations from these synthetic data. The second part presents the results of the two synthetic tests. A discussion about their implications on palaeopole studies based on real magnetic field data is presented in the third part of the paper. The conclusions of our study together with suggestions for future directions are presented in the fourth part of the paper. The expression of the vector Spherical Harmonics (SH) and the derivation of the key mathematical formula of the methodology section are given in Appendix.

## 2 METHODOLOGY

The basic principle of the two synthetic tests is the comparison of the true palaeopole location with the palaeopole locations recovered through an analysis of the induced magnetic field. In both synthetic tests, we generate a magnetization distribution, **M**, hereinafter termed the input or the induced magnetization. The induced magnetic field, **B**, observed at point **r** outside of the sources, is then calculated by means of the following formula (see, e.g. Blakely 1996, eq. 5.3)

$$ \mathbf{B}(\mathbf{r}) = -\nabla \left[ \frac{\mu_0}{4\pi} \int_{\Omega'} \mathbf{M}(\mathbf{r}') \nabla' \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} \right) d\Omega' \right], \quad (1) $$

with **r** = (r, θ, φ) the position vector of the observation point and r, θ and φ the radius, colatitude and longitude, respectively, in the sphere's reference system, **r**' = (r', θ', φ') the position vector of the magnetized point, dΩ' = (r')² sinθ'dθ'dφ'dr' and μ₀ the magnetic permeability of vacuum. We simulate the retrieval of palaeopole locations from the observed induced magnetic field by using **B** to infer the synthetic palaeopole locations.

### 2.1 First synthetic test: the effect of assumptions on magnetization's null space

#### 2.1.1 Generating the induced magnetic field

In the first synthetic test, the input magnetization, **M**, is induced by a centred dipole that magnetizes a susceptibility distribution, χ,

and is, therefore, given by (see, e.g. eqs 1.15–1.17 of Langel & Hinze 1998)

$$ \mathbf{M}(\mathbf{r}') \approx \frac{\chi(\mathbf{r}')}{\mu_0} \mathbf{B}_{\text{inducing}}(\mathbf{r}'), \quad (2) $$

where **B**_{inducing} is given by (see, e.g. Blakely 1996, eq. 4.14)

$$ \mathbf{B}_{\text{inducing}}(\mathbf{r}') = \frac{\mu_0}{4\pi} \frac{m}{r'^3} [3(\hat{\mathbf{m}} \cdot \hat{\mathbf{r}}')\hat{\mathbf{r}}' - \hat{\mathbf{m}}], \quad r' \neq 0, \quad (3) $$

with **m** = m·**m** the magnetic dipole moment. Eq. (2) is an approximation in that it takes only the inducing field and not the true ambient field into account, but their difference is negligible. The susceptibility distribution, χ, is considered to be confined to a spherical shell of infinitesimal thickness, on the surface of which it varies laterally according to:

$$ \chi(\mathbf{r}') = \chi_0 e^{-\left(\frac{|\mathbf{r}' - \mathbf{r}_c|}{\sigma}\right)^2}, \quad (4) $$

where |r' - r_c| denotes the arc length between r' and r_c, with r_c being the position vector of the maximum susceptibility, χ₀ and σ is a measure of the spread. This distribution has been chosen because it allows us to regulate the spectral content of **M** and **B**, since larger σ result in smoother susceptibility distributions.

Having generated **M**, we calculate **B** by means of eq. (1) on a Gauss-Legendre grid and obtain its SH coefficients up to a sufficiently large l_max, so that its full spectral content is represented. In the following, we present our methodology for obtaining palaeopole locations by means of the SH coefficients of **B**.

#### 2.1.2 Estimating palaeopole locations

In order to study the retrieval of palaeopole locations from measurements of the induced magnetic field **B**, we need to properly constrain what part of the magnetization can be recovered from such measurements. Thanks to the decomposition of the magnetization in vector SH, we can split **M** into its null space and into its part that gives rise to the observed magnetic field. This decomposition reads (see eq. 28 by Gubbins et al. 2011)

$$ \mathbf{M} = \mathcal{I} + \mathcal{E} + \mathcal{T}, \quad (5) $$

with

$$ \mathcal{I} = \sum_{m,l} g_{\mathcal{I},l}^{m} \mathbf{Y}_{l,l-1}^{m,c} + h_{\mathcal{I},l}^{m} \mathbf{Y}_{l,l-1}^{m,s}, \quad (6) $$

$$ \mathcal{E} = \sum_{m,l} g_{\mathcal{E},l}^{m} \mathbf{Y}_{l,l+1}^{m,c} + h_{\mathcal{E},l}^{m} \mathbf{Y}_{l,l+1}^{m,s}, \quad (7) $$

$$ \mathcal{T} = \sum_{m,l} g_{\mathcal{T},l}^{m} \mathbf{Y}_{l,l}^{m,c} + h_{\mathcal{T},l}^{m} \mathbf{Y}_{l,l}^{m,s}, \quad (8) $$

where {g_{\mathcal{T},l}^{m}, h_{\mathcal{T},l}^{m}}, {g_{\mathcal{E},l}^{m}, h_{\mathcal{E},l}^{m}} and {g_{\mathcal{T},l}^{m}, h_{\mathcal{T},l}^{m}} are the vector SH coefficients of order m and degree l of I, E and T, respectively, and {Y_{l,l-1}^{m,c}, Y_{l,l-1}^{m,s}}, {Y_{l,l+1}^{m,c}, Y_{l,l+1}^{m,s}} and {Y_{l,l}^{m,c}, Y_{l,l}^{m,s}} are the three types of real vector SH (their expressions are given in Appendix A, eqs A1–A3).

Vector SH have been used to distinguish between magnetic fields of internal and external origin (see, e.g. Backus 1996; Mayer & Maier 2006; Gerhards 2011). The study by Gubbins et al. (2011) introduces this formalism in the study of magnetization. By construction, I is the only part of **M** that generates an observable magnetic field outside of the magnetized layer, E describes the part of **M** that generates a magnetic field that is confined inside or below the magnetized layer and T is the toroidal part of **M**, which
is zero in the electrically insulating regions outside the magnetized layer (see Gubbins et al. 2011). Since $\mathcal{I}$ is the only part of $\mathbf{M}$ that gives rise to the magnetic field outside of the magnetized layer, it is hereinafter termed the visible magnetization. The sum $\mathcal{E} + \mathcal{T}$ spans the null space of the magnetization and is therefore not recoverable through magnetic field measurements without further assumptions (see, e.g. Gerhards 2015, for such a study).

In order to simulate the situation of imposing arbitrary constraints on the null space of the magnetization, we set $\mathcal{E} + \mathcal{T} = 0$ and calculate the palaeopoles based only on $\mathcal{I}$. For the geometry of the source chosen here, that is, infinitesimally thin spherical shell with laterally varying magnetization of finite spectral content, it can be shown that the vector SH coefficients of $\mathcal{I}$ are given by (see Appendix B for the derivation and eqs 32 and 33 by Gubbins et al. 2011, for a similar expression)

$$ \left\{ \begin{matrix} g_{\mathcal{I},l}^m \\ h_{\mathcal{I},l}^m \end{matrix} \right\} = \frac{R}{\mu_0} \sqrt{\frac{2l+1}{l}} \left( \frac{R}{r''} \right)^{l+1} \left\{ \begin{matrix} g_l^m \\ h_l^m \end{matrix} \right\}, \quad (9) $$

where $\{g_{\mathcal{I},l}^m, h_l^m\}$ are the SH coefficients of the induced magnetic field in units of nT, $R$ is their reference radius, $r''$ is the radius of the magnetized spherical shell and $\{g_{\mathcal{I},l}^m, h_{\mathcal{I},l}^m\}$ are the vector SH coefficients of $\mathcal{I}$ in units of nA. By using these coefficients, we calculate $\mathcal{I}$ on a given grid and obtain a palaeopole location for each point of this grid by means of the standard formulae used in palaeomagnetism (see, e.g. Butler 1992, chap. 7). We note here that we report the south pole of the inducing dipole as the palaeopole.

For the numerical implementation of this test case, we arbitrarily set, without loss of generality, the inducing, centred dipole to point to the North Pole and to have a magnetic moment $m$ equal to $7.94 \times 10^{22} \text{ Am}^2$ (close to the current dipole moment of the Earth, see Constable (2007). We set the radius of the spherical shell equal to Mars' mean radius, that is, $R = 3393.5 \text{ km}$ and take $r'' = R - 30 \text{ km}$. We take the maximum of the susceptibility distribution, $\chi_0$, to be equal to 1 and to be located at $(\theta_c, \phi_c) = (50^\circ, 150^\circ)$ (see Fig. 1a), and we set $\sigma$ equal to 1000 km. Finally, we choose the SH truncation degree $l_{max} = 30$, which is found to be compatible with the choice of $\sigma = 1000 \text{ km}$. This again is without loss of generality. What is of importance here is that $l_{max}$ is large enough to represent fully the spectral content of the magnetic signal generated by the magnetization distribution. Had we opted for a less smooth susceptibility distribution, we would have to increase $l_{max}$. The results from this synthetic test are shown in the third section, together with the results from the second synthetic test.

## 2.2 Second synthetic test: the effect of measurements' limited spatial resolution

### 2.2.1 Generating the induced magnetic field

For our second synthetic test, we consider a dipolar magnetization source:

$$ \mathbf{M}(\mathbf{r}') = \delta_{r',r''} \delta_{\theta',\theta''} \delta_{\phi',\phi''}, \quad (10) $$

where $\mathbf{r}' = (r'', \theta'', \phi'')$ is the position vector of the dipole, and $\delta$ is Kronecker's delta. We opt for this source because it gives rise to a magnetic field of very high spectral content and it has zero null space, since according to eq. (3) there exists no other magnetization $\mathbf{m}' \neq \mathbf{m}$ that yields the same magnetic field as $\mathbf{m}$. Similarly to the previous synthetic test, we introduce $\mathbf{M}$ into eq. (1), calculate $\mathbf{B}$ on a Gauss-Legendre grid and obtain its SH coefficients up to a certain $l_{max}$.

### 2.2.2 Estimating palaeopole locations

This time, the part of the magnetization that can be recovered from measurements of the induced magnetic field, $\mathbf{B}$, is the part whose spectral content is constrained by the measurements. Splitting the magnetization into the constrained and unconstrained parts, gives

$$ \mathbf{M} = \mathcal{I} + \mathcal{I}', \quad (11) $$

where

$$ \mathcal{I} = \sum_{l=1}^{l_{max}} \sum_{m} g_{\mathcal{I},l}^{m} \mathbf{Y}_{l,l-1}^{m,c} + h_{\mathcal{I},l}^{m} \mathbf{Y}_{l,l-1}^{m,s}, \quad (12) $$

and

$$ \mathcal{I}' = \sum_{l=l_{max}+1}^{\infty} \sum_{m} g_{\mathcal{I},l}^{m} \mathbf{Y}_{l,l-1}^{m,c} + h_{\mathcal{I},l}^{m} \mathbf{Y}_{l,l-1}^{m,s}. \quad (13) $$

In order to examine the effect of making arbitrary assumptions about the part of the spectral content that is not constrained by magnetic field measurements, we infer palaeopoles based only on a truncated version of the magnetization source. We do so by neglecting $\mathcal{I}'$, and we infer palaeopoles based only on $\mathcal{I}$. As previously, we obtain the vector SH coefficients of $\mathcal{I}$ by means of eq. (9), we calculate $\mathcal{I}$ on a given grid, and we obtain palaeopoles for each point of this grid.

For the numerical implementation of this test case, we arbitrarily set, again without loss of generality, $l_{max} = 110$, $\mathbf{m}_d = (10^{12}, 10^{12}, 10^{12}) \text{ Am}^2$ and $\mathbf{r}' = (R - d, 125^\circ, 295^\circ)$, with $R = 3393.5 \text{ km}$ and $d = 30 \text{ km}$.

## 3 RESULTS

The susceptibility distribution used in the first synthetic test is shown in Fig. 1(a). Stereographic projections of the North Pole (latitude between degrees $60^\circ$ and $90^\circ$) and South Pole (latitude between degrees $-60^\circ$ and $-90^\circ$) are shown, to the left and right, respectively, along with a Robinson projection of all latitudes in the centre of Fig. 1(a). Note that latitude is shown here, while coordinates in the text refer to colatitude. The corresponding induced magnetization, as calculated from eq. (2), is shown in Fig. 1(b). As expected, its East component is zero. For illustration purposes, we show in Fig. 1(c) the induced magnetic field at $R + 300 \text{ km}$. As already mentioned, we base our analysis on the SH coefficients of this induced magnetic field, which are independent of the observational altitude, since we deal with noise-free and perfectly distributed synthetic data. The visible magnetization, obtained through eqs (6) and (9), is shown in Fig. 1(d). Its vector components differ significantly from those of the input magnetization, both in terms of intensity and geometry. Still, the visible magnetization exactly reproduces the observed magnetic field (shown in Fig. 1c). Finally, the real and the estimated palaeopoles are shown in Fig. 1(e). Again, stereographic projections of the North Pole (latitude between degrees $60^\circ$ and $90^\circ$, and isolines every $10^\circ$) and South Pole (latitude between degrees $-60^\circ$ and $-90^\circ$, and isolines every $10^\circ$) are shown, in the left- and right-hand panels, respectively, along with a Robinson projection of all latitudes in the centre of the figure. The red circle corresponds to the palaeopoles estimated from all the gridpoints of the input magnetization. As expected, they all lie at the South Pole, which is the location of the true palaeopole. The green triangle corresponds to the palaeopole derived from the gridpoint of the visible magnetization that lies at the susceptibility maximum. We observe that it does not coincide with the real palaeopole location. Moreover, we plot the palaeopoles derived from the gridpoints of
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The results of the first synthetic test. (a) The susceptibility distribution. (b)–(d) The magnetization induced by a centred dipole that points to the geographic North Pole, the induced magnetic field at 300 km altitude, and the corresponding visible magnetization. (e) The spatial distribution of palaeopoles derived from the input magnetization (red circle), from the point of the visible magnetization that lies at the maximum of the susceptibility (green triangle), and from the points of the visible magnetization that lie within a circle of 1000 km around the maximum of the susceptibility (blue crosses). Stars are the palaeopoles recovered by means of the method by Parker, 1991 (see Section 5 for details). See the text for details on the projections.</figcaption></figure>

the visible magnetization that lie within a radius of 1000 km around
the susceptibility maximum (blue crosses). We observe that they
are all different from the correct palaeopole position and that they
show a large dispersion. This demonstrates that the estimated mag-

netization distribution, under the assumption of zero null space,
does not correspond any more to a magnetization that has been
induced by a single inducing dipole, as it is the case for the input
magnetization.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a)–(c) The magnetic moment of a single dipole, the magnetic field induced by this dipole truncated at *l*<sub>max</sub> = 110 and evaluated at 300 km altitude, and the corresponding visible magnetization. (d) The palaeopole location derived from this dipole (red circle), from the point of the visible magnetization that lies at the location of the dipole (green triangle), and from the palaeopole locations derived from the points of the visible magnetization that lie within a circle of 200 km around the location of the dipole (blue crosses). See the text for details on the projections.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** (a) The magnetic field at Mars mean radius generated by a dipole located at ($\theta = 7°$ and $\phi = 32°$) and at 150 km depth with $I = 70°$, $D = -10°$ and $m = 5 \times 10^{14}$ Am$^2$, over a spherical cap with $\theta_0 = 18°$ centred at the North Pole. (b) The lithospheric magnetic field of Mars at the same location as predicted by the model of Morschhauser *et al.* (2014). For both (a) and (b), projection is Stereographic (latitude between degrees 60° and 90°). (c) The SH power spectrum of the first 300 degrees of the magnetic field shown in the first row (blue curve, where the bracket shows its part that is not yet constrained by the available magnetic field models), and the R-SCHA power spectrum of the lithospheric magnetic field shown in the second row (red curve).</figcaption></figure>

can be uniquely recovered. We suggest that a fruitful avenue would be to focus on finding such cases that are at the same time relevant to the problem of palaeopoles estimation and on applying them on real magnetic field measurements (see Gerhards 2016) for a study in this direction using synthetic data). Eqs (5)–(8) provide a framework to search for conditions under which $\mathcal{I}$ can be inverted for **M**. However, we expect any geophysically plausible assumptions relevant to palaeopole studies to hold only over local and regional spatial scales. For this reason, eqs (5)–(8) should be addressed in a localized reference system. Such calculations go beyond the scope of this paper but are direct perspectives of this study.

A method based on a geophysically plausible assumption relevant to palaeopole studies, the assumption of unidirectional magnetization over local spatial scales, has been introduced by Parker (1991)

and has been recently applied by Thomas *et al.* (2017) on Mars and Oliveira & Wieczorek (2017) on Moon. We test their method within the framework of our first synthetic test, described in Section 2.1, since we expect that the assumption of unidirectional magnetization is valid there if we focus on a sufficiently small area over the sphere. For this purpose, we invert all three components of the induced magnetic field estimated at 300 km altitude, for a grid of 931 dipoles covering a circular area of 4°, 2°, 1° and 0.5° radius, centred on the maximum of the susceptibility distribution [for more details on this method see Parker (1991); Thomas *et al.* (2017); Oliveira & Wieczorek (2017)]. Fig. 1(e) shows the estimated best-fit palaeopoles for all four cases, each estimated palaeopole is shown by one star. According to this figure, for an area of 0.5° radius, the palaeopole is recovered correctly (see green star with black outline, almost
overlapping with the red circle). This is because on such a small area the magnetization direction is indeed almost uniform. However, as the area gets larger, the estimated palaeopoles move away from the correct palaeopole location (see yellow, black and red stars, all with black outlines, for the areas of 1°, 2° and 4° radius, respectively). This demonstrates the main conclusion of our study: palaeopole location estimation is correct if and only if the assumptions about the sources, whether implicit or explicit, are accurate.

## ACKNOWLEDGEMENTS

We wish to thank two anonymous reviewers for their constructive comments that significantly improved our paper. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), within the Schwerpunktprogramm 1488-Planetary Magnetism under grants LE2477/3-2 (FV and VL) and GR3751/1-1 (AM, MG and PT). It is based on synthetic data that can be reproduced by means of the equations given in the paper and on the lithospheric magnetic field model of Mars by Morschhauser et al. (2014). This is IPGP publication number: 3887.

## REFERENCES

* Arkani-Hamed, J., 2001. Paleomagnetic pole positions and pole reversals of Mars, *Geophys. Res. Lett.*, **28**(17), 3409–3412.

* Arkani-Hamed, J. & Boutin, D., 2014. Analysis of isolated magnetic anomalies and magnetic signatures of impact craters: evidence for a core dynamo in the early history of the moon, *Icarus*, **237**, 262–277.

* Backus, G., Parker, R. & Constable, C., 1996. *Foundations of Geomagnetism*, Cambridge Univ. Press.

* Baratchart, L., Hardin, D.P., Lima, E.A., Saff, E.B. & Weiss, B.P., 2012. Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions, *Inverse Probl.*, **29**, 015004, doi:10.1088/0266-5611/29/1/015004.

* Berguig, M.C., Hamoudi, M., Cohen, Y. & Thébault, E., 2008. Lunar magnetic pole positions deduced from high albedo magnetic anomalies, *Asian J. Earth Sci.*, **1**(2), 62–76.

* Biswas, S. & Ravat, D., 2005a. Why Meaningful paleopoles can't be determined without special assumptions from Mars global surveyor data?, in *36th Annual Lunar and Planetary Science Conference*, p. 2192.

* Biswas, S. & Ravat, D., 2005b. Non-uniqueness of the modeled magnetization vectors used in determining paleopoles on Mars, *Master thesis*, Southern Illinois University Carbonale.

* Blakely, R.J., 1996. *Potential Theory in Gravity and Magnetic Applications*, Cambridge Univ. Press.

* Boutin, D. & Arkani-Hamed, J., 2006. Pole wandering of Mars: evidence from paleomagnetic poles, *Icarus*, **181**(1), 13–25.

* Butler, R.F., 1992. *Paleomagnetism: Magnetic Domains to Geologic Teranes*, Blackwell Scientific Publications.

* Constable, C., 2007. Dipole moment variation, in *Encyclopedia of Geomagnetism and Paleomagnetism*, pp. 159–161, Springer, The Netherlands.

* Cournède, C., Gattacceca, J. & Rochette, P., 2012. Magnetic study of large Apollo samples: possible evidence for an ancient centered dipolar field on the Moon, *Earth planet. Sci. Lett.*, **331**, 31–42.

* Ditty, M.L. & Ravat, D., 2014. Mars paleopole and magnetization estimates from magnetization edge effects: a new method, in *45th Annual Lunar and Planetary Science Conference*, p. 2075.

* Ditty, M.L. & Ravat, D., 2015. Finding Mars paleopoles from magnetization edge effects to determine the history of Mars' core dynamo, *Master's thesis*, University of Kentucky.

* Frawley, J.J. & Taylor, P.T., 2004. Paleo-pole positions from martian magnetic anomaly data, *Icarus*, **172**(2), 316–327.

* Gerhards, C., 2011. Spherical decompositions in a global and local framework: theory and an application to geomagnetic modeling, *Int. J. Geomath.*, **1**(2), 205–256.

* Gerhards, C., 2015. On the unique reconstruction of induced spherical magnetizations, *Inverse Prob.*, **32**(1), doi:10.1088/0266-5611/32/1/015002.

* Gerhards, C., 2016. Improved reconstruction of dipole directions from spherical magnetic field measurements, preprint (arXiv:1608.01216).

* Gubbins, D., Ivers, D., Masterton, S.M. & Winch, D.E., 2011. Analysis of lithospheric magnetization in vector spherical harmonics, *Geophys. J. Int.*, **187**, 99–117.

* Hood, L.L. & Zakharian, A., 2001. Mapping and modeling of magnetic anomalies in the northern polar region of mars, *J. geophys. Res.*, **106**(E7), 14601–14619.

* Hood, L.L., Richmond, N.C., Harrison, R.J. & Lillis, K.P., 2007. East-west trending magnetic anomalies in the Southern Hemisphere of Mars: modeling analysis and interpretation, *Icarus*, **191**(1), 113–131.

* Langel, R.A. & Hinze, W.J., 1998. *The Magnetic Field of the Earth's Lithosphere: The Satellite Perspective*, Cambridge Univ. Press.

* Langlais, B. & Purucker, M., 2007. A polar magnetic paleopole associated with Apollinaris Patera, Mars, *Planet. Space Sci.*, **55**(3), 270–279.

* Langlais, B., Purucker, M.E. & Mandea, M., 2004. Crustal magnetic field of Mars, *J. geophys. Res.*, **109**, E02008, doi:10.1029/2003JE002048.

* Maus, S. & Haak, V., 2003. Magnetic field annihilators: invisible magnetization at the magnetic equator, *Geophys. J. Int.*, **155**(2), 509–513.

* Mayer, C. & Maier, T., 2006. Separating inner and outer Earth's magnetic field from CHAMP satellite measurements by means of vector scaling functions and wavelets, *Geophys. J. Int.*, **167**(3), 1188–1203.

* Milbury, C. & Schubert, G., 2010. Search for the global signature of the martian dynamo, *J. geophys. Res.*, **115**, E10010, doi:10.1029/2010JE003617.

* Milbury, C., Schubert, G., Raymond, C.A., Smrekar, S.E. & Langlais, B., 2012. The history of Mars' dynamo as revealed by modeling magnetic anomalies near Tyrhrenus Mons and Syrtis Major, *J. geophys. Res.*, **117**, E1007, doi:10.1029/2012JE004099.

* Morschhauser, A., Lesur, V. & Grott, M., 2014. A spherical harmonic model of the lithospheric magnetic field of Mars, *J. geophys. Res.*, **119**(6), 1162–1188.

* Oliveira, J.S. & Wieczorek, M.A., 2017. Testing the axial dipole hypothesis for the moon by modeling the direction of crustal magnetization, *J. geophys. Res.*, **122**(2), 383–399.

* Parker, R.L., 1991. A theory of ideal bodies for seamount magnetism, *J. geophys. Res.*, **96**(B10), 16101–16112.

* Plattner, A. & Simons, F.J., 2015. High-resolution local magnetic field models for the martian south pole from mars global surveyor data, *J. geophys. Res.*, **120**(9), 1543–1566.

* Quesnel, Y., Langlais, B. & Sotin, C., 2007. Local inversion of magnetic anomalies: implication for Mars' crustal evolution, *Planet. Space Sci.*, **55**(3), 258–269.

* Runcorn, S.K., 1975. An ancient lunar magnetic dipole field, *Nature*, **253**, 701–703.

* Spenke, K.F., 2005. Martian magnetic paleopoles: a geostatistical approach, *Geophys. Res. Lett.*, **32**, L09201, doi:10.1029/2005GL022840.

* Spenke, K.F. & Baker, L.L., 2000. Magnetization, paleomagnetic poles, and polar wander on Mars, *Icarus*, **147**(1), 26–34.

* Takahashi, F., Tsunakawa, H., Shimizu, H., Shibuya, H. & Matsushima, M., 2014. Reorientation of the early lunar pole, *Nature Geosci.*, **7**, 409–412.

* Thébault, E., Schott, J.J. & Mandea, M., 2006. Revised spherical cap harmonic analysis (R-SCHA): validation and properties, *J. geophys. Res.*, **111**, B01102, doi:10.1029/2005JB003836.

* Thomas, P., Grott, M., Morschhauser, A. & Vervelidou, F., 2017. Paleopole reconstruction of Martian magnetic field anomalies, in *48th Annual Lunar and Planetary Science Conference*.

* Tsunakawa, H., Takahashi, F., Shimizu, H., Shibuya, H. & Matsushima, M., 2015. Surface vector mapping of magnetic anomalies over the Moon using Kaguya and Lunar Prospector observations, *J. geophys. Res.*, **120**(6), 1160–1185.

* Vervelidou, F. & Thébault, E., 2015. Global maps of the magnetic thickness and magnetization of the Earth's lithosphere, *Earth Planets Space*, **67**(1), 1–19, doi:10.1186/s40623-015-0329-5.
APPENDIX A: REAL VECTOR SPHERICAL HARMONICS

The expressions of the real vector Spherical Harmonics (SH) employed in this study are (see eqs 16–27 of Gubbins et al. 2011, for the respective expressions in terms of complex vector SH)

$$
\mathbf{Y}_{l,l-1}^{m(c,s)} = \frac{1}{r^{l-1}\sqrt{l(2l+1)}} \nabla[r^l \mathbf{Y}_l^{m(c,s)}(\theta, \phi)]. \quad (A1)
$$

$$
\mathbf{Y}_{l,l+1}^{m(c,s)} = \frac{r^{l+2}}{\sqrt{(l+1)(2l+1)}} \nabla \left[ \frac{1}{r^{l+1}} \mathbf{Y}_l^{m(c,s)}(\theta, \phi) \right] \quad (A2)
$$

and

$$
\mathbf{Y}_{l,l}^{m(c,s)} = - \frac{1}{\sqrt{l(l+1)}} \mathbf{r} \times \nabla \mathbf{Y}_{l}^{m(c,s)}(\theta, \phi), \quad (A3)
$$

with

$$
Y_l^{m,c} = P_l^m(\cos \theta) \cos(m\phi), \quad (A4)
$$

and

$$
Y_l^{m,s} = P_l^m(\cos \theta) \sin(m\phi) \quad (A5)
$$

being the Schmidt semi-normalized real SH of degree *l* and order *m*, whose norm is given by

$$
\oint \mathbf{Y}_l^{m(c,s)} \mathbf{Y}_{l'}^{m'(c,s)} d\Omega = \frac{4\pi}{2l+1} \delta_{ll'} \delta_{mm'}. \quad (A6)
$$

APPENDIX B: CONVERSION FORMULA

In this section, we derive a formula that converts any set of magnetic field SH coefficients to a set of visible magnetization SH coefficients,
under the assumption that the magnetization is confined to a spherical layer of infinitesimal thickness (see Gubbins *et al.* 2011, for an
alternative derivation that leads to a similar formula).

We denote the layer's radius by $r'' = R - d$, where $R$ is the reference radius of the magnetic field SH coefficients. Then, **M** reads:

$$
\mathbf{M}(\mathbf{r}') = \delta_{r'r''} \mathbf{g}(\theta', \phi'), \quad (B1)
$$

where $\delta$ is Kronecker's delta and $\mathbf{g}$ represents the laterally varying part of $\mathbf{M}$. We have (e.g. Blakely 1996, eq. 5.2)

$$
V(\mathbf{r}) = \frac{\mu_0}{4\pi} \int_{\Omega'} \mathbf{M}(\mathbf{r}') \nabla' \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} \right) d\Omega', \quad (B2)
$$

with $V$ the magnetic potential at point $\mathbf{r}$ due to the magnetization distribution $\mathbf{M}$ that lies in the magnetized volume $\Omega'$. Introducing eq. (B1) into eq. (B2), we obtain

$$
V(\mathbf{r}) = \frac{\mu_0}{4\pi} \int_{\partial\Omega_{r''}} \mathbf{g}(\theta', \phi') \nabla' \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} \right) \Big|_{\mathbf{r}'=\mathbf{r}''} ds, \quad (B3)
$$

with $ds = (r'')^2 \sin\theta' d\theta' d\phi'$, and $\vartheta \Omega_{r''}$ the surface of volume $\Omega$ for which $r = r''$. According to eq. (5), $\mathbf{g}$ can be split into its $\mathcal{I}$, $\mathcal{E}$ and $\mathcal{T}$ parts. This gives

$$
V(\mathbf{r}) = \frac{\mu_0}{4\pi} \int_{\partial\Omega_{r''}} \mathcal{I}(\theta', \phi') \nabla' \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} \right) \Big|_{\mathbf{r}'=\mathbf{r}''} ds \quad (B4)
$$

$$
+ \frac{\mu_0}{4\pi} \int_{\partial \Omega_{r''}} \mathcal{E}(\theta', \phi') \nabla' \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} \right) \Big|_{\mathbf{r}'=\mathbf{r}''} ds \quad (B5)
$$

$$
+ \frac{\mu_0}{4\pi} \int_{\partial \Omega_{r''}} \mathcal{T}(\theta', \phi') \nabla' \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} \right) \Big|_{\mathbf{r}'=\mathbf{r}''} ds. \quad (B6)
$$

Using eqs (6)–(8) and (A1)–(A3), and replacing the term $\nabla' (\frac{1}{|\mathbf{r}-\mathbf{r}'|})|_{\mathbf{r}'=\mathbf{r}''}$ by its SH expansion, we find that the terms (B5) and (B6) are zero,
as expected. For the remaining term, (B4), we obtain:

$$
\begin{align*}
V(\mathbf{r}) &= \frac{\mu_0}{4\pi} \int_{\partial\Omega_{r''}} \mathcal{I}(\theta', \phi') \nabla' \left( \frac{1}{|\mathbf{r} - \mathbf{r}'|} \right) \Big|_{\mathbf{r}'=\mathbf{r}''} ds \\
&= \frac{\mu_0}{4\pi} \int_0^\pi \int_0^{2\pi} \sum_{l'} \sum_{m'} \left[ g_{\mathcal{I},l'}^{m',c}(\theta', \phi') + h_{\mathcal{I},l'}^{m',s}(\theta', \phi') \right] \left( \sqrt{\frac{l'}{2l' + 1}} \right) ds d\Omega_{r''}.
\end{align*}
$$
$$
\begin{aligned}
& \cdot \frac{1}{r^2} \sum_l \sum_m l \left(\frac{r''}{r}\right)^{l-1} \left[Y_l^{m,c}(\theta, \phi) Y_l^{m,c}(\theta', \phi') + Y_l^{m,s}(\theta, \phi) Y_l^{m,s}(\theta', \phi')\right] ds \\
& + \frac{\mu_0}{4\pi} \int_0^\pi \int_0^{2\pi} \sum_{l'} \sum_{m'} \left[g_{l,l'}^{m'} \frac{\partial Y_{l'}^{m',c}(\theta', \phi')}{\partial \theta'} + h_{l,l'}^{m'} \frac{\partial Y_{l'}^{m',s}(\theta', \phi')}{\partial \theta'}\right] \left(\sqrt{\frac{1}{l'(2l' + 1)}}\right) ds \\
& \cdot \frac{1}{rr''} \sum_l \sum_m \left(\frac{r''}{r}\right)^l \left[Y_l^{m,c}(\theta, \phi) \frac{\partial Y_l^{m,c}(\theta', \phi')}{\partial \theta'} + Y_l^{m,s}(\theta, \phi) \frac{\partial Y_l^{m,s}(\theta', \phi')}{\partial \theta'}\right] ds \\
& + \frac{\mu_0}{4\pi} \int_0^\pi \int_0^{2\pi} \sum_{l'} \sum_{m'} \left[g_{l,l'}^{m'} \frac{\partial Y_{l'}^{m',c}(\theta', \phi')}{\partial \phi'} + h_{l,l'}^{m'} \frac{\partial Y_{l'}^{m',s}(\theta', \phi')}{\partial \phi'}\right] \left(\sqrt{\frac{1}{l'(2l' + 1)}}\right) \frac{1}{\sin \theta'} ds \\
& \cdot \frac{1}{rr'' \sin \theta'} \sum_l \sum_m \left(\frac{r''}{r}\right)^l \left[Y_l^{m,c}(\theta, \phi) \frac{\partial Y_l^{m,c}(\theta', \phi')}{\partial \phi'} + Y_l^{m,s}(\theta, \phi) \frac{\partial Y_l^{m,s}(\theta', \phi')}{\partial \phi'}\right] ds.
\end{aligned}
$$

Further on, integrating over $\theta'$ and $\phi'$, while accounting for well-known properties of the SH, we obtain:

$$V(\mathbf{r}) = \mu_0 \sum_l \sum_m \left(\sqrt{\frac{l}{2l+1}}\right) \left(\frac{r''}{r}\right)^{l+1} \left(g_{l,l}^m Y_l^{m,c} + h_{l,l}^m Y_l^{m,s}\right). \qquad (B7)$$

The magnetic potential $V(\mathbf{r})$ can also be expanded in SH:

$$V(\mathbf{r}) = R \sum_l \sum_m \left(\frac{R}{r}\right)^{l+1} \left(g_l^m Y_l^{m,c} + h_l^m Y_l^{m,s}\right). \qquad (B8)$$

Combining eqs (B7) and (B8), we can establish a relationship between the coefficients {$g_l^m, h_l^m$} of $V$ and the coefficients {$g_{l,l}^m, h_{l,l}^m$} of $\mathcal{I}$:

$$\mu_0 \sum_l \sum_m \sqrt{\frac{l}{2l+1}} \left(\frac{r''}{r}\right)^{l+1} \left(g_{l,l}^m Y_l^{m,c} + h_{l,l}^m Y_l^{m,s}\right) = R \sum_l \sum_m \left(\frac{R}{r}\right)^{l+1} \left(g_l^m Y_l^{m,c} + h_l^m Y_l^{m,s}\right). \qquad (B9)$$

Multiplying both sides of eq. (B9) with the real SH, integrating over the surface of a sphere and accounting for the orthogonality properties of the SH functions, we obtain the following relationship:

$$\left\{ \begin{matrix} g_{l,l}^m \\ h_{l,l}^m \end{matrix} \right\} = \frac{R}{\mu_0} \sqrt{\frac{2l+1}{l}} \left(\frac{R}{r''}\right)^{l+1} \left\{ \begin{matrix} g_l^m \\ h_l^m \end{matrix} \right\}, \qquad (B10)$$

where $R$ is given in metres, {$g_l^m, h_l^m$} in units of nT and {$g_{l,l}^m, h_{l,l}^m$} in units of nA, being the coefficients of a vertically integrated magnetization.