A Brief Note on the Computation of Silent from Nonsilent Contributions of Spatially Localized Magnetizations on a Sphere

Abstract. Any square-integrable vector field f over a sphere S can be decomposed into three unique contributions: one being the gradient of a function harmonic inside the sphere (denoted by f+), one being the gradient of a function harmonic in the exterior of the sphere (denoted by f−), and one being tangential and divergence-free (denoted by fdf ). In geomagnetic applications this is of relevance because, if we consider f to be identified with a magnetization, only the contribution f+ can generate a non-vanishing magnetic field in the exterior of the sphere. Thus, we call f− and fdf “silent” and f+ “nonsilent”. If f is known to be spatially localized in a subregion of the sphere, then f+ and f− are coupled due to their potential field nature. In this short paper, we derive an approach that makes use of this coupling in order to compute the contribution f− from knowledge of the contribution f+.


Introduction
The non-uniqueness of the reconstruction of a magnetization M from magnetic field data B is well-known (e.g., [1,5]). Recently this has been investigated in more detail in several publications: for the Euclidean setup with applications in SQUID microscopy in [4,15], for spherical geometries with applications in geomagnetism and planetary magnetism in [3,10,11,12,19,20,14]. The studies above can be divided into those that focus on dealing with M in spectral domain (cf. [12,19,20,14]) and those dealing with spatially localized M (cf. [3,4,10,11,15]). We focus on spatially localized M but use computations in spectral domain. To be precise, throughout the course of the paper, we understand spatial localization as strict spatial localization, i.e., there exists a subregion Γ ⊂ S such that M(x) = 0 for x ∈ S \ Γ.
More precisely, let a square-integrable (vectorial) magnetization M be given on the unit sphere S. Its so-called Hardy-Hodge decomposition takes the form where M + can be expressed as M + = ∇Φ + (with Φ + being harmonic inside the sphere, i.e., in the ball B), M − can be expressed as M − = ∇Φ − (with Φ − being harmonic outside the sphere, i.e., in R 3 \ B), and M df is tangential and divergence-free (more details can be found, e.g., in [2,4,10] and Section 2.2). In general, only the part M + can be reconstructed uniquely from knowledge of the corresponding magnetic field B in the exterior of the sphere S. Therefore, M − and M df are called "silent" while M + is called "nonsilent". If we assume that M is spatially localized, then both M + and M − are uniquely determined by knowledge of B in the exterior. Inverting B subject to localization constraints on M might be one way to obtain the contributions M + and M − (as has been done, e.g., in [4,10,11]). However, the short paper at hand more generally aims at obtaining M − from knowledge of M + under the assumption that (an otherwise unknown) M is spatially localized, without considering any connection to a magnetic field B. The connection to a magnetic field merely serves as one possible application of the approach: M + could be obtained by inversion of B via spectral methods typically used in the geomagnetic community without localization constraints (e.g., [12,19,20,14]). Then we use the localization constraints and the approach proposed in Section 3 to infer M − from this particular M + . Namely, we construct linear operators L 1 , L 2 such that the divergence-free contribution M df is still not determined uniquely, unless one makes further assumptions (e.g., M being of induced type (cf. [11,19,14]), meaning that the magnetization takes the form M = χB 0 , where χ can be interpreted as a scalar-valued susceptibility and B 0 as an ambient magnetic field, which is known or of which at least the general structure is known). Throughout the course of the paper, we typically use the notation f = f + + f − + f df for vector fields on the sphere. We only identify f with a magnetization M when we want to highlight the connection to applications in geomagnetism. Some required notations on vector spherical harmonics and vector field decompositions as well as some auxiliary results on the relation between magnetization and magnetic field are described in Section 2. The desired relation (1.2) between f + and f − is explained in Section 3. The main assertions are reflected in Definition 3.3 and Corollary 3.4. Some numerical examples are presented in Section 4. In the conclusion in Section 5 we additionally provide some outlook for possible further studies.

Basic Notations and Existing Results
First, we define a specific type of vector spherical harmonics that suits the Hardy-Hodge decomposition mentioned in the introduction. Afterwards, we recapitulate the corresponding vector field decomposition and its implications for inverse magnetization problems.

