A brief note on the computation of silent from nonsilent contributions of spatially localized magnetizations on a sphere

Abstract

Any square-integrable vector field $\mathbf {f}$ over a sphere $\mathbb {S}$ can be decomposed into three unique contributions: one being the gradient of a function harmonic inside the sphere (denoted by $\mathbf {f}_+$), one being the gradient of a function harmonic in the exterior of the sphere (denoted by $\mathbf {f}_-$), and one being tangential and divergence-free (denoted by $\mathbf {f}_{df}$). In geomagnetic applications this is of relevance because, if we consider $\mathbf {f}$ to be identified with a magnetization, only the contribution $\mathbf {f}_+$ can generate a non-vanishing magnetic field in the exterior of the sphere. Thus, we call $\mathbf {f}_-$ and $\mathbf {f}_{df}$ “silent” and $\mathbf {f}_+$ “nonsilent”. If $\mathbf {f}$ is known to be spatially localized in a subregion of the sphere, then $\mathbf {f}_+$ and $\mathbf {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 $\mathbf {f}_-$ from knowledge of the contribution $\mathbf {f}_+$.

1 Introduction

The non-uniqueness of the reconstruction of a magnetization $\mathbf {M}$ from magnetic field data $\mathbf {B}$ is well-known (e.g., Backus et al. 1996; Blakely 1995). Recently this has been investigated in more detail in several publications: for the Euclidean setup with applications in SQUID microscopy in Baratchart et al. (2013) and Lima et al. (2013), for spherical geometries with applications in geomagnetism and planetary magnetism in Baratchart and Gerhards (2017), Gerhards (2016), Gerhards (2019), Gubbins et al. (2011), Vervelidou and Lesur (2018), Vervelidou et al. (2017) and Lesur and Vervelidou (2020). The studies above can be divided into those that focus on dealing with $\mathbf {M}$ in spectral domain (cf. Gubbins et al. 2011; Vervelidou and Lesur 2018; Vervelidou et al. 2017; Lesur and Vervelidou 2020) and those dealing with spatially localized $\mathbf {M}$ (cf. Baratchart and Gerhards 2017; Baratchart et al. 2013; Gerhards 2016, 2019; Lima et al. 2013). We focus on spatially localized $\mathbf {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 $\Gamma \subset \mathbb {S}$ such that $\mathbf {M}(x)=0$ for $x\in \mathbb {S}{\setminus }{\overline{\Gamma }}$.

More precisely, let a square-integrable (vectorial) magnetization $\mathbf {M}$ be given on the unit sphere $\mathbb {S}$. Its so-called Hardy–Hodge decomposition takes the form

$$\begin{aligned} \mathbf {M}=\mathbf {M}_++\mathbf {M}_-+\mathbf {M}_{df}, \end{aligned}$$

(1.1)

where $\mathbf {M}_+$ can be expressed as $\mathbf {M}_+=\nabla \Phi _+$ (with $\Phi _+$ being harmonic inside the sphere, i.e., in the ball $\mathbb {B}$), $\mathbf {M}_-$ can be expressed as $\mathbf {M}_-=\nabla \Phi _-$ (with $\Phi _-$ being harmonic outside the sphere, i.e., in $\mathbb {R}^3{\setminus }{\overline{\mathbb {B}}}$), and $\mathbf {M}_{df}$ is tangential and divergence-free (more details can be found, e.g., in Baratchart et al. (2018), Baratchart et al. (2013), Gerhards (2016), and Sect. 2.2). In general, only the part $\mathbf {M}_+$ can be reconstructed uniquely from knowledge of the corresponding magnetic field $\mathbf {B}$ in the exterior of the sphere $\mathbb {S}$. Therefore, $\mathbf {M}_-$ and $\mathbf {M}_{df}$ are called “silent” while $\mathbf {M}_+$ is called “nonsilent”. If we assume that $\mathbf {M}$ is spatially localized, then both $\mathbf {M}_+$ and $\mathbf {M}_-$ are uniquely determined by knowledge of $\mathbf {B}$ in the exterior.

Inverting $\mathbf {B}$ subject to localization constraints on $\mathbf {M}$ might be one way to obtain the contributions $\mathbf {M}_+$ and $\mathbf {M}_-$ (as has been done, e.g., in Baratchart et al. (2013), Gerhards (2016) and Gerhards (2019)). However, the short paper at hand more generally aims at obtaining $\mathbf {M}_-$ from knowledge of $\mathbf {M}_+$ under the assumption that (an otherwise unknown) $\mathbf {M}$ is spatially localized, without considering any connection to a magnetic field $\mathbf {B}$. The connection to a magnetic field merely serves as one possible application of the approach: $\mathbf {M}_+$ could be obtained by inversion of $\mathbf {B}$ via spectral methods typically used in the geomagnetic community without localization constraints (e.g., Gubbins et al. 2011; Vervelidou and Lesur 2018; Vervelidou et al. 2017; Lesur and Vervelidou 2020). Then we use the localization constraints and the approach proposed in Sect. 3 to infer $\mathbf {M}_-$ from this particular $\mathbf {M}_+$. Namely, we construct linear operators $L_1$, $L_2$ such that

$$\begin{aligned} L_1[\hat{\mathbf {M}}_+]=L_2[\hat{\mathbf {M}}_-,\hat{\mathbf {M}}_{df}], \end{aligned}$$

(1.2)

where $\hat{\mathbf {M}}_+$, $\hat{\mathbf {M}}_-$, $\hat{\mathbf {M}}_{df}$ denote the vectors of Fourier coefficients of $\mathbf {M}_+$, $\mathbf {M}_-$, $\mathbf {M}_{df}$, respectively. Solving (1.2) for $\hat{\mathbf {M}}_-$ and $\hat{\mathbf {M}}_{df}$ yields the desired result. However, while $\mathbf {M}_-$ is determined uniquely, the divergence-free contribution $\mathbf {M}_{df}$ is still not determined uniquely, unless one makes further assumptions (e.g., $\mathbf {M}$ being of induced type (cf. Gerhards 2019; Vervelidou and Lesur 2018; Lesur and Vervelidou 2020), meaning that the magnetization takes the form $\mathbf {M}=\chi \mathbf {B}_0$, where $\chi$ can be interpreted as a scalar-valued susceptibility and $\mathbf {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 $\mathbf {f}=\mathbf {f}_++\mathbf {f}_-+\mathbf {f}_{df}$ for vector fields on the sphere. We only identify $\mathbf {f}$ with a magnetization $\mathbf {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 Sect. 2. The desired relation (1.2) between $\mathbf {f}_+$ and $\mathbf {f}_-$ is explained in Sect. 3. The main assertions are reflected in Definition 3.3 and Corollary 3.4. Some numerical examples are presented in Sect. 4. In the conclusion in Sect. 5 we additionally provide some outlook for possible further studies.

2 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.

2.1 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 (Plattner and Simons 2014; Simons et al. 2006), spherical cap harmonics (Thébault et al. 2006), or spherical wavelets Freeden et al. 1998; Freeden and Windheuser 1997; Holschneider et al. 2003). 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 Sect. 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 $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ to be expanded as

$$\begin{aligned} \mathbf {f}=\sum _{l=0}^\infty \sum _{m=-l}^l\hat{f}_{l,m}^- \mathbf {y}_{l,m}^-+ \hat{f}_{l,m}^+ \mathbf {y}_{l,m}^++ \hat{f}_{l,m}^{df} \mathbf {y}_{l,m}^{df}, \end{aligned}$$

(2.1)

where equality is meant in $\text {L}^2(\mathbb {S},\mathbb {C}^3)$-sense. Analogously, any scalar-valued function $f\in \text {L}^2(\mathbb {S})$ can be expressed as

$$\begin{aligned} f=\sum _{l=0}^\infty \sum _{m=-l}^l\hat{f}_{l,m} Y_{l,m}. \end{aligned}$$

(2.2)

The involved coefficients $\hat{f}_{l,m}^-$, $\hat{f}_{l,m}^+$, $\hat{f}_{l,m}^{df}$, $\hat{f}_{l,m}$ are defined in Definition 2.2.

Definition 2.1

For degree $l\in \mathbb {N}_0$ and order $m=-l,\ldots ,l$, we define three types of vector spherical harmonics (e.g. Freeden and Schreiner 2009; Gubbins et al. 2011):

$$\begin{aligned}&\mathbf {y}_{l,m}^-(x)=\mathbf {y}_{l,m}^-(\theta ,\phi )=\frac{1}{\sqrt{(l+1)(2l+1)}}\nabla _y\left( \frac{1}{|y|^{l+1}}Y_{l,m}(y)\right) \bigg |_{y=x}, \end{aligned}$$

(2.3)

$$\begin{aligned}&\mathbf {y}_{l,m}^+(x)=\mathbf {y}_{l,m}^+(\theta ,\phi )=\frac{1}{\sqrt{l(2l+1)}}\nabla _y\left( |y|^lY_{l,m}(y)\right) \Big |_{y=x}, \end{aligned}$$

(2.4)

$$\begin{aligned}&\mathbf {y}_{l,m}^{df}(x)=\mathbf {y}_{l,m}^{df}(\theta ,\phi )=-\frac{i}{\sqrt{l(l+1)}}\text {L}_\mathbb {S}Y_{l,m}(x). \end{aligned}$$

(2.5)

For the special case $l=0$, we need to define $\mathbf {y}_{0,0}^+(x)=\mathbf {y}_{0,0}^{df}(x)=0$ in order to avoid inconsistencies. Generally, $x=x(\theta ,\phi )=(\cos (\phi )\sin (\theta ),\sin (\phi )\sin (\theta ),\cos (\theta ))^T$ denotes a vector on the unit sphere $\mathbb {S}=\{x\in \mathbb {R}^3:|x|=1\}$ with longitude $\phi \in [0,2\pi )$ and co-latitude $\theta \in [0,\pi ]$. Occasionally, we write $t=\cos (\theta )\in [-1,1]$ to denote the polar distance. By $\nabla$ we denote the the Euclidean gradient, while $\nabla _\mathbb {S}$ denotes the surface gradient and $\text {L}_\mathbb {S}=\hat{\mathbf {n}}\times \nabla _\mathbb {S}$ the surface curl gradient with respect to the sphere $\mathbb {S}$ ($\hat{\mathbf {n}}$ represents the unit normal vector pointing into the exterior of $\mathbb {S}$). For the involved scalar spherical harmonics, we use the definition

$$\begin{aligned} Y_{l,m}(x)=\sqrt{\frac{2l+1}{2-\delta _{m,0}}}P_l^m(\cos (\theta )) e^{im\phi }, \end{aligned}$$

(2.6)

where $P_l^m$ denote the Schmidt-nomarlized associated Legendre functions. That is, we have $\int _\mathbb {S}|Y_{l,m}(y)|^2d\omega (y)=4\pi$ and $\int _{-1}^1 P_l^m(t)P_{l'}^{m}(t)dt=\frac{2(2-\delta _{m,0})}{2l+1}\delta _{l,l'}$.

Definition 2.2

For a vectorial square-integrable function $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$, the Fourier coefficients in (2.1) are defined via

$$\begin{aligned} \hat{f}_{l,m}^-&=\langle \mathbf {f},\mathbf {y}_{l,m}^-\rangle _{\text {L}^2(\mathbb {S},\mathbb {C}^3)}= \frac{1}{4\pi }\int _\mathbb {S}\mathbf {f}(y)\cdot \mathbf {y}_{l,m}^-(y)^* d\omega (y), \\ \hat{f}_{l,m}^+&=\langle \mathbf {f},\mathbf {y}_{l,m}^+\rangle _{\text {L}^2(\mathbb {S},\mathbb {C}^3)}=\frac{1}{4\pi }\int _\mathbb {S}\mathbf {f}(y)\cdot \mathbf {y}_{l,m}^+(y)^* d\omega (y), \\ \hat{f}_{l,m}^{df}&=\langle \mathbf {f},\mathbf {y}_{l,m}^{df}\rangle _{\text {L}^2(\mathbb {S},\mathbb {C}^3)}=\frac{1}{4\pi }\int _\mathbb {S}\mathbf {f}(y)\cdot \mathbf {y}_{l,m}^{df}(y)^* d\omega (y). \end{aligned}$$

The collection of all Fourier coefficients is denoted by vectors ${\hat{\mathbf {f}}}_+=(\hat{f}_{0,0}^+,\hat{f}_{1,-1}^{+},\hat{f}_{1,0}^+,\hat{f}_{1,1}^+,\hat{f}_{2,-2}^{+},\ldots )^T$, ${\hat{\mathbf {f}}}_-=(\hat{f}_{0,0}^-,\hat{f}_{1,-1}^{-},\hat{f}_{1,0}^-,\hat{f}_{1,1}^-,\hat{f}_{2,-2}^{-},\ldots )^T$, ${\hat{\mathbf {f}}}_{df}=(\hat{f}_{0,0}^{df},\hat{f}_{1,-1}^{df},\hat{f}_{1,0}^{df},\hat{f}_{1,1}^{df},\hat{f}_{2,-2}^{df},\ldots )^T$. The case case $l=0$ is again special in the sense that $\hat{f}_{0,0}^+=\hat{f}_{0,0}^{df}=0$ in all cases. For scalar-valued functions $f\in \text {L}^2(\mathbb {S})$, we analogously set

$$\begin{aligned} \hat{f}_{l,m}&=\frac{1}{4\pi }\int _\mathbb {S}f(y)Y_{l,m}(y)^* d\omega (y). \end{aligned}$$

The collection of Fourier coefficients is denoted by the vector $\hat{f}=(f_{0,0},f_{1,-1},f_{1,0},f_{1,1},f_{2,-2},\ldots )^T$.

The vector spherical harmonics (2.3)–(2.5) correspond (up to a possibly varying sign or normalization factor) to $\mathbf {Y}_{l,l+1}^m$, $\mathbf {Y}_{l,l-1}^m$, $\mathbf {Y}_{l,l}^m$ in Gubbins et al. (2011), Vervelidou and Lesur (2018) and Lesur and Vervelidou (2020) and $\tilde{\mathbf {y}}_{l,m}^{(1)}$, $\tilde{\mathbf {y}}_{l,m}^{(2)}$, $\tilde{\mathbf {y}}_{l,m}^{(3)}$ in Freeden and Gerhards (2012), Freeden and Schreiner (2009), Gerhards (2016) and Gerhards (2019).

Definition 2.3

Throughout this paper, $\Gamma \subset \mathbb {S}$ denotes an open subregion of the sphere with smooth boundary and ${\overline{\Gamma }}\not =\mathbb {S}$. A function $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ is said to be spatially localized in $\Gamma$ if $\mathbf {f}(x)=0$, for almost all $x\in \mathbb {S}{\setminus }{\overline{\Gamma }}$, i.e., $\overline{\text {supp}(\mathbf {f})}\subset \Gamma$.

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 Sect. 4 on the numerical implementation.

Definition 2.4

For $k\in \mathbb {N}_0$, Sobolev spaces of scalar functions are defined as follows:

$$\begin{aligned} \text {H}_k&=\left\{ f\in \text {L}^2(\mathbb {S}):\Vert f\Vert _{\text {H}_k}^2=\sum _{l=0}^\infty \sum _{m=-l}^l\left( l+\frac{1}{2}\right) ^{2k}\left| \hat{f}_{l,m}\right| ^{2}<\infty \right\} . \end{aligned}$$

These Sobolev spaces can equivalently be defined via the corresponding spaces of sequences of Fourier coefficients. Then they are denoted by lower case letters:

$$\begin{aligned} \text {h}_k&=\left\{ \hat{f}=\left( \hat{f}_{0,0},\hat{f}_{1,-1},\hat{f}_{1,0},\hat{f}_{1,1},\hat{f}_{2,-2},\ldots \right) :\hat{f}_{l,m}\in \mathbb {C},\,\Vert \hat{f}\Vert _{\text {h}_k}^2=\sum _{l=0}^\infty \sum _{m=-l}^l\left( l+\frac{1}{2}\right) ^{2k}\left| \hat{f}_{l,m}\right| ^{2}<\infty \right\} . \end{aligned}$$

Sobolev spaces of vector functions are defined as follows: $\mathbf {H}_k=\{\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3):\hat{\mathbf {f}}_+,\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df}\in \text {h}_k\}$. In the special case $k=0$, we get $\text {H}_0=\text {L}^2(\mathbb {S})$, $\mathbf {H}_0=\text {L}^2(\mathbb {S},\mathbb {C}^3)$, and $\text {h}_0=\ell _2$ (the latter denoting the space of square-summable sequences).

2.2 Vector field decompositions

The crucial ingredient of this paper is the Hardy–Hodge decomposition: Any vector field $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ can be decomposed uniquely in the form

$$\begin{aligned} \mathbf {f}=\mathbf {f}_++\mathbf {f}_-+\mathbf {f}_{df}, \end{aligned}$$

(2.7)

where $\mathbf {f}_{df}$ is a tangential and divergence-free vector field, $\mathbf {f}_+=\nabla \Phi _+|_\mathbb {S}$ is the gradient of a function $\Phi _+\in \text {L}^2({\overline{\mathbb {B}}})$ that is harmonic in the interior of the unit ball $\mathbb {B}=\{x\in \mathbb {R}^3:|x|<1\}$ and $\mathbf {f}_-=\nabla \Phi _-|_\mathbb {S}$ is the gradient of a function $\Phi _-\in \text {L}^2(\mathbb {R}^3{\setminus }{\mathbb {B}})$ that is harmonic in the exterior of the unit ball $\mathbb {R}^3{\setminus }{\overline{\mathbb {B}}}$. We note that throughout this paper, a vector field $\mathbf {f}$ on the sphere is said to be “divergence-free” if $\nabla _\mathbb {S}\cdot P_{tan}[\mathbf {f}]=0$ (with $P_{tan}$ denoting the orthogonal projection onto the tangential contribution of $\mathbf {f}$). With the vector spherical harmonics from Sect. 2.1 and the definition of Fourier coefficients in Definition 2.2, it becomes clear that the contributions $\mathbf {f}_+$, $\mathbf {f}_-$, and $\mathbf {f}_{df}$ of a vector-valued function $\mathbf {f}$ can be represented by

$$\begin{aligned}&\mathbf {f}_-=\sum _{l=0}^\infty \sum _{m=-l}^l\hat{f}_{l,m}^- \mathbf {y}_{l,m}^-, \end{aligned}$$

(2.8)

$$\begin{aligned}&\mathbf {f}_+=\sum _{l=0}^\infty \sum _{m=-l}^l\hat{f}_{l,m}^+ \mathbf {y}_{l,m}^+, \end{aligned}$$

(2.9)

$$\begin{aligned}&\mathbf {f}_{df}=\sum _{l=0}^\infty \sum _{m=-l}^l\hat{f}_{l,m}^{df} \mathbf {y}_{l,m}^{df}. \end{aligned}$$

(2.10)

This leads us to the definition of the following function spaces:

$$\begin{aligned} \mathbf {H}_+&=\overline{\text {span}\{\mathbf {y}_{l,m}^+,l\in \mathbb {N}_0,m=-l,\ldots ,l\}}^{\Vert \cdot \Vert _{\text {L}^2(\mathbb {S},\mathbb {C}^3)}} \end{aligned}$$

(2.11)

$$\begin{aligned} \mathbf {H}_-&=\overline{\text {span}\{\mathbf {y}_{l,m}^-,l\in \mathbb {N}_0,m=-l,\ldots ,l\}}^{\Vert \cdot \Vert _{\text {L}^2(\mathbb {S},\mathbb {C}^3)}} \end{aligned}$$

(2.12)

$$\begin{aligned} \mathbf {H}_{df}&=\overline{\text {span}\{\mathbf {y}_{l,m}^{df},l\in \mathbb {N}_0,m=-l,\ldots ,l\}}^{\Vert \cdot \Vert _{\text {L}^2(\mathbb {S},\mathbb {C}^3)}} \end{aligned}$$

(2.13)

They allow an orthogonal decomposition $\text {L}^2(\mathbb {S},\mathbb {C}^3)=\mathbf {H}_+\oplus \mathbf {H}_-\oplus \mathbf {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 $\mathbf {f}=\mathcal {E}+\mathcal {I}+\mathcal {T}$ in Gubbins et al. (2011), Vervelidou and Lesur (2018), Lesur and Vervelidou (2020) and Vervelidou et al. (2017) or as $\mathbf {f}={\tilde{\mathbf {f}}}^{(1)}+{\tilde{\mathbf {f}}}^{(2)}+{\tilde{\mathbf {f}}}^{(3)}$ in Gerhards (2016) and Gerhards (2019). More details on this decomposition and its connection to magnetizations can be found in Baratchart et al. (2018) and Baratchart et al. (2013) for the Euclidean case and in Gerhards (2016) 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 $\mathbf {B}=-\nabla \Phi$ with the magnetic potential

$$\begin{aligned} \Phi (x)=\frac{\mu _0}{4\pi }\int _\mathbb {S}\mathbf {M}(y)\cdot \frac{x-y}{|x-y|^3}\,d\omega (y) \end{aligned}$$

be known for all $x\in \mathbb {R}^3{\setminus }{\overline{\mathbb {B}}}$. Furthermore, be $\mathbf {M}=\mathbf {M}_++\mathbf {M}_-+\mathbf {M}_{df}$ the Hardy–Hodge decomposition of the magnetization $\mathbf {M}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$.

• (a) The contribution $\mathbf {M}_+$ is determined uniquely by the knowledge of $\mathbf {B}$, while $\mathbf {M}_-,\mathbf {M}_{df}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ can be chosen arbitrarily.

• (b) Let it be known in advance that the magnetization $\mathbf {M}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ is spatially localized in $\Gamma \subset \mathbb {S}$. Then $\mathbf {M}_+$ and $\mathbf {M}_-$ are determined uniquely by the knowledge of $\mathbf {B}$, while $\mathbf {M}_{df}\in \text {L}^2(\mathbb {S},\mathbb {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 $\mathbf {B}$ in order to determine the contribution $\mathbf {M}_-$ of the magnetization $\mathbf {M}$, the following corollary of the theorem above is more relevant.

Corollary 2.6

Let the function $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ be spatially localized in $\Gamma \subset \mathbb {S}$. Then $\mathbf {f}_-$ is uniquely determined by the knowledge of $\mathbf {f}_+$.

3 Relating the fourier coefficients of $\mathbf {f}_+$ and $\mathbf {f}_-$

This very brief section contains the main assertion of the paper. The goal is to connect the vectors of Fourier coefficients $\hat{\mathbf {f}}_+$, $\hat{\mathbf {f}}_-$, and $\hat{\mathbf {f}}_{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 $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ is spatially localized in $\Gamma$. We observe that

$$\begin{aligned} \int _{\mathbb {S}{\setminus }\Gamma } \mathbf {f}(y)\cdot \mathbf {y}_{l,m}^-(y)^* \,d\omega (y)=0, \qquad \int _{\mathbb {S}{\setminus }\Gamma } \mathbf {f}(y)\cdot \mathbf {y}_{l,m}^+(y)^* \,d\omega (y)=0, \end{aligned}$$

(3.1)

for all $l\in \mathbb {N}_0$, $m=-l,\ldots ,l$. The conditions (3.1) can be equivalently expressed in the form

$$\begin{aligned} P_+[\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}]=0,\qquad P_-[\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}]=0, \end{aligned}$$

(3.2)

where $P_+$ denotes the orthogonal projection onto $\mathbf {H}_+$ and $P_-$ the orthogonal projection onto $\mathbf {H}_-$, and $\chi _{\mathbb {S}{\setminus }\Gamma }$ denotes the characteristic function with $\chi _{\mathbb {S}{\setminus }\Gamma }(x)=1$, if $x\in \mathbb {S}{\setminus }\Gamma$, and $\chi _{\mathbb {S}{\setminus }\Gamma }(x)=0$, if $x\in \Gamma$. While it is immediately clear that any function $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ that is spatially localized in $\Gamma \subset \mathbb {S}$ satisfies (3.2), the other way around is not necessarily true. But, if a function $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ satisfies (3.2), it follows that $\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}$ is divergence-free and, obviously, $\mathbf {f}-\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}$ is spatially localized in $\Gamma$. This implies the following theorem.

Theorem 3.1

If $\mathbf {f}\in \text {L}^2(\mathbb {S},\mathbb {C}^3)$ satisfies (3.2), then the knowledge of $\mathbf {f}_+$ uniquely determines $\mathbf {f}_-$.

Proof

Corollary 2.6 implies that the contribution $(\mathbf {f}-\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f})_+$ of $\mathbf {f}-\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}$ uniquely determines the contribution $(\mathbf {f}-\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f})_-$. Since (3.2) yields that $\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}$ is divergence-free, we obtain that $\mathbf {f}_+$ and $(\mathbf {f}-\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f})_+$ coincide as well as $\mathbf {f}_-$ and $(\mathbf {f}-\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f})_-$, which already proves the assertion of the theorem. $\square$

