# A new proposal for spherical cap harmonic modelling

E. Thébault, J. J. Schott, M. Mandea, J. P. Hoffbeck

► To cite this version:

E. Thébault, J. J. Schott, M. Mandea, J. P. Hoffbeck. A new proposal for spherical cap harmonic modelling. Geophysical Journal International, 2004, 159, pp.83-103. 10.1111/j.1365-246X.2004.02361.x. insu-03600274

HAL Id: insu-03600274

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

Submitted on 7 Mar 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A new proposal for spherical cap harmonic modelling

E. Thébault,<fnref n="1" />,* J. J. Schott,<fnref n="1" /> M. Mandea<fnref n="2" /> and J. P. Hoffbeck<fnref n="3" />

Accepted 2004 May 20. Received 2004 May 11; in original form 2003 June 27

## SUMMARY

The geomagnetic field above the surface of the Earth in the current-free region may be expressed as the gradient of a scalar potential solving Laplace's equation. For regions with a fairly dense coverage of data at different altitudes, a regional model ought to offer a better spatial resolution of the regional field over the volume under study than a global field expanded in spherical harmonics (SH). The spherical cap harmonics analysis (SCHA) is an attractive regional modelling tool having close relationship with global SH. With the SCHA adopted so far, difficulties arise in upward continuation and in establishing a relationship between global and local Gauss coefficients. Such a relationship would be useful, for instance, for introducing prior constraint on an inverse problem dealing with the estimation of local Gauss coefficients based upon a local data set. In this paper, we show that these difficulties are overcome if the SCHA modelling is formulated as a boundary value (BV) problem in a cone bounded radially by the surface of the Earth and an upper surface suitable for satellite data, and bounded laterally in order to encompass a specific region of study. Although the example worked out here applies only to a limited class of fields, which verifies some special flux condition, the ideas behind this formalism are quite general and should offer a new way of processing data in a bounded region of space.

**Key words:** magnetic field, regional modelling, spherical cap harmonic analysis.

## 1 INTRODUCTION

Regional modelling is in theory a powerful method for detailed description of potential fields over areas where an appropriate dense set of data is available. The advent of satellite magnetic measurements brought a major breakthrough for the reconstruction of a regional magnetic field at short wavelengths. For a limited portion of the Earth, spherical harmonics (SH), well known for global modelling, are no longer suitable, because they are no longer orthogonal over the restricted area. Common techniques, like polynomial modelling in latitude and longitude or rectangular harmonic analysis (Alldredge 1981), have been used successfully before the availability of satellite data but the resulting models could not be properly upward or downward continued (Haines 1990). The spherical cap harmonic analysis (SCHA) proposed by Haines (1985) is an attractive regional modelling. Its formalism looks like a natural extension of the spherical harmonic analysis. The method is claimed to be valid over any spherical cap at any altitude above the surface of the Earth. On the basis of these assertions, SCHA has been widely used for crustal anomalies modelling (De Santis et al. 1989), field modelling (Hwang & Chen 1997), or even regional secular variation modelling (Korte & Haak 2000). Nevertheless, SCHA practitioners encounter two kind of difficulties. First, it is broadly admitted that its behaviour goes worst as the cap goes smaller. Numerical experiments demonstrate that the convergence of SCHA is dramatically slow for relatively small caps and an insufficient expansion produces unrealistic oscillations in the results when interpolating on a dense grid. When dealing with a small cap, a larger expansion is required, most of the time larger than permitted by the number of available data. This is especially true for the regional secular variation and explains why it is meaningless to model it over small caps with SCHA (De Santis et al. 1997). However, the second and most serious problem, stressed by De Santis & Falcone (1995) and also mentioned by Lowes (Langel & Hintz 1998), comes from the failure of modelling correctly the radial dependence. None of the basis functions proposed by Haines is appropriate for adjusting fields decreasing as $r^{-n}$ where n is a real integer. Therefore, it is impossible, using this formalism, to solve properly a boundary value (BV) problem where the boundary conditions are a global field derived from an SH model (either a complete model or a model containing only the highest degrees in the expansion). In this paper, we propose a new SCHA expansion that removes this drawback. The proper resolution of the BV problem on the conical domain, depicted on Fig. 1, enables us to establish a relationship between global and local Gauss coefficients. We finally succeed in reconstructing the corresponding local field by these relationships at any altitude and colatitude within Ω. However, as a

<footnote n="1">¹Ecole et Observatoire des Sciences de la Terre, 5 rue rené Descartes, 67000 Strasbourg, France. E-mail: Erwan.Thebault@eost.u-strasbg.fr</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, 4 place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote n="3">³E.N.S Cachan, Université de Rennes I, 35042 Rennes Cedex, France</footnote>

<footnote>*Now at: GeoForschungs Zentrum of Potsdam, Telegrafenberg, 14473 Potsdam, Germany. E-mail: Erwan.Thebault@gfz-potsdam.de.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Definition of the domain of study Ω bounded by the terrestrial surface $r = a$ and the upper surface $r = ae^S$ with $S \neq 0$. The real scalar $S$ is chosen such that the data recorded by satellites lie inside the volume. $\partial\Omega_{\theta_0}$ is the boundary $\theta = \theta_0$, and $\partial\Omega_a$ and $\partial\Omega_b$ respectively are the lower and upper caps.</figcaption></figure>

result of the adopted approach, which makes use of information on the field and not on the potential, the modelling proposed here is suitable only for a restricted class of fields.

## 2 DEFINITION OF THE PROBLEM

The domain Ω, on which potentials and fields are dealt with, is depicted on Fig. 1. Although it is not strictly speaking correct, we will call this domain a conical domain or even a cone for the sake of brevity. Assuming that the sources are located outside the cone and that no current crosses the surfaces defined by $\partial_{\theta_0}\Omega$, $\partial_a\Omega$ and $\partial_b\Omega$, the field is curl- and divergence-free. Hence, the main problem describing the magnetic potential must satisfy Laplace's equation under some boundary conditions. The boundary conditions will constrain the expression and properties of the harmonic functions building up a base in the space of the functions that verify Laplace's equation on Ω. Classical boundary conditions are of Dirichlet, Neumann or mixed type. In the present study, we consider Neumann-like conditions as they seem to be the most natural ones (because they deal with the field itself, which is the measurable quantity in practice). Let us write these conditions that lead to the following BV problem:

$$ \Delta V_c = 0, \tag{1} $$

$$ \left. \frac{\partial V_c}{\partial \nu} \right|_{\partial \theta_0 \Omega} = \left. \frac{1}{r} \frac{\partial V_c}{\partial \theta} \right|_{\partial \theta_0 \Omega} = F(r, \varphi), \tag{2} $$

$$ \left. \frac{\partial V_c}{\partial \nu} \right|_{\partial a \Omega} = - \left. \frac{\partial V_c}{\partial r} \right|_{\partial a \Omega} = -G_1(\theta, \varphi), \tag{3} $$

$$ \left. \frac{\partial V_c}{\partial \nu} \right|_{\partial b \Omega} = \left. \frac{\partial V_c}{\partial r} \right|_{\partial b \Omega} = G_2(\theta, \varphi), \tag{4} $$

where $\partial/\partial\nu$ stands for the outer normal derivative to the boundary $\partial\Omega$ of the conical domain. The boundary $\partial\Omega$ is split up into three parts: $\partial_{\theta_0}\Omega$ is the lateral boundary $\{\theta = \theta_0, a \le r \le b\}$, $\partial_a\Omega$ the lower cap $\{\theta \le \theta_0, r = a\}$ and $\partial_b\Omega$ the upper cap $\{\theta \le \theta_0, r = b\}$ (see Fig. 1).
The divergence-free condition imposes a null flux through the surfaces:

$$ \Phi = \int \int_{\partial\theta\Omega} F ds - \int \int_{\partial\phi\Omega} G_1 ds + \int \int_{\partial\phi\Omega} G_2 ds = 0. \quad (5) $$

To solve the boundary value problem, we may decompose it into two subproblems having homogeneous Neumann boundary conditions. The decomposition, which leads to two independent problems, is straightforwardly guided by the initial problem stated in eqs (1) to (4).

Let us define two potential $V_1$ and $V_2$ solving separately problems 1 and 2 settled below. The BV problem 1 may be defined by the following equations:

$$ \Delta V_1 = 0, \quad (6a) $$

$$ \left. \frac{\partial V_1}{\partial \theta} \right|_{\partial \Omega_{\theta}} = F(r, \varphi), \quad (6b) $$

$$ \left. \frac{\partial V_1}{\partial r} \right|_{\partial \Omega_a} = 0, \quad (6c) $$

$$ \left. \frac{\partial V_1}{\partial r} \right|_{\partial \Omega_b} = 0. \quad (6d) $$

As, $\vec{di}\vec{v}\vec{B} = 0$, there can be no net flux into or out of the boundary $\partial\Omega$. From (6c) and (6d) it becomes:

$$ \int \int_{\partial \Omega_{\theta}} F ds = 0 \quad (7) $$

and problem 2 defined by:

$$ \Delta V_2 = 0, \quad (8a) $$

$$ \left. \frac{\partial V_2}{\partial \theta} \right|_{\partial \Omega_{\theta}} = 0, \quad (8b) $$

