Non-singular spherical harmonic expressions of geomagnetic vector and gradient tensor fields in the local north-oriented reference frame

General expressions of magnetic vector (MV) and magnetic gradient tensor (MGT) in terms of the firstand second-order derivatives of spherical harmonics at different degrees/orders are relatively complicated and singular at the poles. In this paper, we derived alternative non-singular expressions for the MV, the MGT and also the third-order partial derivatives of the magnetic potential field in the local north-oriented reference frame. Using our newly derived formulae, the magnetic potential, vector and gradient tensor fields and also the third-order partial derivatives of the magnetic potential field at an altitude of 300 km are calculated based on a global lithospheric magnetic field model GRIMM_L120 (GFZ Reference Internal Magnetic Model, version 0.0) with spherical harmonic degrees 16–90. The corresponding results at the poles are discussed and the validity of the derived formulas is verified using the Laplace equation of the magnetic potential field.


Introduction
Compared to the magnetic vector and scalar measurements, magnetic gradients lead to more robust models of the lithospheric magnetic field.The ongoing Swarm mission of the European Space Agency (ESA) provides measurements not only of the vector and scalar data but also an estimate of their east-west gradients (e.g., Olsen et al., 2004Olsen et al., , 2015;;Friis-Christensen et al., 2006).Kotsiaros andOlsen (2012, 2014) proposed to recover the lithospheric magnetic field through magnetic space gradiometry in the same way that has been done for modeling the gravitational potential field from the satellite gravity gradient tensor measurements by the Gravity field and steady-state Ocean Circulation Explorer (GOCE).Purucker (2005), Purucker et al. (2007), Sabaka et al. (2015) and Kotsiaros et al. (2015) also reported efforts to model the lithospheric magnetic field using magnetic gradient information from the satellite constellation.Their results showed that, by using gradient data, the modeled lithospheric magnetic anomaly field has enhanced shorter wavelength content and a much higher quality compared to models built from vector field data.This is because the gradient data can remove the highly time-dependent contributions of the magnetosphere and ionosphere that are correlated between two side-by-side satellites.
The second-order magnetic gradient tensor consists of spatial derivatives highlighting certain structures of the magnetic field (e.g., Schmidt andClark, 2000, 2006).It can be used to detect the hidden and small-scale magnetized sources (e.g., Pedersen and Rasmussen, 1990;Harrison and Southam, 1991) and to investigate the orientation of the lineated magnetic anomalies (e.g., Blakely and Simpson, 1986).Quantitative magnetic interpretation methods such as the analytic signal, edge detection, spatial derivatives, Euler deconvolution, and transformations, all set in a Cartesian coordinate system (e.g., Blakely, 1995;Purucker and Whaler, 2007;Taylor et al., 2014) also require calculating the higherorder derivatives of the magnetic anomaly field and need to be extended to regional and global scales to handle the curvature of Earth and other planets.Ravat et al. (2002) and Ravat (2011) utilized the analytic signal method and the total Published by Copernicus Publications on behalf of the European Geosciences Union.
gradient to interpret the satellite-altitude magnetic anomaly data.Therefore, both the magnetic field modeling and also the geological interpretations require the calculation for the partial derivatives of the magnetic field, possibly at the poles for specific systems of coordinates.Spherical harmonic analysis, established originally by Gauss (1839), is generally used to model the global magnetic internal fields of Earth and other terrestrial planets (e.g., Maus et al., 2008;Langlais et al., 2010;Thébault et al., 2010;Finlay et al., 2010;Lesur et al., 2013;Sabaka et al., 2013;Olsen et al., 2014).Series of spherical harmonic functions themselves made of Schmidt semi-normalized associated Legendre functions (SSALFs) (e.g., Blakely, 1995;Langel and Hinze, 1998) are fitted by least squares to magnetic measurements, giving the spherical harmonic coefficients (i.e., the Gaussian coefficients) defining the model.Kotsiaros andOlsen (2012, 2014) presented the MV (magnetic vector) and the MGT (magnetic gradient tensor) using a spherical harmonic representation and, of course, their expressions are singular as they approach the poles.Even if there are satellite data gaps around the poles, it is advisable to use non-singular spherical harmonic expressions for the MV and the MGT in case airborne or shipborne magnetic data are utilized (e.g., Golynsky et al., 2013;Maus, 2010).A rotation of the coordinate system is always possible to avoid the polar singularity, but this solution is very ineffective for large data sets.
In this paper, following Petrovskaya and Vershkov (2006) and Eshagh (2008Eshagh ( , 2009) ) for the gravitational gradient tensor in the local north-oriented, orbital reference and geocentric spherical frames, the non-singular expressions in terms of spherical harmonics for the MV, the MGT and the thirdorder derivatives of the magnetic potential field in the specially defined local-north-oriented reference frame (LNORF) are presented.In the next section, the traditional expressions of the MV and the MGT are first stated, some necessary propositions are then proved and, lastly, new non-singular expressions are derived.In Sect.3, the new formulae are tested using the global lithospheric magnetic field model GRIMM_L120 (GFZ Reference Internal Magnetic Model, version 0.0) (Lesur et al., 2013) and compared with the results by traditional formulae.Finally, some conclusions are drawn and further applications are also discussed.

