Retrieving lithospheric magnetisation distribution from magnetic field models

SUMMARY We investigate to which extent the radially averaged magnetisation of the lithosphere can be recovered from the information content of a spherical harmonic model of the generated magnetic ﬁeld combined with few simple hypotheses. The results obtained show ﬁrstly that a hypothesis of magnetisation induced by a ﬁeld of internal origin, even over a localised area, is not suﬃcient to recover uniquely the radially averaged magnetisation and, secondly, that this magnetisation can be recovered when a constant magnetisation direction is assumed. An algorithm to recover the magnetisation direction and distribution is then described and tested over a synthetic example. It requires to introduce a cost function that vanishes when estimated in a system of coordinate with its Z axis aligned with the magnetisation direction. Failing to ﬁnd a vanishingly small value for the cost function is an indication that a constant magnetisation direction is not a valid hypothesis for the studied magnetic ﬁeld model. The range of magnetisation directions that are compatible with the magnetic ﬁeld model and a given noise level, can also be estimated. The whole process is illustrated by analysing a local, isolated maximum of the Martian magnetic ﬁeld.


INTRODUCTION
Rock magnetisation, particularly its remanent part, carries important information about a planet history. For the Earth it is used to study the evolution of its magnetic field, but also its tectonic and thermal history, the intrusion of dikes, volcanism ... etc... The list is not exhaustive and applications linked to exploration geophysics are not included. Unfortunately there is no method available to estimate uniquely the 3-D rock magnetisation distribution from measurements of the magnetic field it generates over a region of interest. Despite the wealth of possible applications, progresses are slow in this field of research. Studying the rock magnetisation therefore still relies on significant field work and collecting small samples that have to be studied later in laboratories. Inverting magnetic survey data to derive rock magnetisation remains a major challenge even when considering only the vertically integrated (or the radially averaged) magnetisation of a layer, drastically reducing this way the complexity of the problem.
It has been known for a long time that the distribution of magnetisation that can be recovered from magnetic measurements above a magnetised layer is non-unique. A typical example is the weakness of the magnetic field due to a spherical shell of uniform susceptibility magnetised by a potential field whose sources are inside the sphere. This problem has been studied first by Runcorn (1975) in the linear approximation to explain the weakness of the lunar magnetic field. It has then be studied by Jackson et al. (1999), and Lesur & Jackson (2000) for spheroidal shell and fully non-linear theory, respectively. Another important result has been obtained by Maus & Haak (2003). They show that even if the inducing field is as simple as a dipole inside the Earth, the radially integrated susceptibility is not uniquely defined when the associated magnetic field is known above the Earth. Numerous other results have been obtained in marine geosciences -e.g Parker et al. (1987). In particular it has been shown that a hypothesis of constant magnetisation direction, combined with imposing positiveness to the model parameters, allows to recover uniquely the direction of magnetisation (Parker 1991).
The technique has been applied not only to seamounts, but also to study magnetic signals from planet Mars (Thomas et al. 2018) and the Moon (Oliveira & Wieczorek 2017). Although the directions derived using this technique are in principle correct, there is no simple way to estimate the accuracy of the obtained direction, and no possibility to test the hypothesis of constant magnetisation direction. Recently, Gerhards (2016Gerhards ( , 2019 provided results on unicity but these results are not immediately applicable to geomagnetic studies. However, Gubbins et al. (2011) made a major step forward to solve the problem of finding the vertically in-Retrieving lithospheric magnetisation distribution 3 tegrated magnetisation on a sphere from spherical harmonic representation of the magnetic field it generates. They consider the case where the magnetisation can be decomposed in vector spherical harmonics, and identified the part of the magnetisation that contributes to the magnetic field (hereinafter the visible magnetisation), and the part that does not contribute (hereinafter the hidden magnetisation).
To progress further, it remains to define which simple and realistic hypotheses have to be set to be able to define in a unique way the whole magnetisation that is associated with the observed magnetic field. By whole magnetisation we refer here to the sum of visible and hidden magnetisation. We point out that the simplest hypothesis where the hidden magnetisation vanishes is, most of the time, a poor choice. This is because the visible magnetisation alone is not compatible with physically acceptable hypotheses such as assuming that the magnetisation is induced by an internal field (Vervelidou et al. 2017b). Although this latter hypothesis has often been used, we show in section 3 that it is not sufficient to define uniquely the whole magnetisation. This has already been demonstrated in Vervelidou & Lesur (2018). The results obtained in this work are repeated here for completeness. Two other cases are presented: -First, we consider in section 4 the case of a localised distribution of magnetisation induced by a field of internal origin and we show here that a localised susceptibility distribution cannot be uniquely estimated from surface magnetic field measurements assuming an internal dipolar inducing field. The non-uniqueness is however not as severe as for the global case.
-Second, in section 5, is presented the case where the radially averaged magnetisation direction is constant. We show that, under this hypothesis, both the direction of magnetisation and its distribution can be recovered from knowledge of the generated magnetic field. Furthermore, we set a criterion that allows to test the validity of the unique magnetisation direction assumption. We point out that a constant magnetisation direction is generally not an acceptable hypothesis unless working over a very small area. This case is therefore a limit case of localised distribution of internally induced magnetisation for a very small area. These latter results are applied for simple synthetic examples in section 6. Finally in section 7, we present and discuss results of magnetisation estimation for a localised area on Mars. Section 8 concludes this paper. We start however, in the next section, by shortly recalling results from Gubbins et al. (2011) that are presented here using a slightly different normalisation.

FIELD MODEL
We consider a simple model for the Earth's magnetised lithosphere made of a spherical shell of constant thickness d. Following Gubbins et al. (2011), the vector field of magnetisation M on a sphere of radius r = a is decomposed in three components: (2.1) that are defined by: Here (θ, φ) are two angles defining a point on a sphere of radius r. If the orthonormal system of coordinates XY Z is such that Z is aligned with the Earth rotation axis, and X points towards the Greenwich meridian, these angles are the usual colatitude and longitude. The Y m (θ, φ) are the Schmidt semi-normalized real spherical harmonics used in geomagnetism where negative orders, m < 0, are associated with sin(|m|φ) terms whereas null or positive orders, m ≥ 0, are associated with cos(|m|φ) terms. ∇ is the usual gradient operator, whereas ∇ h is its horizontal component. We use the normalisation: (2.6) The magnetic field generated in the lithosphere can also be described above the Earth surface with vector spherical harmonics. Our choice of normalisation gives: where theg m are the Gauss coefficients for the lithospheric field. This way of describing the magnetic field is strictly equivalent to the usual: whereṼ (r, θ, φ) is the magnetic potential for fields of internal origin.
For any volume V of magnetised rocks, the relation linking the magnetisation to the magnetic potential it generates is:Ṽ where H(r, r ) = 1/|r − r | is the usual Green function with r defining the position (r, θ, φ), and r the integration point (r , θ , φ ). The ∇ operator applies to the prime quantities. We are interested in calculating the magnetic potential generated above the Earth surface, so r > r and the Green function expressed in term of spherical harmonics is: (2.10) Therefore multiplying both sides of equation (2.9) by Y m (θ, φ) and integrating over the unit sphere, gives:g where a is the radius of the magnetised layer of thickness d. It is also the reference radius of the Gauss coefficients. With the orthogonality of the vector spherical harmonics, the surface integral over the unit sphere S 1 gives 4π e m . We note that the surface integral leaves out the contributions of i m and t m . Assuming that e m are independent of the radius, the radial integral becomes: which finally gives:g This relation (2.13) differs from Gubbins et al. (2011) not only by the normalisation, but also by the way the radial dependence is handled that is here the same as in Vervelidou et al. also implies that the visible magnetisation can be recovered only for the wavelengths for which the lithospheric magnetic field is known.

RECOVERING THE GLOBAL INDUCED MAGNETISATION
Let us assume that the lithospheric field, in equation (2.7), is defined from spherical harmonic degree = 1 up to degree L. The path we follow to recover the whole magnetisation is to derive first a susceptibility model from the visible magnetisation and, in a second step, estimate the hidden magnetisation from this model of susceptibility (Vervelidou & Lesur 2018).
We consider a known magnetic field of internal origin B i that induces a magnetisation M in the rocks. We assume further that these rocks have a scalar susceptibility χ that can be decomposed in a combination of spherical harmonics. Then the generated magnetisation is aligned with the inducing field and is given by: where the g m and the χ m are the spherical harmonic coefficients for the inducing field and the susceptibility respectively. These are defined relative to the same reference radius r = a as in section 2. Multiplying both sides byŶ m , −1 (θ, φ), and integrating over the spherical surface S 1 of radius unity leads to a relation linking the spherical harmonic coefficients of the visible magnetisation to those of the susceptibility and the inducing field: In the same way, multiplying equation (3.1) byŶ m , +1 (θ, φ) (resp.Ŷ m , (θ, φ)) and integrating over the spherical surface S 1 leads to relations for the hidden part of the magnetisation M i (resp. M t ): where the E m,m ,m , , are the Elsasser integrals: If the susceptibility coefficients χ m can be derived from equation (3.2) then, the g m being known, it is straightforward to derive the hidden part of the magnetisation with equations (3.4). However, the magnetic field model of the lithosphere has a given maximum spherical harmonic degree L that is also, via equation (2.13), the maximum degree of the model of visible magnetisation M e (θ, φ). Therefore there are a maximum of L(L + 2) independent equations (3.2) that can be used to resolve the spherical harmonic coefficients of susceptibility χ m . If the maximum spherical harmonic degree of the (known) inducing field of internal origin B i is L i , then the exclusion rules of the Gaunt integrals (Moon 1979) are such that susceptibility spherical harmonic coefficients of degrees up to L χ = L + L i have to be estimated in order to have all equations (3.2) exactly verified. There is no reason to exclude the spherical harmonic Y 0 0 (θ, φ) for the susceptibility, and therefore there are (L χ + 1) 2 unknown coefficients χ m . This number is necessarily greater than the number of equations (3.2), and there is a significant null-space to this inverse problem where by null-space we mean a combination of susceptibility coefficients that cannot be resolved from the information provided by the magnetic lithospheric field model. The dimension of this null-space drops to (2L + 4) -i.e. (2L χ + 2), when it is assumed that the inducing field is dipolar -i.e. L i = 1, and in that case equations (3.2) reduce to: We made in equation (3.6) a judicious choice of system of coordinates such that the spherical harmonic expansion of the inducing field is The above estimation of dimensions is sufficient to show that a hypothesis of induced, radially average, magnetisation is not a requirement strong enough to define in a unique way the whole magnetisation distribution from a global model of the magnetic field generated in the lithosphere, even when the inducing field is a simple dipole. The structure of the null-space for susceptibility is actually well known. It includes the constant susceptibility layer case (Runcorn 1975) -i.e. χ 0 0 , the χ ±Lχ Lχ coefficients that do not contribute to any of the equations (3.6), and the (2L χ − 1) annihilators described in Maus & Haak (2003). More details are given in the supplementary information provided in Vervelidou & Lesur (2018). In this same work the missing information necessary to resolve uniquely the susceptibility distribution (and therefore the whole magnetisation) is obtained through the synthetic model of radially integrated susceptibility defined by Hemant & Maus (2005).
It is worth mentioning that the above equation (3.6) has similarities with the equations given (2011) has already been identified and set aside to study the problem of characterising the susceptibility.

RECOVERING LOCALISED INDUCED MAGNETISATION
The hypothesis of the previous section is worth investigating further because it is in numerous cases a very good first order approximation of the magnetisation direction. An induced, radially averaged, magnetisation has very often been assumed to interpret local magnetic models of the lithospheric field.
Here, we investigate whether there are favorable settings, when working on a localised area, that reduce the null space. For this we use Slepian functions that have been introduced for geo-potential field studies by Wieczorek & Simons (2005); , and for the particular case of vector field modelling by Plattner & Simons (2014. This technique defines a system of representation of a quantity on the sphere strictly equivalent to a spherical harmonic representation up to a maximum degree L, but where only a finite number of parameters, smaller than (L + 1) 2 , is required to describe accurately this quantity on a small local area. The technique is particularly efficient when working with spherical cap .
The way we proceed is to define the susceptibility distribution on a spherical cap through a Retrieving lithospheric magnetisation distribution 9 limited set of Slepians, and investigate how much they overlap with the null-space for susceptibility defined in section 3, as a function of the aperture and the position of the spherical cap.
An orthogonality between the selected Slepians and the null-space (i.e.-no overlap) implies that the susceptibility distribution can be defined in a unique way inside the cap from the model of visible magnetisation (-and therefore from the magnetic field model). This way to proceed is possible because both Slepians and null-space are just combinations of spherical harmonics and their orthogonality can be established by calculating their scalar products.
The basis functions for the null space are estimated by building the linear system defined by equations (3.6). This defines the matrix H linking the susceptibility Gauss coefficients to the visible magnetisation Gauss coefficients. The null space is then spanned by the eigenvectors of the product H t H associated with null (or vanishingly small) eigenvalues.
Results are presented in figures (1, 2). Each point on these figures is a scalar product between a Slepian and a basis function of the null-space, both being normalised to unity. The scalar products are estimated for all possible combination of selected Slepians and basis functions.
The set of Slepian functions used is such that the susceptibility is well represented on the cap, an remains negligibly small outside the cap. The number of Slepians used is smaller that the number of Gauss coefficients describing the visible magnetisation, and figures (1, 2) show that most of them are not orthogonal to the null space.
In figure (1) is investigated the value of the scalar products as a function of the cap aperture angle. Red crosses are set for the scalar products involving the part of the susceptibility nullspace associated with a constant susceptibility distribution, whereas black crosses are set for the remaining of the null-space. We observe that as the aperture angle decreases, the scalar products get smaller. However they never all become vanishingly small. We conclude that working on a small area never resolves the non-uniqueness of the susceptibility although, the scalar products being small, the contribution of the unknown susceptibility becomes smaller.
We emphasise that these results are strictly valid only in the mathematical framework we have set. In particular, working with Slepians does not allow to have truly local distributions of susceptibility but only localised distributions that never become exactly zero outside the region of interest. In case of local distributions the results from Gerhards (2016) hold.
In figure (2) is investigated the dependence relative to the position of the spherical cap in co-latitude. We expect to observe a dependence because the null-space of the susceptibility  Overall, the part of susceptibility that cannot be resolved from the data becomes smaller as the area of investigation reduces, particularly for caps situated away from the magnetic equator.
Therefore, it is likely that for very small areas the unresolved part of the susceptibility can be simply neglected, and the whole magnetisation reconstructed. We show however in the next section that for very small areas better approximations are possible leading to full knowledge of the whole magnetisation.

MAGNETISATION DIRECTION
We consider in this section the case of constant magnetisation direction. This hypothesis leads to a new set of equations linking the observed magnetic field to the visible and hidden magnetisation that allow, as shown below, to recover the whole magnetisation slightly smoothed in space. We work here again at global scale, with spherical harmonics, as in sections (2,3).
We point out however that a constant magnetisation direction is a hypothesis that is generally not valid in real cases, unless the area of investigation is localised.
An inducing magnetic field constant in space can be represented as a magnetic field B e made of a combination of Y m , −1 vector spherical harmonics of maximum degree L i = 1 . The magnetisation induced by such a field is: (5.1) The q m 1 are the Gauss coefficients for the inducing field. They are defined on the sphere of reference radius a. Multiplying both sides byŶ m , −1 (θ, φ), and integrating over the sphere S 1 leads to a relation linking the spherical harmonic coefficients of the visible magnetisation to those of the susceptibility and the inducing field: where again we made a choice of coordinate system such that the spherical harmonic expansion of the inducing field is B e = q 0 1Ŷ 0 1,0 (θ, φ). We note that in such a model, the inducing field direction is exclusively oriented along the Z-axis of the coordinate system.
The equations (5.3) shows that: (i) If the maximum degree of the spherical harmonic model of the lithospheric magnetic field and associated visible magnetisation is L, then the maximum spherical harmonic degree of susceptibility contributing to the visible magnetisation is L−1 for an inducing field constant in space. We conclude that the number of unknown susceptibility coefficients to be estimated is L 2 .
(ii) The spherical harmonic coefficients of the susceptibility χ ± −1 being not defined, all spherical harmonic coefficients of the visible magnetisation having their orders ±m equal to necessarily vanish -i.e. e ± = 0 for all = 1, 2 · · · L. It follows that the number of independent equations (5.3), that in principle should be L(L + 2) as in section 3, reduces here to L(L + 2) − 2L = L 2 .
These results imply that the L 2 one-to-one relations (5.3) between the visible magnetisation e m and the susceptibility coefficients allows to estimate in an unique way the products q 0 1 χ m for all l = 0, 1, · · · L − 1; there is therefore no null-space for the susceptibility. On the other hand, the direction of magnetisation can be found by searching the orientation of the system of coordinates for which all e ± , = 1, 2 · · · L vanish. This information can be reduced to a single scalar cost function K L 1 (e m ) that has to be zero when the system of coordinates has its Z-axis perfectly aligned with the magnetisation direction. As an example one may choose: where the denominator of the fraction is independent of the orientation of the coordinate system. Since the equation (2.13) sets a one-to-one relationship between magnetic field and magnetisation coefficients, the direction of magnetisation can also be found directly from the

Retrieving lithospheric magnetisation distribution 13
Gauss coefficients of the lithospheric magnetic field using the cost function K a i (g m ).
The products q 0 1 χ m −1 being known from the relations (5.3), we look now for relations linking these products to the hidden magnetisation. We start again from equation (5.1) but now multiply both sides byŶ m , +1 (θ, φ) (resp.Ŷ m , (θ, φ)) and integrate over the spherical surface S 1 . This leads to the relations (5.5) for M i (resp. M t ): The susceptibility model derived from the visible magnetisation goes up to degree L − 1; these new relations (5.5) allow to reconstruct the hidden part of the magnetisation at least up to spherical harmonic degree L − 2 for the i m and L − 1 for the t m . We conclude that if a lithospheric magnetic field model is defined up to spherical harmonic degree L, and that a hypothesis of constant magnetisation direction is valid, then the whole magnetisation can be recovered up to degree L − 2.

APPLICATION TO SYNTHETIC DATA
The process that allows recovering the whole, radially averaged, magnetisation of a body from the magnetic field it generates is illustrated here by applying it to a synthetic data set. The required assumption is that the magnetisation direction is constant over the volume of interest, that for this specific example is the whole Earth crust. This case is described in section 5.
The starting model is made of three dipoles defined in Table (1), all with their magnetisation vector directions defined by the same two angles: the angle θ = 56.145 • between the magnetisation direction and the Z axis of the coordinate system, and the angle φ = 116.565 • between the magnetisation projection on the XY plane and the X axis. This angle is counted positive towards the Y axis. We use here the usual coordinate system where the X, Y, Z axis are the Earth rotation axis for Z, the axis perpendicular to Z pointing towards the Greenwich meridian for X, and the Y direction completes the right-handed orthogonal coordinate system. The location of the dipoles, in particular their depths, and their magnetisation strengths do not correspond to realistic sources but they have been chosen such that they generate a lithospheric magnetic field whose strength is typical for Earth. Its vertical down component is mapped over the Earth on figure (3a). The Gauss coefficientsg m parameterising this field are then estimated up to spherical harmonic degree L = 200. This defines the synthetic data set. When mapped at the Earth surface, this model differs slightly from the figure (3a) because the maximum spherical harmonic degree is too low to represent the sharp features of the field generated by the dipoles. Finally it should be noted that the dipoles generating the field are all magnetised in the same direction, and therefore there is no need to restrict to a localised area the magnetic field; we work here at global scale, using spherical harmonics as our basis functions.
The magnetisation is recovered by first estimating its direction. To test if a direction defined by a pair of angles (θ, φ) is the direction of magnetisation, it is necessary to evaluate the magnetic field Gauss coefficients in the coordinate system that has its Z axis aligned with the (θ, φ) direction. This direction is the magnetisation direction if the cost function K 200 1 (g m ) vanishes. Figure (  Here, we apply our methodology to infer the magnetisation direction, and thus a palepole, using a recent lithospheric magnetic field model of Mars (Langlais et al. 2019). We focus on a known isolated maximum of the magnetic field strength. Isolated anomalies have been considered ideal cases for paleopole estimation since they are thought to reflect clearly the direction of the inducing magnetic field. The one we choose is situated at 143 • co-latitude, -3 • longitude East (see figure 5a). This maximum has been first identified in the model of Morschhauser et al. (2014) and studied in Thomas et al. (2018). In this latter paper, the magnetisation direction is estimated using the Parker (1991) algorithm and found to have a declination D = 188 • and inclination I = 52 • . Using the usual cartesian planetocentric X, Y, Z coordinate system with the Z axis aligned with the planet rotation axis, these declination and inclination values correspond approximately to the angles θ = 75 • , φ = −178 • (or θ = 105 • , φ = 2 • for the antiparallel direction) where θ is the angle between magnetisation direction and the Z axis, while φ is the angle made between the X axis and the projection on the XY plane of the magnetisation direction.
In order to apply our methodology we proceed as follows. We consider a spherical cap of 4.3 • aperture angle, centred on the maximum of the isolated magnetic field anomaly. We generate a global synthetic data set whose values inside this cap are exactly those of the spherical harmonic model of Langlais et al. (2019). Outside this cap the potential field is linearly decreased to zero over a two degree rim. We then process this data set to estimate global spherical harmonic coefficients up to spherical harmonic degree 200. The increase of the maximum spherical harmonic degree from 134, as in the original model, to 200 is necessary to reduce the edge effects of the localisation. The set of Gauss coefficients serve as input to our calculations for the estimation of the underlying magnetisation. Similar results may have been obtained using Slepians to localised the magnetic field model, but we found the approach described here computationally much more efficient for our study. In alignment with our methodology, we assume that over the cap a constant magnetisation direction is a valid hypothesis.
As in the synthetic case of section 6, we search for the direction of magnetisation by minimising the cost function K 200 1 (g m ) defined in equation (5.4); it is mapped in figure (5b). The direction of magnetisation that minimises K 200 1 (g m ) is found for θ = 133.25 • , φ = 6.5 • . Given this direction, the magnetisation distribution can be computed by estimating successively the e m , the products q 0 1 χ m , the i m , t m and the whole magnetisation. These quantities are presented in figure (5) and (6). One can notice in figure (6) the weakness of the toroidal part of the magnetisation M t compared to the synthetic case. It is also interesting to see that the maximum of M e and M i contributions in the X and Z directions do not coincide. Finally, the M e and M i contributions along the Y direction cancel each other to give ultimately a weak whole magnetisation in this direction (see figure 5).
The direction of magnetisation we found corresponds to a local declination and inclination of 36.7 • and −78.4 • respectively -or 216.7 • , 78.4 • for the antiparallel direction. The former numbers are those corresponding to the patterns displayed in figure (5c), choosing the oppo-site sens would leads to a product q 0 1 χ with opposite sign, to give ultimately the same whole magnetisation patterns. As opposed to the Parker method (Parker 1991), there is no need here to impose a positive magnetisation strength.
An important question is the precision of the direction obtained. In this Martian example case, no direction of magnetisation leads to a vanishing small value of the cost function K 200 1 (g m ). This is an indication that the magnetic field model is not compatible with a unique direction of magnetisation, either because the magnetisation direction is truly not constant, or because the magnetic field model is noisy. The functions k (g m ), defined in equation (5.4), and the power spectrum of the magnetic field are plotted as functions of for different orientations of magnetisation (i.e. different coordinate systems) in figure (7). It can be seen that the k (g m ) values are large for all when K 200 1 (g m ) is large; an example is given here for θ = 124 • and φ = 149 • (green crosses). On the other hand, there is a large set of directions for which K 200 1 (g m ) is small, including the direction found by Thomas et al. (2018) (blue crosses), or the planetocentric X, Y, Z coordinate system (black crosses). The values for the optimum direction (θ = 133.25 • , φ = 6.5 • ) are indicated with red crosses. These directions all have k (g m ) values that fall very rapidly with to negligible values. Therefore it is the long wavelength content of the model that controls the direction that minimises the cost function. The power of the model (shown by the solid black line in figure (7)) is weak for these long wavelengths because the magnetic field is only described over a small cap. For the same reasons, the Gauss coefficients for small spherical harmonic degrees are likely to have large variances. It follows that the optimum magnetisation direction cannot be estimated with precision in this case. All directions of magnetisation associated with small values of the cost function (equation 5.4) are such that the associated magnetisation distributions generate very similar magnetic field maps.
The proper way to identify the acceptable directions of magnetisation is to estimate the variance v K of the K 200 1 (g m ) function and to check if a null value is inside one or two standard deviation interval of the K 200 1 (g m ) calculated value. The way the K a i (g m ) variance is estimated is described in appendix-A. Here we assumed that the Gauss coefficients of the original Martian magnetic field model (Langlais et al. 2019) all have the same standard deviation sg m = 0.05nT.
In figure (8) the quantity defined in equation (7.1) is mapped as a function of the magnetisation direction parametrised by (θ, φ). from top to bottom, the X,Y and Z directions. The generated magnetic field is mapped in figure (5) Positive values indicate magnetisation directions where a null value is inside the two standard deviation interval of the K 200 1 (g m ) calculated value. We note that the direction of magnetisation obtained by Thomas et al. (2018) is inside this area of positive values. Therefore, the magnetisation distribution obtained assuming this direction of magnetisation is difficult to distinguish from the one obtained using the optimal direction. Overall, we note the complexity of the patterns in figure (8). It results not only from the localisation and rotation processes required to estimate the K 200 1 (g m ) values, but also from the magnetic field model itself. We expect that the range of possible magnetisation directions reduces for more complex magnetic field distribution and for larger spherical caps. SH degree R(l) X 1.E-3 theta=000, phi=000 theta=133, phi=006 theta=105, phi=002 theta=124, phi=149 Figure 7. Estimated values of the k (g m ) (see equation 5.4, no units) depending on the orientation of the coordinate system. This orientation is provided through two angles, θ and φ, giving the angle made by the new Z axis with the original Z axis, and the angle made by its projection on the XY plane with the X axis, respectively. If the k (g m ) values all vanish for a given coordinate system, its Z axis is then the magnetisation direction generating the observed field.
The power spectrum of the lithospheric field model R( ) (nT 2 ), multiplied by 1.E-3, is also provided.
As a final remark, we note that the approach presented here consists of finding the magnetisation direction by fitting a model of the magnetic field in the spherical harmonic system of representation. Whatever is the direction of magnetisation chosen it is always possible to set magnetisation coefficients that match exactly, through equation (2.13) and (5.3), theg m for all , m as long as m = ± . In contrast, theg ± have to be zero because of equation (5.3).
Therefore the magnetisation direction that provides the best fit to the magnetic field model defines a coordinate system that minimises theg ± . In space domain, this translates in a very specific way of fitting the magnetic field values. This is nonetheless what should in principle be achieved through Parker's method where the magnetisation direction is found using a classic least-squares fit to the observed magnetic field.

CONCLUSION
Being able to derive the rocks magnetisation from the magnetic field it generates is a challenge when dealing with magnetic survey data. To tackle this challenge we proceed on from the results obtained by Gubbins et al. (2011) and investigate how much can be recovered of the whole radially averaged magnetisation of the lithosphere using few simple hypotheses  (5b)). A null K 200 1 (g m ) value is associated with a system of coordinate that has its Z axis aligned with the magnetisation direction. and knowing a magnetic field model of the generated field defined up to spherical harmonic degree L. The results obtained show firstly that a hypothesis of magnetisation induced by a known field of internal origin, even over a localised area, is not sufficient to recover uniquely the radially averaged magnetisation and, secondly, that this magnetisation can be recovered when a constant magnetisation direction is assumed.
The first of these results has already been described in Vervelidou & Lesur (2018), however, we extend it here by studying the case of localised area. We show that as the spherical cap defining the area of investigation becomes smaller, and for caps located close to the magnetic poles, the null space of the magnetisation reduces. For sufficiently small areas, in particular for those close to the magnetic poles, one may want to neglect this part of susceptibility that cannot be recovered, to retain an averaged magnetisation distribution not too different from the true distribution.
The second result is that assuming a constant magnetisation direction allows to recover both this direction and the distribution of radially averaged magnetisation slightly smoothed in space (assuming a magnetic field model is built up to the spherical harmonic degree L, the magnetisation is recovered only up to degree L−2). There is, under the hypothesis of constant magnetisation direction, no restriction on the size of the area of investigation; it can be the whole planet. However, assuming a constant magnetisation direction is generally valid only over a small area.
To find the direction of magnetisation, we have set a cost function K a i (g m ) (see equation (5.4)). This function must be zero when the system of coordinates in which the Gauss coefficients are defined has its Z axis aligned with the direction of magnetisation. When the constant direction of magnetisation is not a valid hypothesis for the studied magnetic field model, then the cost function does not vanish, but can get small in the broad direction of magnetisation of the area investigated. Of course when the variance of the Gauss coefficientsg m is known, it is possible to estimate the variance of the cost function, and in this way determine the directions of magnetisation for which a zero value is inside one, two or three standard-deviations from the actual cost function value. It is therefore possible to provide information on the range of magnetisation directions that may be compatible with the magnetic field model. This possibility to investigate the whole range of acceptable magnetisation directions is a clear advantage of the technique presented here over the Parker (1991) algorithm, and the way it has been applied by Thomas et al. (2018). We expect, therefore, that our methodology can enable progress as far as the study of past dynamo planetary magnetic fields are concerned.
During all numerical tests we ran to recover the magnetisation distribution using our methodology, we noticed that the part of the magnetisation that is most difficult to recover with accuracy is the large scale component. This is a natural effect since large scale magnetised structures tend to produce weak magnetic field and are therefore poorly constrained by the input data. This effect is particularly difficult to handle when working on small areas since, often, the large scale component of the local magnetic field model depends on the technique used to build the model -e.g. Slepians, wavelets, etc.... However, working in spherical geometry over small area is possibly not the best way to go forward; plane geometries are preferable.
For Earth studies, given the complex spatial distribution of magnetisation that can reach very small scales, it is preferable to work with local magnetic surveys for which a plane geometry is generally used. The algorithm we have presented here is not directly applicable for such geometry. Nonetheless, the results obtained by Gubbins et al. (2011) on spherical geometry, have been extended to plane geometry in Gubbins et al. (2017). We have therefore good hope to be able to find a robust algorithm to derive magnetisation orientation and distribution, from a magnetic field model defined in plane geometry. for the production of the maps and figures of the manuscript. This publication has the IPGP number XXXXX.

APPENDIX A: ESTIMATING THE VARIANCE OF THE COST FUNCTION
Let us introduce a perturbation δg m of the magnetic field Gauss coefficients defined in a system of coordinates with its Z axis in the (θ, φ) direction. In this same system of coordinates, the K a i (g m ) function is then also subject to a perturbation: The variance v K of K a i (g m ) is then obtained by: where C δg is the covariance matrix of the δg m . This covariance matrix is in the system of coordinates in which the K a i (g m ) is estimated. Typically, if we have a covariance matrix C 0 δg associated with a global spherical harmonic model of a lithospheric magnetic field defined in the usual planetocentric system of coordinates, the C δg matrix is estimated using: where R (θ,φ) is the operator that estimate the Gauss coefficients in a system of coordinate with its Z axis aligned with the (θ, φ) direction, and P (θ ,φ ) is an operator that defines the magnetic field model inside the spherical cap centred on (θ , φ ) from the original global spherical harmonic model. The super-script t indicates the transpose operator.
In section (7), we used the equations (A.3), (A.4) to compute the variance of K a i (g m ), assuming the covariance matrix of the Martian lithospheric field C 0 δg is diagonal.