Remark 3.2

To be more precise on the implication that (3.2) renders $\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}$ to be divergence-free, we remind the reader that in Sect. 2.2 it has been stated that any square-integrable vector field $\mathbf {f}$ can be expressed in the form $\mathbf {f}=\mathbf {f}_++\mathbf {f}_-+\mathbf {f}_{df}$. The conditions $P_+[\mathbf {f}]=0$ and $P_-[\mathbf {f}]=0$ imply that $\mathbf {f}_+=\mathbf {f}_-=0$, so that $\mathbf {f}=\mathbf {f}_{df}$. The latter is divergence-free by definition. Thus, (3.2) yields that $\chi _{\mathbb {S}{\setminus }\Gamma }\,\mathbf {f}$ is divergence-free. Furthermore, if $\mathbf {f}$ is spatially localized in some $\Gamma \subset \mathbb {S}$, the property (3.2) clearly holds true. We can now add any tangential and divergence-free vector field $\mathbf {g}$ that is spatially localized in $\mathbb {S}{\setminus }{\overline{\Gamma }}$ to obtain some ${\bar{\mathbf {f}}}=\mathbf {f}+\mathbf {g}$ that still satisfies $P_+[\chi _{\mathbb {S}{\setminus }\Gamma }\,{\bar{\mathbf {f}}}]=0$ and $P_-[\chi _{\mathbb {S}{\setminus }\Gamma }\,{\bar{\mathbf {f}}}]=0$. However, ${\bar{\mathbf {f}}}$ is obviously not spatially localized in $\Gamma$ anymore, which is why (3.2) does not necessarily imply spatial localization of $\mathbf {f}$ in $\Gamma$.