Methodology
In this section, the traditional expressions of MV and MGT are presented and their numerical problems are stated.Then, based on some necessary mathematical derivations, new expressions are given.

Traditional expressions
The scalar potential V of Earth's magnetic field in a sourcefree region can be expanded in the truncated series of spher-ical harmonics at the point P (r, θ , ϕ) with the geocentric distance r, co-latitude θ and longitude ϕ (e.g., Backus et al., 1996): where a = 6371.2km is the radius of Earth's magnetic reference sphere; P m l (cos θ ) (or P m l for simplification) is the SSALF of degree l and order m; L is the maximum spherical harmonic degree; and g m l and h m l are the geomagnetic harmonic coefficients describing Earth's internal sources.
If considered in the LNORF {x, y, z} (e.g., Olsen et al., 2010), where the z axis points downward in the geocentric radial direction, the x axis points to the north, and the y axis towards the east (that is, a right-handed system).At the poles, we define that the x axis points to the meridian of 180 • E (or 180 • W) at the North Pole and of 0 • at the South Pole, which will be discussed in Sect.3. Therefore, the three components of the MV can be expressed as The MGT can be written as (e.g., Kotsiaros and Olsen, 2012) The expressions for V , B z and B zz can be calculated stably even for very high spherical harmonic degrees and orders by using the Holmes and Featherstone (2002a) scheme.However, there exist the singular terms of 1/sinθ and 1/sin 2 θ in Eqs.(2b), (4b), (4d) and (4e) when the computing point approaches to the poles.Moreover, some expressions contain the terms of first-and second-order derivatives of SSALFs, such as Eqs.( 2a) and ( 4a)-(4d).Nevertheless, the up to second-order derivatives for very high degrees and orders of SSALFs can be recursively calculated by the Horner algorithm (Holmes and Featherstone, 2002b).These algorithms are relatively complicated and thus we want to use alternative expressions to avoid the singular terms and also the partial derivatives of SSALFs.It should be stated that our work differs from those presented by Petrovskaya and Vershkov (2006) and Eshagh (2009) in the LNORF and also the associated Legendre functions (ALFs).Nonetheless, the following mathematical derivations are carried out based on their studies on gravity fields.

Mathematical derivations
To deal with the singular terms and first-and second-order derivatives of the SSALFs, some useful mathematical derivations are introduced and proved in the following.

Derivation of d P m
l /dθ Based on Eq. (Z.1.44)in Ilk (1983), and the relation between the ALFs and the SSALFs is thus, the first-order derivative of the SSALFs can be deduced as a l,m = 0.5 b l,m = −0.5 where and δ is the Kronecker delta.