$$ \left. \frac{\partial V_2}{\partial r} \right|_{\partial \Omega_a} = G_1(\theta, \varphi), \quad (8c) $$

$$ \left. \frac{\partial V_2}{\partial r} \right|_{\partial \Omega_b} = G_2(\theta, \varphi), \quad (8d) $$

using the same reasoning as with problem 1 leads to the flux condition:

$$ - \int \int_{\partial \Omega_a} G_1 ds + \int \int_{\partial \Omega_b} G_2 ds = 0. \quad (9) $$

The splitting of the initial BV problem into subproblems 1 and 2 implicitly imposes additional flux conditions not included in the complete problem. Therefore, the functions $F$, $G_1$ and $G_2$ of the initial problem have to be modified in order to fulfil also the flux constraints of problems 1 and 2. This question will be dealt with in Section 5.1.

The solution of the problem 2 has already been given by Haines (1985). It is now widely known that this problem involves the associated Legendre function with real degrees $n$ and integer orders $m$. Because the derivative of the associated Legendre functions have been set to zero on the boundary $\theta_0$, in Haines's notation this would correspond to the $k-m$ even Legendre functions (see Section 4 for further details). For a complete development of $V_2$, we refer to Haines's theory (Haines 1985, 1988).

On the other hand, the solution of problem 1 is quite unusual in geomagnetism and requires more attention, therefore, the focus will be on its resolution. Clearly, $V_c$, the solution of the initial BV problem, is the sum of $V_1$ and $V_2$, provided that $F$, $G_1$ and $G_2$ satisfy the conditions (7) and (9).

### 3 SOLUTION FOR THE POTENTIAL $V_1$

Let us recall that the Laplacian in spherical coordinates is

$$ \Delta = \left( \frac{\partial^2}{\partial r^2} + \frac{2}{r} \frac{\partial}{\partial r} \right) + \frac{1}{r^2 \sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial}{\partial \theta} \right) + \frac{1}{r^2 \sin^2 \theta} \frac{\partial^2}{\partial \phi^2} $$

and that the usual way of solving $\Delta V = 0$ is to search for solutions of the form $V_1 = R(r)\Theta(\theta)\Phi(\varphi)$ where $r$ is the radial distance, $\theta$ the colatitude and $\varphi$ the longitude. This decomposition leads to the well-known triplet of ordinary differential equations:

$$ \frac{d}{dr} \left[ r^2 \frac{dR(r)}{dr} \right] = \kappa R(r), \quad (10a) $$
$$ \sin \theta \frac{1}{\Theta} \frac{d}{d\theta} \left( \sin \theta \frac{d\Theta}{d\theta} \right) + \kappa \sin^2 \theta = m^2, \qquad (10b) $$

$$ \frac{d^2\Phi}{d\varphi^2} + m^2\Phi = 0. \qquad (10c) $$

## 3.1 Radial function R

Let $\lambda = -\kappa$. Eq. (10a) may be written in the form:

$$ -\frac{d}{dr} \left[ r^2 \frac{dR(r)}{dr} \right] = \lambda R(r). \qquad (11) $$

This equation associated with the boundary conditions imposed in eqs (6c) and (6d) forms a Sturm-Liouville problem leading to the construction of an orthogonal basis on the interval $[a, b]$, for the inner product:

$$ \int_a^b f.g dr = \langle f, g \rangle_{[a,b]}. $$

It is worth stressing the following point. Usually, Sturm-Liouville problems arise from eqs (10b) and (10c), which, consequently, fix the values of $m$ and the eigenvalue $\lambda = n(n + 1)$. In the present case, as it will be shown below, the set of admissible eigenvalues $\lambda$ is no longer determined by eq. (10b), which does not define a Sturm-Liouville problem as a result of the absence of a boundary condition, but, to the contrary, by eq. (10a) associated with the boundary conditions mentioned above.

Eq. (11) belongs to the Euler class of ordinary differential equation (ODE). It may be classically transformed into a linear ODE by means of the change of independent variables (Zwillinger 1989) $r = ae^t$, where $a$ is an arbitrary constant, taken in our specific problem as the mean radius of the Earth:

$$ \frac{d^2 R(t)}{dt^2} + \frac{dR(t)}{dt} + \lambda R(t) = 0. $$

The solution of this linear ordinary differential equation depends on the value of its discriminant.

(i) If $1 - 4\lambda \neq 0$ then the solution, written in terms of $r$, is:

$$ R(r) = A_1 \left(\frac{r}{a}\right)^n + A_2 \left(\frac{a}{r}\right)^{n+1} \qquad (12) $$

with the eigenvalues $\lambda = -n(n + 1)$.

(ii) If $1 - 4\lambda = 0$ then the eigenvalue $\lambda = +1/4$ and thus the solution, in terms of $r$, is:

$$ R(r) = \left[ A_1 \log \left( \frac{r}{a} \right) + A_2 \right] \sqrt{\frac{a}{r}}. \qquad (13) $$

The factor constants $A_1$ and $A_2$ are in general complex. Their value depends on the chosen normalization. These solutions hold for any value of the degree $n$, including $n$ complex. The solution (12), which takes apparently the familiar form, is an instance where $n$ actually belongs to a complex set of numbers as will be shown in the next section.

### Boundary conditions

In the case $1 - 4\lambda \neq 0$, eq. (12) associated with the boundary conditions (6c) and (6d) leads to the following system:

$$ nA_1 - (n + 1)A_2 = 0, $$

$$ nA_1 \frac{b^{n-1}}{a^n} - (n + 1)A_2 \frac{a^{n+1}}{b^{n+2}} = 0. $$

To obtain a non-trivial zero solution in $A_1$ and $A_2$, the determinant must be zero. This condition is fulfilled for

(i) $n = 0, n = -1$, that is $\lambda = -\kappa = 0$. It is easily seen that in this case, $R(r) = R_0$ is a constant. Obviously, a constant function fits the boundary conditions.

(ii) $a^{2n+1} = b^{2n+1}$. This equality leads to a complex set of values $n$, taking the form $n = -\frac{1}{2} + i \frac{p\pi}{\log(b/a)}$ where $p \in \mathbb{N}^*$; using Euler's trigonometric formula and setting $b = ae^S$ to simplify the notation, the radial functions reduce to:

$$ R_p(r) = \sqrt{\frac{a}{r}} \left\{ A_1 \cos \left[ \frac{p\pi}{S} \log \left( \frac{r}{a} \right) \right] + A_2 \sin \left[ \frac{p\pi}{S} \log \left( \frac{r}{a} \right) \right] \right\}. \qquad (14) $$
Taking into account the boundary conditions we get the following relation between $A_1$ and $A_2$:

$$A_1 = \frac{2\pi p}{S} A_2.$$

The normalization of the relevant basis functions $R_0$ and $R_p(r)$, such that $\|R_p\|^2 = \|R_0\|^2 = a$, leads to the set of functions:

$$R_0 = \frac{1}{\sqrt{(e^S - 1)}}, \qquad (15)$$

$$R_p(r) = \frac{(2S)^{-1/2}}{\sqrt{\left(\frac{2\pi p}{S}\right)^2 + 1}} \left\{ \frac{2\pi p}{S} \cos\left[\frac{p\pi}{S} \log\left(\frac{r}{a}\right)\right] + \sin\left[\frac{\pi p}{S} \log\left(\frac{r}{a}\right)\right] \right\} \quad \text{with } p \in \mathbb{N}^*. \qquad (16)$$

$R_0$ and $R_p(r)$ being solutions of a Sturm-Liouville problem, are orthogonal and form a complete basis in the Hilbert space $L^2(a, b)$ for the inner product defined above. Thus, these functions oscillate around zero and have $p$ roots in the interval $[a, ae^S]$.

Let us examine the case $1 - 4\lambda = 0$. For $\lambda = 1/4$, the solution in $R$ has the expression (13). It may be shown that such a function cannot satisfy the boundary condition unless altitudes $a$ and $b$ are identical. In such a case, the cone is reduced to a simple cap. Thus, this solution has to be discarded.

## 3.2 Function $\Theta(\theta)$

Eq. (10b) is widely known as the Legendre equation whose solutions are the $P_n^m$ Legendre functions. The degrees $n$ are related to the eigenvalues $\kappa = n(n + 1)$. Legendre functions can take various representations but it is a classical result that the differential eq. (10b) is a particular case of the hypergeometric equation (Robin 1959; Hobson 1965). We rely on previous works and adopt the same conventions (Haines 1985). Hence:

$$P_n^m = \frac{(-1)^m}{2^m m!} N_n^m \sin^m \theta F \left( m - n, n + m + 1, 1 + m; \sin^2 \frac{\theta}{2} \right). \qquad (17)$$

The factor $N_n^m$ is chosen here to suit the Schmidt normalization:

$$\begin{cases} N_n^0 = 1 & \text{if } m = 0 \\ N_n^m = \frac{\Gamma(n + m + 1)}{\Gamma(n - m + 1)} & \text{if } m \neq 0. \end{cases}$$

It must be noted that this expression is valid only for $n \ge m$. This implies that functions involving $n < m$ are a particular case discussed below. We write explicitly the hypergeometric functions in order to describe each kind of Legendre functions. For any value of its arguments $m - n, n + m + 1, m + 1$ and $\sin^2 \frac{\theta}{2}$, the hypergeometric function can be expanded as:

$$F \left( m - n, n + m + 1, m + 1; \sin^2 \frac{\theta}{2} \right) = \sum_{j=0}^{\infty} \delta_j \sin^{2j} \frac{\theta}{2}, \qquad (18)$$

where $\delta_j$ is calculated by the recurrence relation:

$$\begin{cases} \delta_0 = 0 \\ \delta_{j+1} = \frac{(m - n + j)(n + m + 1 + j)}{(j + 1)(j + m + 1)} \delta_j. \end{cases} \qquad (19)$$

The general solution (19) holds for any complex values of $n$ and $m$ and hence will be relevant for solving problems 1 and 2. In any case, $m$ is a positive integer but noticeable differences occur as a result of the degrees $n$. Remembering the discussion in Section 3.1, in problem 1 two cases have to be analysed.

### 3.2.1 Case 1: $n = -1/2 + ip\pi/S$

Replacing $n$ in eq. (17) gives:

$$K_p^m = \frac{(-1)^m}{2^m m!} N_p^m \sin^m \theta F \left( m + \frac{1}{2} - i \frac{p\pi}{S}, m + \frac{1}{2} + i \frac{p\pi}{S}, m + 1; \sin^2 \frac{\theta}{2} \right). \qquad (20)$$

We changed the notation in order to follow previous works describing these functions (Zhurina & Karmazina 1966). They were first studied independently by Mehler and Heine in 1896, and are also known as the conical functions (Robin 1959). From the expansion of the
hypergeometric series, we have for the recurrence factor $\delta_j$:

$$ \left\{ \begin{array}{l} \delta_0 = 0 \\ \delta_{j+1} = \frac{(m + 1/2 + j)^2 + (p\pi/S)^2}{(j + 1)(j + m + 1)} \delta_j. \end{array} \right. \qquad (21) $$

From eq. (21), we conclude that the functions $F$ are null or positive and increasing functions of $\theta$. Depending on the aperture $\theta_0$ of the cone, they can reach very large values for large $m$ and $p$. This behaviour makes the numerical computation of the functions and their derivatives rather difficult. Note that, contrary to the usual Legendre functions appearing in geomagnetism, the conical functions do not oscillate around zero and are perfectly monotonic. The absence of roots is a consequence of the lack of a boundary condition for eq. (10b).

In the complex case, the gamma ratio can be expanded and simplified. We obtain (Zhurina & Karmazina 1966):

$$ N_p^0 = 1, $$

$$ N_p^m = \frac{\Gamma(1/2 + m + ip\pi/S)}{\Gamma(1/2 - m + ip\pi/S)} = (-1)^m \prod_{j=1}^{m} \left[ \left( \frac{p\pi}{S} \right)^2 + \left( \frac{2j-1}{2} \right)^2 \right]. $$

### 3.2.2 Case 2: $n = 0$ and $n = -1$

These functions are a particular case of the Legendre functions. It should be noted that $P_0^m = P_{-1}^m$ as a special case of the general relationship $P_n^m = P_{-n-1}^m$ (Robin 1959). However, eq. (17) or equivalently:

$$ P_n^m = \frac{(-1)^m}{m!} N_n^m \tan^m \frac{\theta}{2} F \left( -n, n+1, 1+m; \sin^2 \frac{\theta}{2} \right) \qquad (22) $$

no longer holds for $n, m$ integers and $m > n$ as a result of the singularity of $\Gamma(n - m + 1)$ in the expression of $N_n^m$. This case occurs here when $m > 0$. For such values, $P_n^m$, with $n, m$ integers and $-m \le n \le m + 1$, may be defined from eq. (22) with $n = n_0 + \varepsilon$, as the limit $\varepsilon \to 0$ of the function $\frac{P_{n_0+\varepsilon}^m}{\varepsilon}$. Actually, this definition leads to the property $P_n^m = -P_{-n-1}^m$ instead of $P_n^m = P_{-n-1}^m$. Therefore, we define $P_n^m$ in the case $n, m$ integer, $-m \le n \le m - 1$ and $m \ge 1$, by:

$$ P_n^m = \lim_{\varepsilon \to 0} \frac{P_{n_0-(-1)^m \varepsilon}^m}{\varepsilon}. \qquad (23) $$

It may be shown, knowing the behaviour of $\Gamma(z)$ around its poles, namely $\Gamma(-n + u) \sim \frac{(-1)^n}{un!}$, $n$ being a positive integer, that $P_n^m$ defined by eq. (23) takes the form:

$$ P_n^m = \frac{\Gamma(m+n+1)\Gamma(m-n)}{\Gamma(m+1)} \tan^m \frac{\theta}{2} \left[ 1 + \Gamma(m+1) \sum_{p=1}^{N} \frac{(-1)^p \Gamma(N+p+1)}{\Gamma(N-p+1)\Gamma(m+p+1)\Gamma(p+1)} \sin^{2p} \frac{\theta}{2} \right], $$

where $N = n$ if $n \ge 0$ and $N = -n - 1$ if $n \le -1$.

Thus, for $n = 0$ and equivalently for $n = -1$ we obtain the expression:

$$ P_0^m = P_{-1}^m = \Gamma(m) \tan^m \frac{\theta}{2}. \qquad (24) $$

Thereby, Legendre functions corresponding to the eigenvalue $\lambda = 0$ are in fact trigonometric functions. Eq. (24) holds true in the case $m = 0$, which gives the result $P_0^0 = 1$.

## 3.3 Expression of potential $V_1$

The solution of eq. (10c) is trivial:

$$ \Phi_p^m = G_p^m \cos(m\varphi) + H_p^m \sin(m\varphi) $$

with $\{G_p^m, H_p^m\}$ the Gauss coefficients associated with $\kappa < -1/4$ and

$$ \Phi_0^m = G_0^m \cos(m\varphi) + H_0^m \sin(m\varphi) $$

with $\{G_0^m, H_0^m\}$ the Gauss coefficients associated with $\kappa = 0$.

We finally recover the total potential $V_1$ solving the BV problem 1:

$$ \begin{aligned} V_1(r, \theta, \varphi) = & a \sum_{p \ge 1} \sum_{m > 0} R_p(r) [G_p^m \cos(m\varphi) + H_p^m \sin(m\varphi)] K_p^m(\theta) \\ & + a \sum_{m > 0} R_0 [G_0^m \cos(m\varphi) + H_0^m \sin(m\varphi)] P_0^m(\theta). \end{aligned} \qquad (25) $$
The first term corresponds to the eigenfunctions $R_p(r)$ obtained from the eigenvalues $\kappa < -1/4$, the second to the constant function in $r$ obtained from the eigenvalue $\kappa = 0$.

## 4 SOLUTION FOR THE POTENTIAL $V_2$

The solution of problem 2 has been given by Haines (1985). The boundary condition (8b) leads to a singular Sturm-Liouville problem where the Legendre eq. (10b) is associated to the particular boundary conditions $\Theta(0)$ finite and $\frac{d\Theta(\theta)}{d\theta}|_{\theta_0} = 0$. The first constraint eliminates the Legendre functions of the second kind, whereas the second one is satisfied for a discrete, infinite set of degrees $n$, real and larger than $m$. If we index the degrees $m$ with an integer $k$, they correspond to the $k-m$ even roots in Haines's notation (1985). In order to follow closely Haine's notation, we should write $k = m + 2k'(k' \in \mathbb{N})$. As the second Haines's solution $[k-m]$ odd, corresponding to the boundary condition $P_{n_k}^m(\theta_0) = 0$ is not incorporated into our set of basis functions, we index the roots of $dP_{n_k}^m/d\theta|_{\theta_0} = 0$ merely with ordinary indices belonging to $\mathbb{N}$, independently for $m$. However, it is obvious that the roots $n_k$ are actually functions of $m$ and $k$.

The function $R(r)$ reduces to eq. (12) with $A_1$ and $A_2$ real. The standard definition (17) of the Legendre functions is relevant in this case and the solution $\Phi$ is unchanged. We thus obtain for the potential field taking into account the finite upper spherical boundary of the domain:

$$ V_2(r, \theta, \varphi) = a \sum_{k \ge 0} \sum_{m>0} \left(\frac{a}{r}\right)^{n_k+1} \left[ G_{n_k}^{i,m} \cos(m\varphi) + H_{n_k}^{i,m} \sin(m\varphi) \right] P_{n_k}^m \\ +a \sum_{k \ge 0} \sum_{m>0} \left(\frac{r}{a}\right)^{n_k} \left[ G_{n_k}^{e,m} \cos(m\varphi) + H_{n_k}^{e,m} \sin(m\varphi) \right] P_{n_k}^m. \qquad (26) $$

$\{G_{n_k}^{i,m}, H_{n_k}^{i,m}\}, \{G_{n_k}^{e,m}, H_{n_k}^{e,m}\}$ are the local Gauss coefficients for the potential $V_2$. Let us recall that, apart from the index range adopted here, this expression is identical to the $k-m$ even expansion given by Haines (1985).