If $\mathbf {f}\in \mathbf {H}_2$, then (3.1) and (3.2) are just alternative notations for the equation $L_1[\hat{\mathbf {f}}_+]=L_2[\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df}]$ with operators $L_1$, $L_2$ defined as below. Latter equation, though less intuitive than (3.1) and (3.2), is more explicit in terms of its numerical implementation in Sect. 4.

Definition 3.3

Linear operators $L_1:\text {h}_2\rightarrow \ell _2$ and $L_2:\text {h}_2\otimes \text {h}_2\rightarrow \ell _2$ are defined as follows

$$\begin{aligned} L_1[\hat{f}]&=\begin{pmatrix}-\sum _{l'=0}^\infty \sum _{m'=-l'}^{l'} \hat{f}_{l',m'}\langle \mathbf {y}_{l',m'}^{+},\mathbf {y}_{l,m}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ -\sum _{l'=0}^\infty \sum _{m'=-l'}^{l'} \hat{f}_{l',m'}\langle \mathbf {y}_{l',m'}^{+},\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \end{pmatrix}_{l\in \mathbb {N}_0,m=-l,\ldots ,l}, \\ L_2[\hat{g},\hat{h}]&=\begin{pmatrix}\sum _{l'=0}^\infty \sum _{m'=-l'}^{l'} \hat{g}_{l',m'}\langle \mathbf {y}_{l',m'}^{-},\mathbf {y}_{l,m}^+\rangle _{\mathbb {S}{\setminus }\Gamma }+\hat{h}_{l',m'}\langle \mathbf {y}_{l',m'}^{df},\mathbf {y}_{l,m}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ \sum _{l'=0}^\infty \sum _{m'=-l'}^{l'} \hat{g}_{l',m'}\langle \mathbf {y}_{l',m'}^{-},\mathbf {y}_{lm}^-\rangle _{\mathbb {S}{\setminus }\Gamma }+\hat{h}_{l',m'}\langle \mathbf {y}_{l',m'}^{df},\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \end{pmatrix}_{l\in \mathbb {N}_0,m=-l,\ldots ,l}. \end{aligned}$$