New expressions
Inserting the corresponding mathematical derivations in the last section into Eqs.( 2) and ( 4), and after some simplifications, the new expressions for MV and MGT can be written as g m l cos mλ + h m l sin mφ a zz l,m P m l , where the corresponding coefficients of the SSALFs are given as follows: Furthermore, some other higher-order partial derivatives and their transformations are usually used to image geologic boundaries in magnetic prospecting, such as the higher-order enhanced analytic signal (e.g., Hsu et al., 1996) we also give the third-order partial derivatives of the magnetic potential field as where the corresponding coefficients of the SSALFs are presented as In this way, we avoid computing recursively the SSALFs with singular terms, their first-and second-order derivatives as in the traditional formulae.The cost is only to calculate two additional degrees and orders for the SSALFs at most.It should be noted that, in this study, we use the conventional form of SSALF that if m < 0, then P m l = (−1) |m| P |m| l and if m > l, then P m l = 0.

Numerical investigation and discussion
We test the derived expressions and the numerical implementation in C/C++ by calculating the magnetic potential, magnetic vector and its gradients and also the third-order partial derivatives of the magnetic potential field on a grid with 0.125 • × 0.125 • cell size at the altitude of 300 km relative to Earth's magnetic reference sphere using the lithospheric magnetic field model GRIMM_L120 (version 0.0) defined by Lesur et al. (2013).The magnetic potential, MV, MGT and the third-order partial derivatives of the magnetic potential field in the two polar regions mapped by the lithospheric field model with spherical harmonic degrees 16-90 are shown in Fig. 1 and Fig. 2, respectively.The corresponding statistics around the North Pole and South Pole are, respectively, presented in Tables 1 and 2 (1 Tesla = 10 3 mT = 10 9 nT = 10 12 pT = 10 18 aT), respectively.The relative error is almost equal to the machine's accuracy.Therefore, this feature proves the validity of our derived formulae.In addition, as shown in Figs. 1 and 2, it is obvious that the MGT and also the third-order partial derivatives of the magnetic potential field enhance the lineation and contacts at the satellite altitude.It also reveals some small-scale anomalies, which is very helpful for further geological interpretation.A core field model with spherical harmonic degrees/orders 1-15 is also used for testing and the results, not shown here, indicate the correctness of the formulae in the full range of the spherical harmonic degrees/orders, where the computational stability of the Legendre function with ultrahigh-order is not considered.Furthermore, the computed magnetic fields are smooth near the poles and do not have the singularities, but some components have the dependence on the direction of the reference frame at the poles.As shown in Fig. 3, the magnetic potential V , the B z , B zz and B zzz components at the poles are independent of the direction of the x P and y P axes; while changing with the direction of the x P and y P axes at the poles, the B x , B y , B xz , B yz , B xzz and B yzz components have a period of 360 • and the B xx , B xy , B yy , B xxz , B xyz and B yyz components have a period of 180 • .These variations can be accurately described by a sine or cosine function relating to the horizontal rotation of the reference frame and the differences among these magnetic effects are magnitude, period and initial phase.Therefore, the B x , B y , B xz , B yz , B xx , B xy , B yy , B xzz , B yzz B xxz , B xyz and B yyz components are not smooth at/across the poles.Moreover, to determine the single value at the poles (Figs. 1, 2) we specially define that the x axis points to the meridian of 180 • E (or 180 • W) at the North Pole and of 0 • at the South Pole, that is, the LNORF moving from Greenwich meridian to the poles.
Compared with the traditional formulae in Sect.2.1, there are two advantages of our derived formulae in Sect.2.3.On the one hand, the traditional up to second-order derivatives are removed in the new formulae; therefore, the relatively complicated method by Horner's recursive algorithm (Holmes and Featherstone, 2002b) can be avoided.On the other hand, the singular terms of 1/sinθ and 1/sin 2 θ are removed in the new formulae; consequently, the scale factor of e.g., 10 −280 (Holmes and Featherstone, 2002a, b) is not required when the computing point approaches the poles, and the magnetic fields at the poles can also be calculated in the defined reference frame.In fact, there are differences between the results by our expressions and those by Horner's recursive algorithm; for instance, if using the same model and the parameters as those in Figs. 1 and 2, the differences