The superscripts $i$ and $e$ denote respectively what we usually call the 'internal' and 'external' coefficients. Nevertheless, let us note the quotation markers around internal and external. In the conical geometry, indeed, the terminology is a bit of a misnomer and the significance is quite different from the usual acceptance in spherical geometry (see Section 5.4).

To summarize, the local potential $V_c$ solving the initial problem, expanded in the local basis, is:

$$ V_c = V_1 + V_2. $$

$V_1$ is expanded in a class of harmonic functions on $\Omega$, $R_p K_p^m \cos m\varphi$ for $p > 0$ and $R_0 P_0^m \cos m\varphi$ for $p = 0$ (respectively for $\sin m\varphi$), verifying boundary conditions $\frac{\partial f}{\partial \nu} = 0$ on $\partial_b \Omega$. Hence, $V_1$ verifies the same conditions and because $\int_{\partial\theta_0} \frac{\partial V_1}{\partial \nu} ds = 0$ then $\int_{\partial\theta_0} \frac{\partial V_1}{\partial \nu} ds = 0$. This last property may be easily checked directly:

$$ \int_{\partial\theta_0\Omega} \frac{\partial V_1}{\partial \nu} ds = a \sum_{m=0}^{\infty} \sum_{p=1}^{\infty} \int_{\partial\theta_0\Omega} R_p(r) (G_p^m \cos m\varphi + H_p^m \sin m\varphi) \left. \frac{dK_p^m}{d\theta} \right|_{\theta_0} ds \\ +a \sum_{m=0}^{\infty} \int_{\partial\theta_0\Omega} R_0(r) (G_0^m \cos m\varphi + H_0^m \sin m\varphi) \left. \frac{dP_0^m}{d\theta} \right|_{\theta_0} ds \\ = 2\pi a \sum_{p=1}^{\infty} G_p^0 \left. \frac{dK_p^m}{d\theta} \right|_{\theta_0} \int_a^b R_p(r) dr. $$

Because the functions $R_p$ are orthogonal on $L^2(a, b)$, and because the constant function belongs to the base, the integral on the last line, which is the expression of the inner product $\langle R_p(r), 1 \rangle$ on $L^2(a, b)$, is null for every index $p$. Hence the property.

Similarly, $V_2$ is expanded in a class of harmonic functions on $\Omega$, $(\frac{a}{r})^{n_k+1} P_{n_k}^m \cos m\varphi$ and $(\frac{r}{a})^{n_k} P_{n_k}^m \cos m\varphi$ (respectively for $\sin m\varphi$) verifying the boundary condition $\frac{\partial f}{\partial \nu} = 0$ on $\partial\theta_0\Omega$. Therefore $V_2$ has the same property and $\int_{\partial\theta_0\Omega} \frac{\partial V_2}{\partial \nu} ds + \int_{\partial\theta_0\Omega} \frac{\partial V_2}{\partial \nu} ds = 0$.

The direct proof of this property is as straightforward as the previous one. It relies upon the fact that $P_{n_k}^0$ and $1$ are orthogonal with respect to the inner product:

$$ \langle f, g \rangle = \int_{0}^{\theta_0} f.g \sin \theta d\theta. $$

Thus, the total potential $V_c = V_1 + V_2$ belongs to a class of harmonic potentials on $\Omega$ that possess particular flux property, namely:

$$ \int_{\partial\theta_0\Omega} \frac{\partial V_c}{\partial \nu} ds = 0. \qquad (27) $$

As a result of the total null flux property of harmonic potentials, eq. (27) is equivalent to:

$$ \int_{\partial_a \Omega} \frac{\partial V_c}{\partial \nu} ds + \int_{\partial_b \Omega} \frac{\partial V_c}{\partial \nu} ds = 0. $$
5 RESOLUTION OF A BOUNDARY VALUE PROBLEM: RELATIONSHIPS WITH
GLOBAL GAUSS COEFFICIENTS

In this section, we will solve BV problems 1 and 2 where the boundary functions F, G₁ and G₂ are derived from a global model developed in ordinary SH. This BV problem will provide linear relationships between local and global Gauss coefficients, which are a generalization of the relationships established by De Santis et al. (1999) where the dependence in the radial variate was not taken into account. Such relationships will prove useful in view of the inverse problem, which consists of calculating a model based upon a set of field values given inside Ω.

A global field is written in the geocentric coordinate system (r, θ', φ'), the reference axis of which is the rotation axis of the Earth, while the local field is expanded in the local coordinate system (r, θ, φ), the reference axis of which is the axis of the cone. Therefore, we must first express the global coefficients {$g_n^m, h_n^m$} in a spherical coordinate frame the axis of which coincides with the axis of the cone. We know that each spherical harmonic defined in a system (θ', φ') can be represented by a linear combination of all the SH in the system (θ, φ) of same degree n (De Santis et al. 1996):