It is to note that we always need to assume $\hat{f}_{0,0}=\hat{h}_{0,0}=0$ in order to avoid problems due to the circumstance that we have defined $\mathbf {y}^+_{0,0}=\mathbf {y}^{df}_{0,0}=0$.

The next corollary is a direct consequence of the previous setup and Theorem 3.1.

Corollary 3.4

If $\mathbf {f}\in \mathbf {H}_2$ satisfies $L_1[\hat{\mathbf {f}}_+]=L_2[\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df}]$, then the knowledge of $\mathbf {f}_+$ uniquely determines $\mathbf {f}_-$. Furthermore, if $\mathbf {f}\in \mathbf {H}_2$ is spatially localized in $\Gamma$, then $L_1[\hat{\mathbf {f}}_+]=L_2[\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df}]$ is satisfied.

4 Numerical implementation and examples

The assertion of Corollary 3.4 justifies the computation of $\mathbf {f}_-$ from knowledge of $\mathbf {f}_+$ by solving $L_1[\hat{\mathbf {f}}_+]=L_2[\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df}]$ for $\hat{\mathbf {f}}_-$ and $\hat{\mathbf {f}}_{df}$. The contribution $\mathbf {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: Find $\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df}\in \text {h}_k$, $k\ge 2$, such that

$$\begin{aligned} (\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df})=\underset{\hat{g},\hat{h}\in \text {h}_2}{\text {argmin }}\Vert L_1[\hat{\mathbf {f}}_+]-L_2[\hat{h},\hat{g}]\Vert _{\ell _2}^2+\lambda \Vert (\hat{g},\hat{h})\Vert _{\text {h}_k}^2, \end{aligned}$$