Spherical Harmonics
Spherical harmonics are particularly suitable for considerations on the entire sphere. They are not necessarily the best choice when dealing with subsets of the sphere (then there exist other options such as Slepian functions [16,17], spherical cap harmonics [18], or spherical wavelets [7,9,13]). Although, in the paper at hand, we are actually focusing on magnetizations localized in such subregions, we use spherical harmonics in order to comply with the standard of typical geomagnetic field models (however, the examples later on in Section 4 show that even for fairly simple magnetization models quite high spherical harmonic degrees have to be considered in order to obtain reasonable results). The subsequently defined vector spherical harmonics allow any function f ∈ L 2 (S, C 3 ) to be expanded as where equality is meant in L 2 (S, C 3 )-sense. Analogously, any scalar-valued function f ∈ L 2 (S) can be expressed as The involved coefficientsf − l,m ,f + l,m ,f df l,m ,f l,m are defined in Definition 2.2. Definition 2.1. For degree l ∈ N 0 and order m = −l, . . . , l, we define three types of vector spherical harmonics (e.g. [8,12]): For the special case l = 0, we need to define y + 0,0 (x) = y df 0,0 (x) = 0 in order to avoid inconsistencies. Generally, x = x(θ, φ) = (cos(φ) sin(θ), sin(φ) sin(θ), cos(θ)) T denotes a vector on the unit sphere S = {x ∈ R 3 : |x| = 1} with longitude φ ∈ [0, 2π) and co-latitude θ ∈ [0, π]. Occasionally, we write t = cos(θ) ∈ [−1, 1] to denote the polar distance. By ∇ we denote the the Euclidean gradient, while ∇ S denotes the surface gradient and L S =n × ∇ S the surface curl gradient with respect to the sphere S (n represents the unit normal vector pointing into the exterior of S). For the involved scalar spherical harmonics, we use the definition where P m l denote the Schmidt-nomarlized associated Legendre functions. That is, we have S |Y l,m (y)| 2 dω(y) = 4π and
Some Sobolev spaces that are required later on are gathered in the following definition. They essentially impose a certain smoothness on the functions under consideration and are used for the regularization terms in Section 4 on the numerical implementation.
These Sobolev spaces can equivalently be defined via the corresponding spaces of sequences of Fourier coefficients. Then they are denoted by lower case letters: Sobolev spaces of vector functions are defined as follows: In the special case k = 0, we get H 0 = L 2 (S), H 0 = L 2 (S, C 3 ), and h 0 = 2 (the latter denoting the space of square-summable sequences).