$$
P_n^m(\theta) \begin{Bmatrix} \cos(m\varphi) \\ \sin(m\varphi) \end{Bmatrix} = \sum_{\mu=0}^{n} \left[ \begin{Bmatrix} a_{n,\mu}^m \\ c_{n,\mu}^m \end{Bmatrix} \cos(m\varphi') + \begin{Bmatrix} b_{n,\mu}^m \\ d_{n,\mu}^m \end{Bmatrix} \sin(m\varphi') \right] P_n^\mu(\theta').
$$

As a result, the set of global rotated Gauss coefficients {$\tilde{g}_n^m, \tilde{h}_n^m$} related to the cone reference frame can be found from the set of Gauss coefficients {$g_n^m, h_n^m$} related to the geocentric reference frame by the linear combinations:

$$
\tilde{g}_n^m = \sum_{\mu=0}^{n} (a_{n,m}^\mu g_n^\mu + c_{n,m}^\mu h_n^\mu), \qquad (28)
$$

$$
\tilde{h}_n^m = \sum_{\mu=0}^{n} (b_{n,m}^\mu g_n^\mu + d_{n,m}^\mu h_n^\mu). \qquad (29)
$$

With these identities, a global field can be expanded in the local reference system (r, θ, φ) and its components can be used as boundary conditions G, F₁ and F₂. Before tackling the boundary problem, we have to return to the flux conditions.

5.1 Discussion on fluxes

We already mentioned that when decomposing the initial problem into two subproblems we cannot choose any boundary conditions F(r, φ), G₁(θ, φ), G₂(θ, φ). Let us assume that F, G₁ and G₂ derive from a global model, having its sources inside the sphere of radius r = a.

Expanded in the cone reference frame, its potential is expressed as usual by:

$$
\tilde{V}_{SH} = a \sum_{n=1}^{\infty} \sum_{m=0}^{n} \left(\frac{a}{r}\right)^{n+1} \left[ \tilde{g}_{n}^{m} \cos(m\varphi) + \tilde{h}_{n}^{m} \sin(m\varphi) \right] P_{n}^{m}(\theta), \quad (30)
$$

where $P_n^m$ are the Schmidt's functions with positive integer degrees $n$ and integer or null orders $m$ and $\{\tilde{g}_n^m, \tilde{h}_n^m\}$ the global Gauss coefficients with respect to the cone reference frame. Therefore:

$$
F(r, \varphi) = \tilde{X}_{SH} = \left. \frac{\partial \tilde{V}_{SH}}{r \partial \theta} \right|_{\theta_0}, \qquad (31)
$$

$$
G_1(\theta, \varphi) = -\tilde{Z}_{SH}(r=a) = -\left. \frac{\partial \tilde{V}_{SH}}{\partial r} \right|_a, \qquad (32)
$$

$$
G_2(\theta, \varphi) = \tilde{Z}_{SH}(r=b) = \left. \frac{\partial \tilde{V}_{SH}}{\partial r} \right|_b . \qquad (33)
$$

The flux condition derived from eq. (6a) now writes:

$$
\int \int_{\Omega_{\theta}} \tilde{X}_{SH} d\sigma_{\theta} = - \sum_{n=1}^{\infty} \sum_{m=0}^{n} \int_{a}^{b} \left(\frac{a}{r}\right)^{n+2} r dr \int_{0}^{2\pi} \left[\tilde{g}_{n}^{m} \cos(m\varphi) + \tilde{h}_{n}^{m} \sin(m\varphi)\right] \left.\frac{dP_{n}^{m}(\theta)}{d\theta}\right|_{\theta_{0}} \sin\theta_{0} d\varphi = 0.
$$

One can easily demonstrate that this expression reduces to:

$$
\int \int_{\Omega_0} \tilde{X}_{SH} d\sigma_\theta = 2\pi a \sin \theta_0 \sum_{n=1}^{\infty} \frac{\tilde{g}_n^0 [a - b(a/b)^{n+2}]}{1+n} \left. \frac{dP_n^0(\theta)}{d\theta} \right|_{\theta_0},
$$
so, writing $b = ae^S$ with $S \neq 0$, we conclude that the flux $\Phi$ is:

$$ \Phi = 2\pi a^2 \sin \theta_0 \sum_{n=1}^{\infty} \frac{\tilde{g}_n^0[1 - e^{-S(n+1)}]}{n+1} \frac{dP_n^0(\theta_0)}{d\theta}. \qquad (34) $$

For an arbitrary set of zonal harmonics, the field does not fulfil the divergence-free condition of problem 1. As such a field is physically divergence-free, there must be the same amount of flux into the spherical ends, so similarly it does not satisfy the condition for problem 2. Thus, this global field cannot be used as boundary conditions without modifications.

In order to circumvent this inconsistency it is possible to modify the complete BV problem in various ways. One possibility is to solve the problem for a potential $\tilde{V} = \tilde{V} - U_n$, where $U_n$ is the potential of a centred zonal multipole of degree $n$, the axes of which coincide with the cone axis. Its expression is $U_n = aC(\frac{a}{r})^{n+1}P_n^0(\theta)$, where the constant C has to be adjusted to make the integral of $\hat{F} = -\hat{X}_{SH} = -\frac{1}{r}[\frac{\partial(\tilde{V} - U_n)}{\partial\theta}]_{\theta_0}$ null on the boundary $\partial\Omega_{\theta_0}$, according to the condition (7). It comes out that:

$$ C = \frac{n\Phi}{2\pi a^2(1 - e^{-Sn})\sin\theta_0 \frac{dP_n^0(\theta_0)}{d\theta}}. \qquad (35) $$

If the functions $G_1$ and $G_2$ derive from the same modified potential $\tilde{V} - U_n$, the condition (9) will be obviously met and therefore so will that of eq. (7). Let us note that taking into account the modified potential requires merely the replacement of $\tilde{g}_n^0$ with $\tilde{g}_n^0 - C$ in eq. (30); all the other Gauss coefficients remaining unchanged. A more radical solution would be to suppress the zonal terms in the complete problem. Clearly, if there are no $m = 0$ terms, all flux conditions are met. Let us draw attention however to the fact that the removal of the zonal terms in the cone reference frame does not mean that we remove them in the geocentric reference frame of the Earth. In the following, we actually adopt the solution of removing zonal terms. For the sake of clarity, in the subsequent section, we use the symbol $\hat{\Phi}$ to denote the use of partial (specifically non-zonal) potentials and fields.

The latter discussion drives our attention to the important fact that when applying boundary conditions we do assume that they also solve the flux conditions. In this case, boundary conditions can not be applied arbitrarily and we will consider the partial potential field $\hat{F}$, $\hat{G}$ for theoretical discussion purposes.

## 5.2 Relationships from problem 1

The boundary conditions of problem 1 are now defined by eq. (31). Let us assume that the magnetic field on the boundary $\partial\Omega_{\theta_0}$ is equal to a field modelled by the SH for orders $m \ge 1$:

$$ X_c(r, \theta_0, \varphi) = \frac{1}{r} \frac{\partial V_c}{\partial \theta} \bigg|_{\Omega_\varphi} = \hat{X}_{SH,\theta_0}(r, \theta_0, \varphi). $$

Writing explicitly the expansion yields:

$$ \hat{X}_{SH,\theta_0} = \frac{a}{r} \sum_{n=1}^{\infty} \sum_{m=1}^{\infty} \left(\frac{a}{r}\right)^{n+1} \left[ \hat{g}_n^m \cos(m\varphi) + \hat{h}_n^m \sin(m\varphi) \right] \frac{dP_n^m(\theta)}{d\theta} \bigg|_{\theta_0}. $$

However, from eq. (8b), $X_{c,2}(r, \theta_0, \varphi) = 0$. Thus, $\hat{X}_{SH,\theta_0}(r, \theta_0, \varphi) = X_{c,1}(r, \theta_0, \varphi)$ with the $X_{c,1}$ component expanded in terms of conical functions at the colatitude $\theta = \theta_0$. The identification of the cosine and sine terms gives:

$$ \sum_{p=1}^{\infty} R_p(r) \left. \frac{dK_p^m(\theta)}{d\theta} \right|_{\theta_0} \left\{ \begin{matrix} G_p^m \\ H_p^m \end{matrix} \right\} + R_0 \left. \frac{dP_0^m(\theta)}{d\theta} \right|_{\theta_0} \left\{ \begin{matrix} G_0^m \\ H_0^m \end{matrix} \right\} = \sum_{n=m}^{\infty} \left(\frac{a}{r}\right)^{n+1} \left. \frac{dP_n^m(\theta)}{d\theta} \right|_{\theta_0} \left\{ \begin{matrix} \hat{g}_n^m \\ \hat{h}_n^m \end{matrix} \right\}. $$

Recall that the basis functions $\{R_0, R_p(r)\}$ are orthogonal. Hence, taking the inner product of the equation with respectively $R_p(r)$ and $R_0$ leads to the expressions of the Gauss coefficients:

$$ \left\{ \begin{matrix} G_p^m \\ H_p^m \end{matrix} \right\} = \sum_{n=m}^{\infty} \frac{\left\langle \left(\frac{a}{r}\right)^{n+1}, R_p(r) \right\rangle \frac{dP_n^m(\theta)}{d\theta} \bigg|_{\theta_0}}{\|R_p\|^2 \frac{dK_p^m(\theta)}{d\theta} \bigg|_{\theta_0}} \left\{ \begin{matrix} \hat{g}_n^m \\ \hat{h}_n^m \end{matrix} \right\} \qquad (36) $$

and

$$ \left\{ \begin{matrix} G_0^m \\ H_0^m \end{matrix} \right\} = \sum_{n=m}^{\infty} \frac{\left\langle \left(\frac{a}{r}\right)^{n+1}, R_0 \right\rangle \frac{dP_n^m(\theta)}{d\theta} \bigg|_{\theta_0}}{\|R_0\|^2 \frac{dP_0^m(\theta)}{d\theta} \bigg|_{\theta_0}} \left\{ \begin{matrix} \hat{g}_n^m \\ \hat{h}_n^m \end{matrix} \right\}. \qquad (37) $$
The inner products involved in these equations are:

$$ \left\langle \left(\frac{a}{r}\right)^{n+1}, R_p(r) \right\rangle = \frac{4\pi p [1 - (-1)^p e^{-S(n+1/2)}]}{\sqrt{(2S)(4\pi^2 p^2 + S^2)}} \frac{(n+1)}{\left[ (n+1/2)^2 + \left(\frac{\pi p}{S}\right)^2 \right]}, \\ \left\langle \left(\frac{a}{r}\right)^{n+1}, R_0 \right\rangle = \frac{1 - e^{-nS}}{n\sqrt{(e^S - 1)}}. $$

For a given couple $(m, p)$, $\{G_p^m, H_p^m\}$ and $\{G_0^m, H_0^m\}$ may be computed from eqs (36) and (37) respectively as long as $\frac{dK_p^m(\theta)}{d\theta}|_{\theta_0 \neq 0}$ in eq. (36), which is always true because $p \ge 1$, and $\frac{dP_0^m(\theta)}{d\theta}|_{\theta_0 \neq 0}$ in eq. (37), which would be not true for $m = 0$ (equivalent to including zonal terms in BV). In this latter case, this restriction would imply that the term $G_0^0$ is undefined and the indetermination reflects the fact that the potential is known within an additive constant.

To conclude this section, we may note that the relations (36) and (37) are established only for non-zonal terms. They may be extended to the zonal terms however, provided that the zonal Gauss coefficients $\hat{g}_n^m$ are modified.

## 5.3 Relationships from problem 2

Assume now that the modified boundary functions $\hat{G}_1$ and $\hat{G}_2$ in eqs (8c) and (8d) are equal to the Z component of a global internal field defined in Section 4.1 by eqs (32) and (33). The $\hat{Z}$ component of the global rotated field with order $m \ge 1$ is given by:

$$ \hat{Z}_{SH} = - \sum_{n=1}^{\infty} \sum_{m=1}^{\infty} (n+1) \left(\frac{a}{r}\right)^{n+2} \left[ \hat{g}_n^m \cos(m\varphi) + \hat{h}_n^m \sin(m\varphi) \right] P_n^m. \quad (38) $$

From eqs (6c) and (6d) $Z_{c,1}|_{a,b} = 0$, hence $Z_{c,1}|_{a,b} = Z_{c,2}|_{a,b}$, where $Z_{c,2}$ derives from eq. (26). When identifying the sine and cosine factors on both boundaries $\partial\Omega_a$ and $\partial\Omega_b$, and using the orthogonal properties of $P_{n_k}^m$, we obtain the following equations.

For $r = a$, that is to say on the boundary $\partial\Omega_a$, we have:

$$ \| P_{n_k}^m \|^2 (n_k + 1) \left\{ \begin{matrix} G_{n_k}^{i,m} \\ H_{n_k}^{i,m} \end{matrix} \right\} - n_k \| P_{n_k}^m \|^2 \left\{ \begin{matrix} G_{n_k}^{e,m} \\ H_{n_k}^{e,m} \end{matrix} \right\} = \sum_{n=m}^{\infty} (n+1) \langle P_n^m, P_{n_k}^m \rangle \left\{ \begin{matrix} \hat{g}_n^m \\ \hat{h}_n^m \end{matrix} \right\} $$

and for $r = b$ (on the boundary $\partial\Omega_b$) we get:

$$ \| P_{n_k}^m \|^2 (n_k + 1) \left(\frac{a}{b}\right)^{n_k + 2} \left\{ \begin{matrix} G_{n_k}^{i,m} \\ H_{n_k}^{i,m} \end{matrix} \right\} - n_k \| P_{n_k}^m \|^2 \left(\frac{b}{a}\right)^{n_k - 1} \left\{ \begin{matrix} G_{n_k}^{e,m} \\ H_{n_k}^{e,m} \end{matrix} \right\} \\ = \sum_{n=m}^{\infty} (n+1) \left(\frac{a}{b}\right)^{n+2} \langle P_n^m, P_{n_k}^m \rangle \left\{ \begin{matrix} \hat{g}_n^m \\ \hat{h}_n^m \end{matrix} \right\}. $$

Solving these latter equations for the Gauss coefficients $(G_{n_k}^{i,m}, H_{n_k}^{i,m})$ and $(G_{n_k}^{e,m}, H_{n_k}^{e,m})$, we find, for a couple $(m, n_k)$, the following relationships:

$$ \left\{ \begin{matrix} G_{n_k}^{i,m} \\ H_{n_k}^{i,m} \end{matrix} \right\} = \sum_{n=m}^{\infty} \frac{(n+1) \langle P_n^m, P_{n_k}^m \rangle (1 - e^{-S(n+n_k+1)})}{\| P_{n_k}^m \|^2 (n_k+1) (1 - e^{-S(2n_k+1)})} \left\{ \begin{matrix} \hat{g}_n^m \\ \hat{h}_n^m \end{matrix} \right\} \quad (39) $$

and

$$ \left\{ \begin{matrix} G_{n_k}^{e,m} \\ H_{n_k}^{e,m} \end{matrix} \right\} = \sum_{n=m}^{\infty} \frac{(n+1) \langle P_n^m, P_{n_k}^m \rangle (1 - e^{S(n_k-n)})}{\| P_{n_k}^m \|^2 n_k (e^{S(2n_k+1)} - 1)} \left\{ \begin{matrix} \hat{g}_n^m \\ \hat{h}_n^m \end{matrix} \right\}. \quad (40) $$

For $k = 0$ and $m = 0$ we have $n_k = 0$, so that the external coefficients $\{G_0^0, H_0^0\}$ are undefined in BV problems including modified zonal terms. As in the case of problem 1, this indetermination reflects the property of potential $V_2$ to be defined within an additive constant.

Thus, eqs (36), (37), (39) and (40) provide a linear relationship between the spherical and the conical expansions. They generalize the relations proposed by De Santis et al. (1999) and are independent of the radius $r$.

## 5.4 ‘Internal’ and ‘external’ contributions

We mentioned in Section 4 that the meaning of ‘internal’ and ‘external’ should not be confused with the classical interpretation of intern and extern, which refers to magnetic field sources respectively inside or outside the sphere of radius $a$. In regional modelling, ‘internal’ and ‘external’ coefficients are involved merely for geometrical reasons (see Section 7.1). It is straightforwardly seen that even in the case of a global internal magnetic field only, ‘external’ Gauss coefficients in eq. (40) are definitely not zero in a finite cone; in an infinite cone ($b$ being infinite on Fig. 1), the ‘external’ coefficients would disappear, assuming, of course, that the field modelled would be source-free in this infinite
cone. Even if this crucial question seems to be bound to this present formalism, it seems reasonable to suggest that the same difficulty could be persistent in Haines's method or even in any regional modelling technique (see for example Haines 1990).

The question arising is how to separate the true internal and external contributions in regional modelling? Before discussing this aspect, we should slightly generalize relations (36), (37), (39) and (40) in order to incorporate formally the global external coefficients. The computation steps are very similar to those sketched in Sections 2 and 3 where we have simply to replace $n + 1$ by $-n$. Let us assume in this section that both internal $\{\hat{g}_n^{i,m}, \hat{h}_n^{i,m}\}$ and external $\{\hat{g}_n^{e,m}, \hat{h}_n^{e,m}\}$ global Gauss coefficients are present in the right hand side of eqs (36), (37), (39) and (40). The question of distinguishing internal from external can be answered to the extent that we are able to invert the relations (36), (37), (39) and (40). Such an inversion is possible, at least theoretically, and is guaranteed by a remarkable unicity property of harmonic functions (Kellogg 1953). This property may be stated as follow: if two functions, harmonic on a spherical shell $a < r < b$ are equal on the conical domain $\Omega = \{0 < \theta < \theta_0, a < r < b\}$, which is a subset of the shell, they are equal all over the shell. In other words, there is a unique harmonic function on the shell (global potential) equal to the harmonic function (local potential) on the conical domain. Thus, having found the local Gauss coefficients by some inversion of data inside the cone, it should be possible to calculate global coefficients and, hence, to evaluate internal and external contributions. In practice, however, the regional as well as the global Gauss coefficients have to be restricted to a finite set, which raises the problem of truncation level of the series. Moreover, the linear system relating the global coefficients to the local ones is badly conditioned, perhaps as badly as the linear problem consisting in computing directly the global coefficients from data restricted to the conical domain. Therefore investigating the characteristics of the ill conditioning as well as the relationship between the two types of problems will require further theoretical work.

# 6 NUMERICAL ASPECTS

Among the functions we went through, the generalized Legendre functions as well as the conical functions deserve special attention. Although modern computers brought a major breakthrough for the calculation of special functions, their evaluation must be handled with care especially in the case of a small cone.

## 6.1 Associated Legendre's function $P_n^m$

The associated Legendre functions $P_n^m$ are solution of a Sturm-Liouville problem and oscillate around zero. The recurrence factor $\delta_j$ in the hypergeometric representation (19) can be either positive or negative depending on the value of $n_k$. One of the main problems when using the direct series expansion is the bad cancellation of large opposite terms with approximately the same absolute value. The Clenshaw's recurrence formula is a well-known way for bypassing this difficulty (Press et al. 1992). However, even if it improves accuracy for relatively small degrees $n_k$, the computation is also defective when $n_k$ increases (Thébault et al. 2002). If this method, indeed, makes it possible to get rid of the sum of the terms, it gives place to the product of the terms, one very small, the other very large, causing the usual failure in floating point representation. This failure is especially encountered for large degrees $n_k$. Remember that the values of $n_k$ strongly depend on the half angle of the cap $\theta_0$; the smaller the cap the larger the degrees are. Note that, numerically, the terms, small or large, should be considered more in relation to the magnitude of the degrees $n_k$, than simply relative to the geometry and the aperture $\theta_0$ of a spherical cap. For all these reasons, we recommend Olver and Smith's algorithm (Olver & Smith 1983) that uses the extended arithmetic range concept and produces very accurate results compared with any other method. Thorough investigations were undertaken and as far as we tested the algorithm, it gives at least a precision of 15 digits, reducing to 12 for very large $n_k$. The limit of the algorithm is nevertheless reached for large value of $n_k$ and $m \sim n_k$ but, in any case, our interval of application is inside the confidence interval of the program. For instance, the computation of Legendre functions is still stable and accurate for $k = 500, m = 50, \theta_0 = 5^\circ$ and $n_{500(50)} \sim 9008.5$. Compared with the Clenshaw's recurrence method, the accuracy of the Legendre's function computation is highly increased for small caps. It also ensures us to find accurately the roots $n_k$ using eq. (8b), which is essential considering the behaviour of the radial functions $(r/a)^{n_k}$ when coping with altitude.

## 6.2 Mehler or conical functions

On the other hand, the Mehler functions can be expanded with the hypergeometric series using eq. (21). The $\delta_j$ factor is strictly positive, thus no bad cancellation arise. The Clenshaw's recurrence formula is a robust process for the evaluation of these functions and our algorithm is based on this method.

Usually, when dealing with special functions, we can rely on previous numerical work. Unfortunately, these functions are a special case and there is no available algorithm to compute them in the general case (Segura, private communication). A computer program is available only for $m = 0$ and $m = 1$ (Köblig 1981) but it is based on the Mehler-Fock integral representation of the conical functions (Zhurina & Karmazina 1966),

$$K_p^m = \frac{2}{\sqrt{\pi}} \frac{\tan^m(\theta/2)}{\Gamma(m + 1/2)} N_p^m \int_0^{2\pi} \frac{\cosh\{2\pi p[\sin(\theta/2) \sin(\beta)/S] \cos^{2m}(\beta)\}}{[1 - \sin^2(\theta/2) \sin^2(\beta)]^{1/2}} d\beta. \quad (41)$$

With this formalism, the conical functions at any order $m$ must be computed using a recurrence formula involving functions with same degree $n_p$ but different orders $m$. Unfortunately, a cot $\theta$ in the recurrence relation produces a singularity generating propagating inaccuracies
when $\theta$ is small. Thus, the recurrence relation is practically of no use. Even if we preferred the hypergeometric expansion using the Clenshaw recurrence to the integral representation, eq. (41) represents, nevertheless, an interesting alternative way of computation of the conical functions. In particular, we found it useful to check our algorithm. In order to carry out these checks, we computed eq. (41) directly using a numerical integration based on a Gauss quadrature. We then compared results given by eq. (41) with the hypergeometric representation (20) computed with the Clenshaw recurrence. The relative error was less than $10^{-8}$ for $m \in [0, 20]$, $p \in [1, 30]$ and $\theta \in [0, 25]$, showing the efficiency of the evaluation. However, this check is limited because, the direct use of the integral representation (41) is time consuming and not accurate when $p/S$ and $\theta$ are large. Indeed, as the function under the integral sign is sharply concentrated in a few peaks, it is therefore unsuitable for numerical integration. Hence, we actually restricted ourselves to the interval $\theta \in [0, 25]$. Although usual recurrence relations for Legendre functions can be used to further check the results for $\theta \ge 25^\circ$, beyond this value even the Clenshaw method must be used wisely and it is likely that a deeper work should be done to adapt the extensive arithmetic process used by Olver and Smith to the computation of the conical functions.

# 7 RECONSTRUCTION OF A GLOBAL FIELD

The reliability of the formalism described above was tested by comparing an IGRF model, rotated first in the cap reference frame and then expanded in SH for $m \ge 1$, with a reconstructed model over a cone. We started arbitrarily from a DGRF90 global model (of course, any global model would have been suitable). The lower cap of the cone was chosen in order to encompass France whose centre of coordinates is $C(\varphi' = 2^\circ E, \theta' = 47^\circ N)$ and the aperture was set to $\theta_0 = 5^\circ$.

Before applying eqs (36), (37), (39) and (40) to obtain the local coefficients, we first rotated the global Gauss coefficients $\{G_n^m, h_n^m\}$ using eqs (28) and (29) in order to get rotated values $(\tilde{G}_n^m, \tilde{h}_n^m)$. The method was settled long ago by Bernard *et al.* (1969) but a numerical process based on least squares has been more recently proposed by De Santis *et al.* (1996) and has proved its efficiency.

Once the coefficients had been rotated, we removed the zonal contributions of the global field to satisfy the flux condition $\Phi = 0$ (see Section 5.1) for both independent problems 1 and 2; thus, $(\hat{G}_n^0, \hat{h}_n^0) = 0$, following the convention adopted in Section 5.1. The boundaries at $r = a$ and $r = b$ were chosen in order to include typical altitudes of satellite data. We set $a$ the terrestrial radius and $b = a + 1030$ km.

Then, eqs (36), (37), (39) and (40) were applied. In theory, the sums over $p$ and $k$ in the local expansions of $V_1$ (25) and $V_2$ (eq. 26) are infinite. For the numerical computations we truncated the sum to $p_{max} = k_{max} = 11$. This gives a total of 292 Gauss coefficients among which 20 are Gauss coefficients corresponding to $\{G_0^m, H_0^m\}$, 130 to $\{G_p^m, H_p^m\}$, and 130 to $\{G_{n_k}^{i,m}, H_{n_k}^{i,m}\}$. Because of their relatively small magnitude, most of the $\{G_{n_k}^{e,m}, H_{n_k}^{e,m}\}$ were at a first sight insignificant and we kept only 12 of them. We will see in Section 7.1 that removing some ‘external’ coefficients especially affects the model with increasing altitude. The components of the reconstructed field were obtained from $V_1$ and $V_2$. The final field is $\vec{B}_c = \vec{B}_{c,1} + \vec{B}_{c,2}$.

The method obviously only holds within the cone domain $\Omega$. Our purpose being to use the method over France, we especially focused on the altitudes 0, 4, 400 and 800 km, where we have respectively at our disposal: repeat station data, aeromagnetic data, Magsat and Ørsted satellite data. We also sketched the contour charts at altitudes 20 and 100 km to support the discussion. Figs 2, 3 and 4 show the reconstructed field components $X_c, Y_c$ and $Z_c$ respectively in the cone and Figs 5, 6 and 7, the residuals between the local $X_c, Y_c, Z_c$ components and the $\tilde{X}_{SH}, \tilde{Y}_{SH}, \tilde{Z}_{SH}$ components from the rotated DGRF90 model computed over the cone using spherical harmonic expansion with $m \ge 1$. Note that the contour lines outside the caps (dotted line circles) are meaningless.

## 7.1 Discussion

Broadly, we see that for each component, the residuals are maximum at the surface of the Earth, which corresponds to the lower limit of the Sturm-Liouville problem 1. This area is prone to Gibbs phenomena as a result of the truncation at $p_{max}$ in the expansion of $r^{-n-1}$ of the basis functions $R_p(r)$. This means that, in order to use terrestrial data, we must use more parameters $\{G_p^m, H_p^m\}$ than actually needed at an upper altitude in order to avoid, as much as possible, the surface edge effects. Similarly, the circular anomalies come from the truncation at $k_{max}$ in the expansion of the basis functions $P_{n_k}^m$, because $\theta_0$ is the boundary of the Sturm-Liouville problem 2. In addition, the circle defined by $r = a$ and the colatitude $\theta_0$ is a geometric singularity where edge effects from the boundary $\theta_0$ and $r = a$ overlap, producing a maximum of discrepancies between the local and the global models. The same reasoning can be applied to the upper circle at $r = b$ and $\theta_0$.

Starting from $r = a$, the accuracy of the model improves with altitude until the middle of the cap interval. We can clearly observe that the residuals decrease until altitude $(b-a)/2$. The circular effects also reduce, proving that at the surface of the Earth, the edges effect are partially a result of the $R_p(r)$ functions in the corners. So we get what is theoretically expected: the modelling error is minimum in the vicinity of $\theta = 0$ and altitude $(b-a)/2 \approx 500$ km. The shape of the isolines on the $X, Y, Z$ components at 400 km, close to the optimal altitude, reflects this behaviour.

Quantitatively, the residuals over France never exceed 5 nT for the $X$ components even at the lower limit at $r = a$, which is, as mentioned before, a critical area. The main discrepancies lie along the dotted line circles. For the altitudes shown, the smallest residuals within the circle seen on Fig. 5 are at 400 km and are less than 0.2 nT. In a similar way the $Y$ residuals are less than 5 nT and have a minimum of 0.2 nT at 400 km (Fig. 6). The residuals in $Z$ show a maximum of 2 nT at altitude 0 km (Fig. 7). In any case, when comparing the values of residuals with the values of the field (Figs 5, 6 and 7) we can consider that the field has been successfully reconstructed.

Special attention is devoted to the altitude 800 km. We distinguish organized features for the contour lines of each component. In fact, for discussion purposes, we have intentionally removed the coefficients $\{G_{n_k}^{e,m}, H_{n_k}^{e,m}\}$, which were at first sight insignificant. However, the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Component $X_c$ of the modelled magnetic field for orders $m > 0$. Several altitudes are represented. The dotted circle represents the edges of the spherical cap. The component is sketched in the geographic reference frame.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Same legend as Fig. 2 for the $Y_c$ component of the modelled magnetic field.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same legend as Fig. 2 for the component $Z_c$ of the modelled magnetic field.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Residuals between the *X* component of the global DGRF90 model and the reconstructed model (see Fig. 2) computed for expansions up to $k_{max} = p_{max} = 11$. Contour intervals are 5, 2, 0.5, 0.1 and 0.1 nT for the altitudes: 0, 20, 100, 400 and 800 km, respectively. The residuals are biased when increasing altitude as a result of the arbitrary removal of ‘external’ coefficients.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Residuals between the *Y* component of the global DGRF90 model and the reconstructed model (see Fig. 3). Same legend as Fig. 5.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Residuals between the Z component of the global DGRF90 model and the reconstructed model (see Fig. 4) computed for expansions up to $k_{max} = p_{max} = 11$. Contour intervals are in this case 5, 5, 1, 1 and 1 nT for the altitudes 0, 20, 100, 400 and 800 km, respectively. See Fig. 5 concerning the increasing bias with altitude.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Variation with respect to altitude of the mean square residuals between true data, derived from a DGRF90 model, and SCHA modelling of X, Y and Z. See text for further discussion.</figcaption></figure>

‘external’ contribution is function of $(\frac{r}{a})^{n_k}$ and for small caps, the degrees $n_k$ reach large values. At high altitudes, even the external Gauss coefficients with low magnitude have a significant weight and consequently cannot be removed arbitrarily.

Fig. 8 displays the systematic investigation of the mean square residuals between the reconstructed field and the original field with the radial distance $r$, where the mean is taken within the cone. One hundred models, sampled every 10 km, were computed between the altitudes 0 and 1030 km and compared directly with the initial DGRF90 model for $m \ge 1$.

The most remarkable feature is the almost perfect symmetry of the magnitude of the residuals, especially at the upper and the lower surfaces, subjected to the highest Gibbs phenomenon. Actually, this behaviour, at first sight satisfying, reflects the fact that the reconstructed field is less accurate with increasing altitudes, because residuals are relative residuals. Adding more ‘external’ coefficients would give, indeed, a shape more consistent with the variation of internal potential fields (sources inside the Earth) with altitude.

It should be noted that the large values of mean square residuals reached on the lower and upper caps are the result of discrepancies concentrated along the edges ($\theta = \theta_0, r = a$) or ($\theta = \theta_0, r = b$). Removing a small (say half a degree) strip would drop the mean square residuals to approximately 5 nT proving that the reconstruction is very stable inside the cone. The oscillations and their amplitude indicate the accuracy in which the radial function $r^{-n}$, with $n$ integer, are modelled and are related to the truncature index $p$. Increasing the index $p$ will, of course, reduce the amplitude of the oscillations (as well as the mean square residuals).

The same kind of study was carried out with several sizes of cone between 5° and 25°. If we increase the upper altitude, we need more parameters from the basis of Mehler functions; conversely, increasing the half-angle of the cap require a larger expansion of the basis of Legendre functions.

# 8 FURTHER INQUIRIES

The BV problem settled in Section 1, although the most natural one, because it implies the field (measurable quantity) and not the potential, is not the unique one leading to the construction of bases like those described in Sections 2 and 3. Other BV problems may be imagined. We will show in a forth coming paper a more general family of BV problems and, as a result, a basis harmonic function that removes the problem of the particular flux conditions expression in eq. (27). However, the mathematics sketched in Sections 2 and 3 will remain essentially the same. In particular, whatever the BV problem, the solution may be expanded into two families of basis harmonic functions. One basis involves radial
functions $R_p$ and Mehler functions $K_p^m$, the second one generalized Legendre functions $P_{n_k}^m$. Before tackling the inverse problem, it would be useful to investigate orthogonality properties of the basis functions, which play a prominent role with respect to this aspect. Bearing in mind the common framework of inverse theory in geomagnetism, which involves the field rather than the potential, a relevant inner product would be:

$$ \langle u, v \rangle = \int_{\Omega} \overline{\nabla u} \cdot \overline{\nabla v} d\tau. \qquad (42) $$

It should be noted however, that eq. (42) does not possess all properties of an inner product on the general Sobolev space $H^1(\Omega)$, as $|| \mathbf{u} || = 0$ implies only $u = \text{constant}$. Eq. (42) defines truly an inner product only on specific subspaces of $H^1(\Omega)$. Leaving aside this detail, which is not really inconvenient in practice, it may be shown that the two bases, elaborated in Sections 3 and 4, are orthogonal within themselves but non-orthogonal between them. This inherent non-orthogonality might be an inconvenience in inverse problem dealing with uneven data distribution over a conical domain, which is unfortunately the common situation.

We will, in a forth coming paper, propose two decompositions leading to fully orthogonal functions. However, beside the orthogonality properties, the specific choice of a set of harmonic basis functions should also take into account the ability of the series expansion to reduce Gibbs phenomenon and, hence, to provide satisfying modelling all over the domain. This property is ultimately connected to the ability of the $R_p$ functions to model $r^n$ over the interval $[a, b]$, where $n$ is a negative or positive integer, and to the ability of the generalized Legendre function $P_{n_k}^m$ to represent ordinary Legendre functions $P_n^m$, where degrees $n$ are positive integers. Both orthogonality properties and Gibbs phenomenon reduction depend on the boundary conditions adopted in the Sturm-Liouville problems verified by $R_p$ and $P_{n_k}^m$ respectively, so that a kind of trade-off has to be found for the basis functions to fulfil both types of constraints.

## 9 CONCLUSION

We present here a new mathematical formalism for SCHA modelling, which models properly the radial variation of the field. Setting, in a first step, a BV problem on a conical domain $\Omega$, bounded by a conical surface $\partial\Omega_\theta(\theta = \theta_0)$ and two spherical surfaces $\partial\Omega_a(r = a)$ and $\partial\Omega_b(r = b)$, we show that this problem may be split into two subproblems noted 1 and 2 having complementary boundary conditions the combination of which reconstitutes the boundary conditions of the initial problem. Problem 2 leads to a set of basis functions identical to one proposed by Haines (1985). The solution of problem 1 is the very originality of this new formalism. It relies upon a Sturm-Liouville problem for the radial function $R(r)$, which, in turn, involves colatitude functions called conical or Mehler functions, which are generalized Legendre functions with degrees of the form $n = -1/2 + ip$. The orthogonal basis functions $R_p(r)$ are able to represent properly any functions $r^n$ or $r^{-n-1}$, $n$ integer, which is ultimately the correct variation of a Newtonian potential field with respect to the radial distance to the sources. The efficiency of the formalism has been demonstrated by the modelling of a global field, which remains accurate over the whole domain $\Omega$. Taking into account that a common practice in regional modelling includes the preliminary removal of a global model, this accuracy should increase when modelling residual fields. Solving the BV problem with this formalism provides linear relations between global and local Gauss coefficients, which, contrary to the relationships obtained from Haines's formalism, are independent of the radial distance $r$. These relations will provide sensible prior information on variance, which should be helpful in inverse problems involving sparse data. They are also in the heart of the separation of internal and external sources, in so far as this problem can be solved with regional modelling.

In this paper, the decomposition is based on Neumann conditions only. This choice seems straightforward in the sense that it involves the field, which is the measurable quantity, and not the potential. However, it allows only to model a restricted class of fields, verifying the condition expressed by eq. (27). This restriction is, of course, a major drawback for the application of this modelling to real data (although one might argue that true regional fields, created by lithospheric sources only, should nearly fulfil the flux constraint). The important point is that, even if it applies to a special class of fields, the BV problem approach adopted in this paper leads to the construction of two sets of functions, which build up an appropriate basis for harmonic potentials on the conical domain $\Omega$. Their mathematical formulation will remain essentially the same in the construction of basis functions relevant for the modelling of general Newtonian potentials.

## ACKNOWLEDGMENTS

This is the CNRS contribution no. 2004.15 UMR 7516. The authors are grateful to the referees A. De Santis and F. J. Lowes, whose comments were quite helpful in the clarification of the mathematics and the improvement of the manuscript.

## REFERENCES

Alldredge, L., 1981. Rectangular Harmonic Analysis applied to the geomagnetic field, *J. geophys. Res.*, **86**(B4), 3021–3026.

Bernard, J., Kosik, J.-C., Laval, G., Pellat, R. & Philippon, J.-P., 1969. Représentation optimale du potentiel géomagnétique dans le repère d'un dipole décentré, incliné, *Ann. Geophys.*, **25**(2), 659–665.

De Santis, A. & Falcone, C., 1995. Spherical Cap models of Laplacian potentials and general fields, in *Geodetic Theory Today*, pp. 141–150, ed. Sansò, F., Springer, Berlin.

De Santis, A., Kerridge, D. & Barraclough, D., 1989. A Spherical Cap Harmonic Analysis model of the crustal magnetic anomaly field in Europe observed by Magsat, *Geomagnetism and Palaeomagnetism*, **261**, 1–17.

De Santis, A., Torta, J. & Falcone, C., 1996. A simple approach to the transformation of spherical harmonic models under coordinate system rotation, *Geophys. J. Int.*, **126**, 263–270.

De Santis, A., Falcone, C. & Torta, J., 1997. SHA vs SCHA for modelling secular variation in a small region such as Italy, *J. Geomag. Geoelectr.*, **49**, 359–371.
De Santis, A., Torta, J. & Lowes, F., 1999. Spherical Cap harmonics revisited and their relationship to Ordinary Spherical Harmonics, *Phys. Chem. Earth*, **24**(11–12), 935–941.

Haines, G.V., 1985. Spherical Cap Harmonic Analysis, *J. geophys. Res.*, **90**, 2593–2598.

Haines, G.V., 1988, Computer Programs for Spherical Cap Harmonic Analysis of potential and general fields, *Comput. Geosci.*, **14**, 4, 413–447.

Haines, G.V., 1990. Regional magnetic field modelling: a review, *J. Geomag. Geoelectr.*, **42**, 1001–1018.

Hobson, E., 1965. *The Theory of Spherical and Ellipsoidal Harmonics*, 2nd reprint edn, Chelsea Publishing Co., New York.

Hwang, C. & Chen, S., 1997. Fully normalized Spherical Cap Harmonics: application to the analysis of sea-level data from TOPEX/POSEIDON and ERS-1, *Geophys. J. Int.*, **129**, 450–460.

Kellogg, O.D., 1953. *Foundations of Potential Theory*, Dover, New York.

Köblig, K., 1981. A program for computing the conical functions of the first kind $P^m_{1/2+it}(z)$ for $m = 0$ and $m = 1$, *Computer Physics Communication*, **23**, 51–61.

Korte, M. & Haak, V., 2000. Modelling european magnetic repeat station and survey data by SCHA in search of time-varying anomalies, *Phys. Earth planet. Int.*, **122**, 205–220.

Langel, R. & Hintz, W., 1998. *The magnetic field of the Earth's lithosphere*, Cambridge University Press, Cambridge, UK.

Olver, F. & Smith, J., 1983. Associated Legendre functions on the cut, *Journal of Computational Physics*, **51**, 502–518.

Press, W., Teulolsky, S., Vetterling, W. & Flannery, B., 1992. *Numerical Recipies in C*, 2nd edn, Cambridge University Press, New York.

Robin, L., 1959. Fonctions Sphériques de Legendre et Fonctions Sphéroidales, *Vol. I–III*, Gauthier-Villars, Paris.

Thébault, E., Schott, J. & Mandea, M., 2002. Geomagnetic field modelling on small spherical caps, Communication at EGS XXVII general assembly, EGS02-1-04673 (SE6.07-1WE2A003), Nice, France.

Zhirina, M. & Karmazina, L. N., 1966. *Tables and Formulae for the spherical functions* $P^m_{1/2+it}(z)$, Pergamon Press, New York.

Zwillinger, D., 1989. *Handbook of differential equations*, Academic Press, New York.