(4.1)

for some fixed regularization parameter $\lambda >0$. The numerical setup for solving (4.1) is described in the next section.

4.1 Bandlimited numerical setup

We choose a bandlimit N and restrict ourselves to finding minimizers of (4.1) among all functions which satisfy that bandlimit. In other words, we are looking for $\hat{\mathbf {f}}_-,\hat{\mathbf {f}}_{df}\in \text {h}_2$ with $\hat{f}_{l,m}^-=\hat{f}_{l,m}^{df}=0$ for all $l\ge N+1$, $m=-l,\ldots ,l$. We abbreviate these solutions by (finite) vectors $\hat{\mathbf {f}}_-^N,\hat{\mathbf {f}}_{df}^N\in \mathbb {C}^{(N+1)^2}$ that contain all Fourier coefficients $\hat{f}_{l,m}^-$, $\hat{f}_{l,m}^{df}$ with $l=0,\ldots ,N$, $m=-l,\ldots ,l$. Furthermore, we define matrices $\varvec{\mathcal {L}}_1\in \mathbb {C}^{2(N+1)^2\times (N+1)^2}$, $\varvec{\mathcal {L}}_2\in \mathbb {C}^{2(N+1)^2\times 2(N+1)^2}$ by

$$\begin{aligned} \varvec{\mathcal {L}}_1&=\begin{pmatrix} -\langle \mathbf {y}_{0,0}^+,\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,-1}^+,\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,0}^+,\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} -\langle \mathbf {y}_{N,N}^+,\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ -\langle \mathbf {y}_{0,0}^+,\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,-1}^+,\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,0}^+,\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} -\langle \mathbf {y}_{N,N}^+,\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \\ -\langle \mathbf {y}_{0,0}^+,\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,-1}^+,\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,0}^+,\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} -\langle \mathbf {y}_{N,N}^+,\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ -\langle \mathbf {y}_{0,0}^+,\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,-1}^+,\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,0}^+,\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} -\langle \mathbf {y}_{N,N}^+,\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \\ \vdots &{}\vdots &{}\vdots &{}\ddots &{}\vdots \\ -\langle \mathbf {y}_{0,0}^+,\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,-1}^+,\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,0}^+,\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} -\langle \mathbf {y}_{N,N}^+,\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ -\langle \mathbf {y}_{0,0}^+,\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,-1}^+,\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} -\langle \mathbf {y}_{1,0}^+,\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} -\langle \mathbf {y}_{N,N}^+,\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \end{pmatrix}, \end{aligned}$$

(4.2)

$$\begin{aligned} \varvec{\mathcal {L}}_2&=\begin{pmatrix} \langle \mathbf {y}_{0,0}^-,\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^-,\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \langle \mathbf {y}_{0,0}^{df},\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^{df},\mathbf {y}_{0,0}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ \langle \mathbf {y}_{0,0}^{-},\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^{-},\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \langle \mathbf {y}_{0,0}^{df},\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^{df},\mathbf {y}_{0,0}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \\ \langle \mathbf {y}_{0,0}^-,\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^-,\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \langle \mathbf {y}_{0,0}^{df},\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^{df},\mathbf {y}_{1,-1}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ \langle \mathbf {y}_{0,0}^-,\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^-,\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \langle \mathbf {y}_{0,0}^{df},\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^{df},\mathbf {y}_{1,-1}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \\ \vdots &{}\ddots &{}\vdots &{}\vdots &{}\ddots &{}\vdots \\ \langle \mathbf {y}_{0,0}^-,\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^-,\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \langle \mathbf {y}_{0,0}^{df},\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^{df},\mathbf {y}_{N,N}^+\rangle _{\mathbb {S}{\setminus }\Gamma } \\ \langle \mathbf {y}_{0,0}^-,\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^-,\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \langle \mathbf {y}_{0,0}^{df},\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } &{} \cdots &{} \langle \mathbf {y}_{N,N}^{df},\mathbf {y}_{N,N}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \end{pmatrix}. \end{aligned}$$

(4.3)

Within this bandlimited setup, the minimzers $\hat{\mathbf {f}}_-^N$, $\hat{\mathbf {f}}_{df}^N$ of (4.1) are obtained as solutions to the normal equations

$$\begin{aligned} \left( \varvec{\mathcal {L}}_2^T\varvec{\mathcal {L}}_2+\lambda \varvec{\mathcal {R}}^T\varvec{\mathcal {R}}\right) \begin{pmatrix}\hat{\mathbf {f}}_-^N\\ \hat{\mathbf {f}}_{df}^N\end{pmatrix}=\varvec{\mathcal {L}}_2^T\varvec{\mathcal {L}}_1\hat{\mathbf {f}}_+^N, \end{aligned}$$

(4.4)