Vector Field Decompositions
The crucial ingredient of this paper is the Hardy-Hodge decomposition: Any vector field f ∈ L 2 (S, C 3 ) can be decomposed uniquely in the form where f df is a tangential and divergence-free vector field, We note that throughout this paper, a vector field f on the sphere is said to be "divergence-free" if ∇ S · P tan [f ] = 0 (with P tan denoting the orthogonal projection onto the tangential contribution of f This leads us to the definition of the following function spaces: They allow an orthogonal decomposition L 2 (S, C 3 ) = H + ⊕H − ⊕H df . We want to mention that so far there exists no common notation for the Hardy-Hodge decomposition in the literature. Instead of the notation (2.7), the same decomposition is denoted as f = E + I + T in [12,19,14,20] or as f =f (1) +f (2) +f (3) in [10,11]. More details on this decomposition and its connection to magnetizations can be found in [2,4] for the Euclidean case and in [10] for the spherical case. Also in these papers can be found the following characterizations for the inverse magnetization problem.
Theorem 2.5. Let the magnetic field B = −∇Φ with the magnetic potential be known for all x ∈ R 3 \B. Furthermore, be M = M + +M − +M df the Hardy-Hodge decomposition of the magnetization M ∈ L 2 (S, C 3 ). (b) Let it be known in advance that the magnetization M ∈ L 2 (S, C 3 ) is spatially localized in Γ ⊂ S. Then M + and M − are determined uniquely by the knowledge of B, while M df ∈ L 2 (S, C 3 ) can be chosen arbitrarily as long as it satisfies the localization constraint.
Since the paper at hand aims at avoiding a detour via the magnetic field B in order to determine the contribution M − of the magnetization M, the following corollary of the theorem above is more relevant.
3 Relating the Fourier Coefficients of f + and f − This very brief section contains the main assertion of the paper. The goal is to connect the vectors of Fourier coefficientsf + ,f − , andf df via linear operators L 1 , L 2 (this is achieved by Definition 3.3 and Corollary 3.4). Therefore, let us assume that the function f ∈ L 2 (S, C 3 ) is spatially localized in Γ. We observe that . . , l. The conditions (3.1) can be equivalently expressed in the form where P + denotes the orthogonal projection onto H + and P − the orthogonal projection onto H − , and χ S\Γ denotes the characteristic function with , the other way around is not necessarily true. But, if a function f ∈ L 2 (S, C 3 ) satisfies (3.2), it follows that χ S\Γ f is divergence-free and, obviously, f −χ S\Γ f is spatially localized in Γ. This implies the following theorem.  To be more precise on the implication that (3.2) renders χ S\Γ f to be divergence-free, we remind the reader that in Section 2.2 it has been stated that any square-integrable vector field f can be expressed in the form f = f + + f − + f df . The conditions P + [f ] = 0 and P − [f ] = 0 imply that f + = f − = 0, so that f = f df . The latter is divergence-free by definition. Thus, (3.2) yields that χ S\Γ f is divergence-free. Furthermore, if f is spatially localized in some Γ ⊂ S, the property (3.2) clearly holds true. We can now add any tangential and divergence-free vector field g that is spatially localized in S \ Γ to obtain somef = f + g that still satisfies P + [χ S\Γf ] = 0 and P − [χ S\Γf ] = 0. However,f is obviously not spatially localized in Γ anymore, which is why (3.  It is to note that we always need to assumef 0,0 =ĥ 0,0 = 0 in order to avoid problems due to the circumstance that we have defined y + 0,0 = y df 0,0 = 0. The next corollary is a direct consequence of the previous setup and Theorem 3.1.

Numerical Implementation and Examples
The assertion of Corollary 3.4 justifies the computation of f − from knowledge of f + by solving The contribution f df cannot be determined uniquely but has to be included as an auxiliary quantity. In order to avoid instabilities due to a possibly unbounded inverse of L 2 , we actually consider the following regularized minimization problem: for some fixed regularization parameter λ > 0. The numerical setup for solving (4.1) is described in the next section.
For the sake of further simplicity, we focus on the case that the subset Γ ⊂ S, in which we assume f to be spatially localized, is a spherical cap with the North Pole e 3 = (0, 0, 1) T as center and with (angular) radius Θ 0 ∈ (0, π), i.e., x · e 3 > cos(Θ 0 )}. This is not a very severe restriction for general applications since typical open subregions of the sphere can be embedded in a spherical cap, and by rotation of the coordinate system this spherical cap can be centered in the North Pole. Such spherical caps have the advantage that most entries of the matrices L 1 and L 2 become zero. The precise computation of the entries is performed in Appendix A.
Remark 4.1. The bandlimited approach for the numerical solution above somehow contradicts the initial assumption that f is spatially localized (since no function can be simultaneously (strictly) spatially localized and bandlimited). However, some sort of discretization has to be done for the numerical evaluation of (4.1) and we chose bandlimitedness for convenience. The precise error due to the bandlimitedness is difficult to quantify analytically. However, the numerical examples in Section 4.2 suggest that the error is minor if the Fourier coefficientsf + l,m decay sufficiently fast while it can be quite significant if such a decay is not given.

Examples
We apply the approach (4.4) to two exemplary functions f , both localized in the northern hemisphere, but the second one being localized in a smaller subregion.

Example 1
First we choose f to be localized in the entire northern hemisphere, i.e., a spherical cap with angular radius Θ 0 = π 2 and the North Pole e 3 = (0, 0, 1) T as center. It is of the form and thus could be considered as a magnetization of induced type, where the ambient inducing field is that of a magnetic dipole with known dipole direction d = (0, 0.6, 0.8) T . The radial part of the function f and its contributions f + and f − are indicated in Figure 1. In fact, one can see that although f is spatially localized in the northern hemisphere, this does not need to hold true for its contributions f + and f − . Reconstructions f N − based on knowledge of the Fourier coefficients of f + up to degree N = 100 and for smoothness parameter k = 2 are shown in Figure 2. For the choice λ = 0.5 · 10 −15 we obtained the best reconstruction of the actual contribution f − among the tested regularization parameters. Sightly overregularized and underregularized results are shown for λ = 10 −12 and λ = 0.5 · 10 −17 . Figures 5 and 3 show the results for the same setup as before, but with noisy input datâ f N,ε we denote a random vector where each component is normally distributed with mean zero and and variance one, scaled to an overall noise level ε N / f N + = 2.5 × 10 −3 . Although the overall noise level is quite small, Figure 5 shows that the relative error for each spherical harmonic degree l and order m can be very significant, in particular for the higher degrees. This is also reflected in the reconstructedf N,ε − in Figure 3: For smoothness parameter k = 2, the general features of f − could be reconstructed, but their amplitude is generally underestimated (for λ = 10 −11 ) or the reconstruction is underregularized and noisy features are introduced (for λ = 10 −12 ). Better results have been obtained for smoothness parameter k = 4, but especially in the northern polar region the reconstruction still reveals some undesired contributions.
Additionally, we varied the truncation degree N for the bandlimitation. The results are indicated in Figure 4. It can be seen that the truncation at a lower degree N = 50 for noise-free    data has a similar effect on the reconstruction (in particular in the northern polar region) as the truncation at a higher degree N = 100 for noisy data. Yet, the effect of the decrease of the bandlimit N is not particularly severe. However, Figure 6 indicates that the chosen function f in this example has a rapidly decreasing power spectrum. Thus, one would not expect a major effect in the first place.

Example 2
We choose the following function that reveals a stronger spatial localization than the first example: The function f and its contributions f + and f − are illustrated in Figure 7. Reconstructions f N − for N = 100 and various choices of smoothing parameters k = 2, 3, 4 are provided in Figure 8 (for each smoothing parameter k, the best result among the tested regularization parameters λ is shown). Overall, it can be said that the results provide a reasonable approximation of the actual contribution f − . However, similar to the noisy setup of the first example, the amplitude of the reconstructions either underestimate the true amplitude of the signal (for k = 2) or some undesirable artifacts in the northern polar regions are introduced (for k = 3, 4). The latter might be explained by the stronger damping of Fourier coefficients at high spherical harmonic degrees as k increases, reducing the ability to approximate more localized features. Overall, the reconstruction of f − seems worse than in Example 1, which is most likely due to the slower decay of the power spectrum towards higher degrees (compare Figure 6). In Figure 9, we have indicated the influence of a varying bandlimit N . Each reconstruction shows the best result for that particular bandlimit among all tested parameters λ and k. It can be seen that the increase in bandlimit reduces artifacts in particular in the northern polar region and that it sharpens the reconstructed features. However, considering the power spectrum indicated in Figure 6, the bandlimit that is required for a precise reconstruction seems significantly higher than the range of spherical harmonic degrees that cover the major contributions of the input signal.

Conclusion
In this paper, we have derived a method on how to numerically connect the contributions f + and f − of a general vector field f under the a priori assumption that f is spatially localized in a subdomain of the sphere. The examples in Section 4 illustrate the capability of the proposed approach but they also indicate the inherent problem of combining localization assumptions and spectral approximation. Other methods have been introduced, e.g., in [14,19], that connect f + and f − via purely spectral arguments and that have been applied to planetary magnetic fields. However, these spectral arguments typically rely on an assumption such as f being of induced type, with knowledge of the general structure of the ambient inducing field (e.g., a constant field). Such an assumption is not necessary if we argue via spatial localization, as discussed in the paper at hand. The assumption of spatial localization substitutes the assumption of f being of induced type. It remains to investigate whether the localization assumption can be reasonably applied to certain setups of planetary magnetic fields. For the numerics, this might require the use of more suitable basis functions, such as Slepians that reflect spatial localization better than spherical harmonics.   For the remaining boundary integrals over ∂Γ, we first observe that ∂ ∂ν = ν · ∇ S and that the unit normal vector ν (pointing into the exterior of S \ Γ) takes the form ν(y) = e3−(e3·y)y sin(Θ0) if Γ is a spherical cap with center e 3 = (0, 0, 1) T and angular radius Θ 0 ∈ (0, π). The unit tangential vector τ takes the form τ (y) = − y×e3 sin(Θ0) . By use of recursion relations for the associated Legendre polynomials (e.g., [8]), we obtain that (2−δm,0) sin(Θ0) P m l (cos(Θ 0 ))P m l (cos(Θ 0 )), else.