Conclusions
We develop in this paper the new expressions for the MV, the MGT and the third-order partial derivatives of the magnetic potential field in terms of spherical harmonics.The traditional expressions have complicated forms involving firstand second-order derivatives of the SSALFs and are singular when approaching the poles.Our newly derived formulae do not contain the first-and second-order derivatives of the SSALFs and remove the singularities at the poles.However, our formulae are derived in the spherical LNORF with specific definition at the poles.For a future application to the magnetic data of a satellite gradiometry mission (e.g., Kotsiaros and Olsen, 2014), it is necessary to describe the MV and the MGT in the local orbital or other reference frame, where the new MV and MGT are the linear functions of the MV and the MGT in the LNORF with coefficients related to the satellite track azimuth (e.g., Petrovskaya and Vershkov, 2006) or other rotation angles.The other main purpose of this paper is, in the future, to contribute to the signal processing and the geophysical and geological interpretation of the global lithospheric magnetic field model, especially near polar areas. of Geosciences, Wuhan) (grant no.SMIL-2015-06) and State Key Laboratory of Geodesy and Earth's Dynamics (Institute of Geodesy and Geophysics, CAS) (grant no.SKLGED2015-5-5-EZ).Jinsong Du is sponsored by the China Scholarship Council (CSC).We would like to thank Mehdi Eshagh and another anonymous reviewer for their constructive comments.All projected figures are drawn using the Generic Mapping Tools (GMT) (Wessel and Smith, 1991).
The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association.

Figure 1 .
Figure 1.Lithospheric magnetic potential, magnetic vector and its gradient fields and third-order partial derivatives of the magnetic potential field around the North Pole (0 • ≤ θ ≤ 30 • ) at the altitude of 300 km as defined by the lithospheric magnetic field model GRIMM_L120 (version 0.0) (Lesur et al., 2013) for spherical harmonic degrees 16-90.(a) is magnetic potential (V ); (b), (c) and (d) are three components (B x , B y and B z ) of the magnetic vector; (e), (f), (g), (h), (i) and (j) are six elements (B xx , B xy , B xz , B yy ,B yz and B zz ) of the magnetic gradient tensor; (k), (l), (m), (n), (o) and (p) are six elements (B xxz , B xyz , B xzz , B yyz ,B yzz and B zzz ) of third-order partial derivatives of the magnetic potential field, respectively.The dark green lines are the plate boundaries by Bird (2003).All maps are shown in polar stereographic projections.

Figure 2 .
Figure 2. Lithospheric magnetic potential, magnetic vector and its gradient fields and third-order partial derivatives of the magnetic potential field around the South Pole (150 • ≤ θ ≤ 180 • ) at the altitude of 300 km as defined by the lithospheric magnetic field model GRIMM_L120 (version 0.0) (Lesur et al., 2013) for spherical harmonic degrees 16-90.(a) is magnetic potential (V ), (b), (c) and (d) are three components (B x , B y and B z ) of the magnetic vector; (e), (f), (g), (h), (i) and (j) are six elements (B xx , B xy , B xz , B yy ,B yz andB zz ) of the magnetic gradient tensor; (k), (l), (m), (n), (o) and (p) are six elements (B xxz , B xyz , B xzz , B yyz ,B yzz and B zzz ) of third-order partial derivatives of the magnetic potential field, respectively.The dark green lines are the plate boundaries by Bird (2003).All maps are shown in polar stereographic projections.

Figure 3 .
Figure 3. Limit values of the magnetic potential (V ), magnetic vector (B x , B y and B z ) and its gradients (B xx , B xy , B xz , B yy ,B yz andB zz ) and third-order partial derivatives of the magnetic potential field (B xxz , B xyz , B xzz , B yyz ,B yzz and B zzz ) at the poles when the local reference frames vary from different meridians (the direction of x P axes changing from different meridians to the poles).Red and blue lines indicate the magnetic effects at the North Pole and at the South Pole, respectively.The reference frame is specially defined that the x P axis points to the meridian of 180 • E (or 180 • W) at the North Pole and 0f 0 • at the South Pole and the y P axis points to the meridian of 90 • E at both poles.The values at both poles shown by black dashed arrows are used to plot the maps in Figs. 1 and 2.
. A simple test is that the MGT meets Laplace's equation of the potential field; that is, the trace of the MGT should be equal to zero.Our numerical results show that the amplitudes of B xx + B yy + B zz in the North Pole and South Pole regions are in the range of [−2.012 × 10 −15 pT m −1 to +2.026 × 10 −15 pT m −1 ]