where the (finite) vector $\hat{\mathbf {f}}_+^N\in \mathbb {C}^{(N+1)^2}$ on the right hand side contains all Fourier coefficients $\hat{f}^+_{l,m}$, with $l=0,\ldots ,N$, $m=-l,\ldots ,l$, of the known function $\mathbf {f}_+$. Since we chose the regularization term in (4.1) to penalize the $\Vert \cdot \Vert _{\text {h}_k}$-norm, the matrix $\varvec{\mathcal {R}}\in \mathbb {R}^{2(N+1)^2\times 2(N+1)^2}$ is a diagonal matrix with diagonal entries $R_{i,i}=(l+\frac{1}{2})^k$, for $i=l^2+m+l+1$ or $i=(N+1)^2+l^2+m+l+1$ and $l=0,\ldots , N$, $m=-l,\ldots , l$.

For the sake of further simplicity, we focus on the case that the subset $\Gamma \subset \mathbb {S}$, in which we assume $\mathbf {f}$ to be spatially localized, is a spherical cap with the North Pole $\mathbf {e}_3=(0,0,1)^T$ as center and with (angular) radius $\Theta _0\in (0,\pi )$, i.e., $\Gamma =\{x=x(\theta ,\phi )\in \mathbb {S}:\phi \in [0,2\pi ),\theta \in [0,\Theta _0)\}=\{x\in \mathbb {S}:x\cdot \mathbf {e}_3>\cos (\Theta _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 $\varvec{\mathcal {L}}_1$ and $\varvec{\mathcal {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 $\mathbf {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 Sect. 4.2 suggest that the error is minor if the Fourier coefficients ${\hat{\mathbf {f}}}^+_{l,m}$ decay sufficiently fast while it can be quite significant if such a decay is not given.

4.2 Examples

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

4.2.1 Example 1

First we choose $\mathbf {f}$ to be localized in the entire northern hemisphere, i.e., a spherical cap with angular radius $\Theta _0=\frac{\pi }{2}$ and the North Pole $\mathbf {e}_3=(0,0,1)^T$ as center. It is of the form

$$\begin{aligned} \mathbf {f}(x)&=Q(x)\left( 3(x\cdot \mathbf {d})x-\mathbf {d}\right) , \nonumber \\ Q(x)&=Q(\theta ,\phi )=\left\{ \begin{array}{ll}(\cos (\theta )-0.1)^3\,(\cos (\theta )-1)^2\,(\phi -2\pi )^3\,\phi ^3\,\sin (2\phi ),&{} \text {if }\cos (\theta )>0.1,\\ 0,&{}\text {else},\end{array}\right. \nonumber \\ \end{aligned}$$

(4.5)

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 $\mathbf {d}=(0,0.6,0.8)^T$. The radial part of the function $\mathbf {f}$ and its contributions $\mathbf {f}_+$ and $\mathbf {f}_-$ are indicated in Fig. 1. In fact, one can see that although $\mathbf {f}$ is spatially localized in the northern hemisphere, this does not need to hold true for its contributions $\mathbf {f}_+$ and $\mathbf {f}_-$. Reconstructions $\mathbf {f}_-^N$ based on knowledge of the Fourier coefficients of $\mathbf {f}_+$ up to degree $N=100$ and for smoothness parameter $k=2$ are shown in Fig. 2. For the choice $\lambda =0.5\cdot 10^{-15}$ we obtained the best reconstruction of the actual contribution $\mathbf {f}_-$ among the tested regularization parameters. Sightly overregularized and underregularized results are shown for $\lambda =10^{-12}$ and $\lambda =0.5\cdot 10^{-17}$.

Fig. 1

Radial component of the (unknown) function $\mathbf {f}$ from (4.5) (left), of its (known) contribution $\mathbf {f}_+$ (center), and of its (unknown) contribution $\mathbf {f}_-$ (right)

Fig. 2

Radial component of the (unknown) contribution $\mathbf {f}_-$ of $\mathbf {f}$ from (4.5) (bottom) and the reconstructed contribution $\mathbf {f}_-^N$ for $N=100$ and regularization parameters $\lambda =10^{-12}$ (top left), $\lambda =0.5\cdot 10^{-15}$ (top center), and $\lambda =0.5\cdot 10^{-17}$ (top right)

Fig. 3

Radial component of the (unknown) contribution $\mathbf {f}_-$ of $\mathbf {f}$ from (4.5) (bottom) and the reconstructed contribution $\mathbf {f}_-^{N,\varepsilon }$ for $N=100$, smoothness parameter $k=2$ and regularization parameters $\lambda =10^{-11}$ (top left), $k=2$ and $\lambda =10^{-12}$ (top center), and $k=4$ and $\lambda =10^{-15}$ (top right)

Figure 3 shows the results for the same setup as before, but with noisy input data $\hat{\mathbf {f}}_+^{N,\varepsilon }=\hat{\mathbf {f}}_+^N+\hat{\varvec{\varepsilon }}^N$. By $\hat{\varvec{\varepsilon }}^N\in \mathbb {C}^{(N+1)^2}$ we denote a random vector where each component is normally distributed with mean zero and and variance one, scaled to an overall noise level ${\Vert \hat{\varvec{\varepsilon }}^N\Vert }/{\Vert \hat{\mathbf {f}}_+^N\Vert }=2.5\times 10^{-3}$. Although the overall noise level is quite small, Fig. 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 reconstructed $\hat{\mathbf {f}}_-^{N,\varepsilon }$ in Fig. 3: For smoothness parameter $k=2$, the general features of $\mathbf {f}_-$ could be reconstructed, but their amplitude is generally underestimated (for $\lambda =10^{-11}$) or the reconstruction is underregularized and noisy features are introduced (for $\lambda =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.

Fig. 4

Radial component of the (unknown) contribution $\mathbf {f}_-$ of $\mathbf {f}$ from (4.5) (bottom) and the best reconstructed contribution $\mathbf {f}_-^{N}$ for $N=50$ without noise (top left), for $N=100$ without noise (top center), and for $N=100$ with noise (top right)

Fig. 5

Fourier spectrum representation of the undisturbed $\hat{\mathbf {f}}_+^{N}$ from Example 1 (left), of the disturbed input data $\hat{\mathbf {f}}_+^{N,\varepsilon }$ (center), and of the relative quotient of the noise $\varvec{\varepsilon }^N$ and the undisturbed $\hat{\mathbf {f}}_+^N$ (right)

Additionally, we varied the truncation degree N for the bandlimitation. The results are indicated in Fig. 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, Fig. 6 indicates that the chosen function $\mathbf {f}$ in this example has a rapidly decreasing power spectrum. Thus, one would not expect a major effect in the first place.

4.2.2 Example 2

We choose the following function that reveals a stronger spatial localization than the first example:

$$\begin{aligned} \mathbf {f}(x)=\mathbf {f}(\theta ,\phi )=\left\{ \begin{array}{ll}\begin{pmatrix}(t-0.3)^3\,(1-t)^3\,\sin (2(\phi -\frac{\pi }{2}))\,\sin (3(\phi -\frac{3\pi }{2})) \\ (t-0.3)^3\,(1-t)^3\,\sin (2\pi t)\,\sin (\phi -\frac{\pi }{2})\,\sin (2(\phi -\frac{3\pi }{2})) \\ (t-0.3)^3\,(1-t)^3\,\sin (2\pi t)\,\sin (\phi -\frac{\pi }{2})\,\sin (2(\phi -\frac{3\pi }{2}))\,e^{-2|\phi -\pi |}\end{pmatrix},\\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \text {if }t=\cos (\theta )>0.3,\,|\phi -\pi |<\frac{\pi }{2}, \\ 0,\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \text {else}.\end{array}\right. \nonumber \\ \end{aligned}$$

(4.6)

The function $\mathbf {f}$ and its contributions $\mathbf {f}_+$ and $\mathbf {f}_-$ are illustrated in Fig. 7. Reconstructions $\mathbf {f}_-^N$ for $N=100$ and various choices of smoothing parameters $k=2,3,4$ are provided in Fig. 8 (for each smoothing parameter k, the best result among the tested regularization parameters $\lambda$ is shown). Overall, it can be said that the results provide a reasonable approximation of the actual contribution $\mathbf {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 $\mathbf {f}_-$ seems worse than in Example 1, which is most likely due to the slower decay of the power spectrum towards higher degrees (compare Fig. 6).

In Fig. 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 $\lambda$ 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 Fig. 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.

Fig. 6

Scaled logarithmic power spectrum of the function $\mathbf {f}_+$ from Example 1 (blue line) and from Example 2 (red line)

Fig. 7

Radial component of the (unknown) function $\mathbf {f}$ from (4.6) (left), of its (known) contribution $\mathbf {f}_+$ (center), and of its (unknown) contribution $\mathbf {f}_-$ (right)

Fig. 8

Radial component of the (unknown) contribution $\mathbf {f}_-$ of $\mathbf {f}$ from (4.6) (bottom) and the reconstructed contribution $\mathbf {f}_-^N$ for $N=100$, smoothness parameter $k=2$ and regularization parameter $\lambda =10^{-14}$ (top left), $k=3$ and $\lambda =0.5\times 10^{-17}$ (top center), and $k=4$ and $\lambda =10^{-18}$ (top right)

Fig. 9

Radial component of the (unknown) contribution $\mathbf {f}_-$ of $\mathbf {f}$ from (4.6) (bottom) and the best reconstructed contribution $\mathbf {f}_-^{N}$ for $N=100$ (top left), for $N=140$ (top center), and for $N=200$ (top right)

5 Conclusion

In this paper, we have derived a method on how to numerically connect the contributions $\mathbf {f}_+$ and $\mathbf {f}_-$ of a general vector field $\mathbf {f}$ under the a priori assumption that $\mathbf {f}$ is spatially localized in a subdomain of the sphere. The examples in Sect. 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 Lesur and Vervelidou (2020) and Vervelidou and Lesur (2018), that connect $\mathbf {f}_+$ and $\mathbf {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 $\mathbf {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 $\mathbf {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.

Appendix A: Computation of matrix entries

In order to obtain the entries of the matrices $\varvec{\mathcal {L}}_1$, $\varvec{\mathcal {L}}_2$ when $\Gamma$ is a spherical cap with center $\mathbf {e}_3=(0,0,1)^T$ and angular radius $\Theta _0\in (0,\pi )$, we exemplarily perform the explicit computations for the terms $\langle \mathbf {y}_{l',m'}^+,\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma }$, $\langle \mathbf {y}_{l',m'}^{df},\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma }$. All other terms can be obtained similarly. It is helpful to consider the following alternative representation of the vector spherical harmonics (2.8)–(2.10) (e.g., Freeden and Gerhards 2012; Freeden and Schreiner 2009):

$$\begin{aligned} \mathbf {y}_{l,m}^-&=-\frac{1}{\sqrt{(l+1)(2l+1)}}\left( \hat{\mathbf {n}}(l+1)Y_{l,m}-\nabla _\mathbb {S}Y_{l,m}\right) , \\ \mathbf {y}_{l,m}^+&=\frac{1}{\sqrt{l(2l+1)}}\left( \hat{\mathbf {n}}\,l\,Y_{l,m}+\nabla _\mathbb {S}Y_{l,m}\right) , \\ \mathbf {y}_{l,m}^{df}&=-\frac{i}{\sqrt{l(l+1)}}\text {L}_\mathbb {S}Y_{l,m}. \end{aligned}$$

Using Green’s formula, and denoting the Laplace-Beltrami operator by $\Delta _\mathbb {S}=\nabla _\mathbb {S}\cdot \nabla _\mathbb {S}$ and the normal derivative with respect to the boundary $\partial \Gamma$ by $\frac{\partial }{\partial {\nu }}$, we obtain

$$\begin{aligned}&\langle \mathbf {y}_{l',m'}^+,\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \\&\quad =-\frac{1}{4\pi \sqrt{(l+1)(2l+1)l'(2l'+1)}}\int _{\mathbb {S}{\setminus }\Gamma }l'(l+1)Y_{l',m'}(y)Y_{l,m}(y)^*\\&\qquad -\nabla _\mathbb {S}Y_{l',m'}(y)\cdot \nabla _\mathbb {S}Y_{l,m}^*(y)\,d\omega (y) \\&\quad =-\frac{1}{4\pi \sqrt{(l+1)(2l+1)l'(2l'+1)}}\int _{\mathbb {S}{\setminus }\Gamma }l'(l+1)Y_{l',m'}(y)Y_{l,m}(y)^*\\&\qquad +Y_{l',m'}(y)\Delta _\mathbb {S}Y_{l,m}^*(y)\,d\omega (y) \\&\qquad +\frac{1}{4\pi \sqrt{(l+1)(2l+1)l'(2l'+1)}}\int _{\partial \Gamma }Y_{l',m'}(y)\frac{\partial }{\partial {\nu }}Y_{l,m}(y)^*\,d\sigma (y) \\&\quad =\frac{(l-l')(l+1)}{4\pi \sqrt{(l+1)(2l+1)l'(2l'+1)}}\int _{\mathbb {S}{\setminus }\Gamma }Y_{l',m'}(y)Y_{l,m}(y)^*\,d\omega (y) \\&\qquad +\frac{1}{4\pi \sqrt{(l+1)(2l+1)l'(2l'+1)}}\int _{\partial \Gamma }Y_{l',m'}(y)\frac{\partial }{\partial {\nu }}Y_{l,m}(y)^*\,d\sigma (y), \\&\langle \mathbf {y}_{l',m'}^{df},\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma } \\&\quad =-\frac{i}{4\pi \sqrt{(l+1)(2l+1)l'(l'+1)}}\int _{\mathbb {S}{\setminus }\Gamma }\text {L}_\mathbb {S}Y_{l',m'}(y)\cdot \nabla _\mathbb {S}Y_{l,m}^*(y)\,d\omega (y) \\&\quad =\frac{i}{4\pi \sqrt{(l+1)(2l+1)l'(l'+1)}}\int _{\mathbb {S}{\setminus }\Gamma }Y_{l',m'}(y)\text {L}_\mathbb {S}\cdot \nabla _\mathbb {S}Y_{l,m}^*(y)\,d\omega (y) \\&\qquad -\frac{i}{4\pi \sqrt{(l+1)(2l+1)l'(l'+1)}}\int _{\partial \Gamma }Y_{l',m'}(y)\tau (y)\cdot \nabla _\mathbb {S}Y_{l,m}(y)^*\,d\sigma (y) \\&\quad =-\frac{i}{4\pi \sqrt{(l+1)(2l+1)l'(l'+1)}}\int _{\partial \Gamma }Y_{l',m'}(y)\tau (y)\cdot \nabla _\mathbb {S}Y_{l,m}(y)^*\,d\sigma (y). \end{aligned}$$

The appearing integrals over $\mathbb {S}{\setminus }\Gamma$ can be simplified to

$$\begin{aligned}&\gamma _{l',m',l,m}\\&\quad =\int _{\mathbb {S}{\setminus }\Gamma } Y_{l',m'}(y) Y_{l,m}(y)^*\,d\omega (y) \\&\quad =\int _0^{2\pi }\int _{\Theta _0}^\pi \sqrt{\frac{2l+1}{2-\delta _{m,0}}}P_l^m(\cos (\theta )) e^{-im\phi }\sqrt{\frac{2l'+1}{2-\delta _{m',0}}}P_{l'}^{m'}(\cos (\theta )) e^{im'\phi }\sin (\theta )\,d\theta \,d\phi \nonumber \\&\quad =\left\{ \begin{array}{ll} 0,&{}\quad \text { if }m\not =m', \\ \frac{2\pi \sqrt{(2l+1)(2l'+1)}}{2-\delta _{m,0}}\int _{-1}^{\cos (\Theta _0)}P_l^m(t)P_{l'}^{m}(t)\,dt,&{}\quad \text { else.} \end{array}\right. \end{aligned}$$

For the remaining boundary integrals over $\partial \Gamma$, we first observe that $\frac{\partial }{\partial {\nu }}=\nu \cdot \nabla _\mathbb {S}$ and that the unit normal vector $\nu$ (pointing into the exterior of $\mathbb {S}{\setminus }\Gamma$) takes the form $\nu (y)=\frac{\mathbf {e}_3-(\mathbf {e}_3\cdot y)y}{\sin (\Theta _0)}$ if $\Gamma$ is a spherical cap with center $\mathbf {e}_3=(0,0,1)^T$ and angular radius $\Theta _0\in (0,\pi )$. The unit tangential vector $\tau$ takes the form $\tau (y)=-\frac{y\times \mathbf {e}_3}{\sin (\Theta _0)}$. By use of recursion relations for the associated Legendre polynomials (e.g., Freeden and Schreiner 2009), we obtain that

$$\begin{aligned} \kappa _{l',m',l,m}&=\int _{\partial \Gamma }Y_{l',m'}(y)\frac{\partial }{\partial \nu }Y_{l,m}(y)^*\,d\sigma (y) \\&=\int _0^{2\pi }\sqrt{\frac{2l'+1}{2-\delta _{m',0}}}e^{im'\phi }P_{l'}^{m'}(\cos (\Theta _0))\frac{1-t^2}{\sin (\Theta _0)}\frac{\partial }{\partial t}\sqrt{\frac{2l+1}{2-\delta _{m,0}}} e^{-im\phi }P_{l}^{m}(t)\bigg |_{t=\cos (\Theta _0)}d\phi \nonumber \\&=\left\{ \begin{array}{ll} 0,&{}\quad \text { if }m\not =m', \\ \frac{2\pi \sqrt{(2l+1)(2l'+1)}}{(2-\delta _{m,0})\sin (\Theta _0)}P_{l'}^{m}(\cos (\Theta _0))(1-t^2)\frac{\partial }{\partial t}P_{l}^{m}(t)\Big |_{t=\cos (\Theta _0)},&{}\quad \text { else,} \end{array}\right. \\&=\left\{ \begin{array}{ll} 0,&{}\quad \text { if }m\not =m', \\ \frac{2\pi \sqrt{(2l+1)(2l'+1)}}{(2-\delta _{m,0})\sin (\Theta _0)}P_{l'}^m(\cos (\Theta _0))\Big ((l+1)\cos (\Theta _0)P_{l}^{m}(\cos (\Theta _0)) \\ -\sqrt{(l-m+1)(l+m+1)}P_{l+1}^{m}(\cos (\Theta _0))\Big ),&{}\quad \text { else,} \end{array}\right. \\ \mu _{l',m',l,m}&=\int _{\partial \Gamma }Y_{l',m'}(y)\tau (y)\cdot \nabla _\mathbb {S}Y_{l,m}(y)^*\,d\sigma (y) \\ {}&=\int _0^{2\pi }\sqrt{\frac{2l'+1}{2-\delta _{m',0}}}e^{im'\phi }P_{l'}^{m'}(\cos (\Theta _0))\frac{1}{\sin (\Theta _0)}\frac{\partial }{\partial \phi }\sqrt{\frac{2l+1}{2-\delta _{m,0}}} e^{-im\phi }P_{l}^{m}(\cos (\Theta _0))d\phi \nonumber \\ {}&=\left\{ \begin{array}{ll} 0,&{}\quad \text { if }m\not =m', \\ -\frac{2im\pi \sqrt{(2l+1)(2l'+1)}}{(2-\delta _{m,0})\sin (\Theta _0)}P_{l'}^m(\cos (\Theta _0))P_{l}^{m}(\cos (\Theta _0)),&{}\quad \text { else.} \end{array}\right. \end{aligned}$$

Summarizing, the entries of the matrices $\varvec{\mathcal {L}}_1$ and $\varvec{\mathcal {L}}_2$ look as follows:

$$\begin{aligned} \langle \mathbf {y}_{l',m'}^+,\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma }&=\frac{(l-l')(l+1)\gamma _{l',m',l,m}+\kappa _{l',m',l,m}}{4\pi \sqrt{(l+1)(2l+1)l'(2l'+1)}}, \\ \langle \mathbf {y}_{l',m'}^-,\mathbf {y}_{l,m}^+\rangle _{\mathbb {S}{\setminus }\Gamma }&=\frac{(l-l')\,l\,\gamma _{l',m',l,m}+\kappa _{l',m',l,m}}{4\pi \sqrt{(l'+1)(2l'+1)l(2l+1)}},\\ \langle \mathbf {y}_{l',m'}^+,\mathbf {y}_{l,m}^+\rangle _{\mathbb {S}{\setminus }\Gamma }&=\frac{(l+l'+1)\,l\,\gamma _{l',m',l,m}+\kappa _{l',m',l,m}}{4\pi \sqrt{l(2l+1)l'(2l'+1)}}, \\ \langle \mathbf {y}_{l',m'}^-,\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma }&=\frac{(l+l'+1)(l+1)\gamma _{l',m',l,m}+\kappa _{l',m',l,m}}{4\pi \sqrt{(l+1)(2l+1)(l'+1)(2l'+1)}}, \\ \langle \mathbf {y}_{l',m'}^{df},\mathbf {y}_{l,m}^-\rangle _{\mathbb {S}{\setminus }\Gamma }&=-\frac{i\,\mu _{l',m',l,m}}{4\pi \sqrt{(l+1)(2l+1)l'(l'+1)}}, \\\langle \mathbf {y}_{l',m'}^{df},\mathbf {y}_{l,m}^+\rangle _{\mathbb {S}{\setminus }\Gamma }&=-\frac{i\,\mu _{l',m',l,m}}{4\pi \sqrt{l(2l+1)l'(l'+1)}}. \end{aligned}$$