# Earth magnetic field modeling from Oersted and Champ data

J. L. Le Mouël, P. Shebalin, A. Khokhlov

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

## To cite this version:

J. L. Le Mouël, P. Shebalin, A. Khokhlov. Earth magnetic field modeling from Oersted and Champ data. Earth Planets and Space, 2010, 62, pp.277-286. 10.5047/eps.2009.11.003. insu-03605342

HAL Id: insu-03605342

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

Submitted on 11 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Earth magnetic field modeling from Oersted and Champ data

J. L. Le Mouël¹, P. Shebalin², and A. Khokhlov²

¹Institut de Physique du Globe de Paris, 4, Place Jussieu, 75252, Paris, France

²International Institute of Earthquake Prediction Theory and Mathematical Geophysics, 84/32, Profsoyuznaya str., 117997, Moscow, Russia

(Received July 11, 2008; Revised October 9, 2009; Accepted November 4, 2009; Online published March 4, 2010)

We present a method to model geomagnetic field requiring only a restricted number of measurements on magnetic survey satellite orbits. These points are chosen in an optimal—or close to optimal—manner relying on recent developments in the problem of numerical integration over spheres. The method allows us to compute a series of models at short time intervals, namely 10 days in the present study. At each of these close dates several models are computed from independent sets of data; their redundancy in turn provides a control of results thanks to which the selection of data—for example, as a function of magnetic activity or latitude—may be reduced. We find that the internal low degree Gauss coefficients derived from Oersted and Champ data, respectively, differ from one another by 1 or 2 nT. We then take as a second example of the method application a brief study of the so-called external field. We compare the first-degree axisymmetric field with the $D_{st}$ index.

**Key words:** Geomagnetic field model, satellite, numerical integration over spheres.

## 1. Introduction

The problem of computing a model of magnetic field that first fits ground observations, then satellite observations (since the years 1960) goes back to Gauss in the 1830s. The most recent models using Oersted and Champ data rely on a least squares technique providing the spherical harmonic coefficients, and some of these solve the Euler angles of the sensor attitude at the same time (Olsen et al., 2006). There is probably no need to abandon spherical harmonic expansion, which is so practical for all applications of the models. Nevertheless, even when keeping to this classical method, different options exist. We will describe and use one of these, the choice of which is guided by a few considerations.

First, it is desirable to retain as much data as possible from high (>55°) latitude regions, despite the large disturbances which are present in those regions. Second, the key issue in spherical harmonic expansion modeling, i.e. the computation of Gauss coefficients, both internal and external, is the geographical distribution of data; integral orthogonality properties of the harmonics over the data points must be strongly adhered to, not loosely. Third, we call for a flexible algorithm, using the minimum but sufficient number of data points to compute the internal and external fields, at short time intervals, for night hours or day hours, different universal times, and different conditions of activity, etc.

Excellent models have been computed by different teams according to their own methods: using data from Pogo, Magsat, Oersted and Champ satellites as well as from

ground magnetic observatories (Sabaka et al., 2002) and Sabaka et al. (2004) have constructed comprehensive models covering the time interval (1960–2002). Maps of the internal field (core + crustal) and the external field (ionospheric and magnetospheric) are made available with a short time sampling. Other teams have produced models of the different ingredients of the geomagnetic field from satellite data covering the time interval (2000, 2005) (Maus et al., 2006; Olsen et al., 2006).

We present the general features of our method in the main text reserving more mathematical considerations for Appendix. Its efficiency will be demonstrated in the following sections.

## 2. The Model

### 2.1 Optimal arrays of points

Let us consider a sphere $S$ of unit radius. An optimal array—in the sense to be defined below—of $N$ points on $S$ is made of the points whose colatitude $\theta_k$ and longitude $\varphi_k$, $k = 1, \dots, N$, are given by the formulae below

$$ \theta_k = \arccos(h_k), \quad h_k = -1 + \frac{2k-2}{N-1}, \quad 1 \le k \le N $$

$$ \varphi_k = \left( \varphi_{k-1} + \frac{\text{const}}{\sqrt{N}} \cdot \frac{1}{\sqrt{1-h_k^2}} \right) (\text{mod } 2\pi), \quad 1 < k < N; $$ $$ \varphi_1 = \varphi_N = 0 $$

A subroutine computes the longitudes $\varphi_k$ and colatitudes $\theta_k$ of the optimal array, $N$ being given (Appendix). In this paper, we will take $N = 1000, 2000, 3000, 5000, 10000$.

Let the surface harmonic functions

$$ Y_{[n,m,c]}(\theta, \varphi) = P_n^m(\cos \theta) \cos m\varphi $$

$$ Y_{[n,m,s]}(\theta, \varphi) = P_n^m(\cos \theta) \sin m\varphi $$

<footnote>Copyright © The Society of Geomagnetism and Earth, Planetary and Space Sciences (SGEPSS); The Seismological Society of Japan; The Volcanological Society of Japan; The Geodetic Society of Japan; The Japanese Society for Planetary Sciences; TERRAPUB.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Examples of optimal array {$Q_l$} (black points) and close to optimal {$P_l$} (red points) extracted from 130 days of data centered on at 2003.0 (a) $N = 1000$, (b) $N = 3000$.</figcaption></figure>

be ranked in the usual lexicographic order, and $u_j$ be the corresponding jth harmonic, internal or external (see Eq. (2)), in the corresponding series. The gradients $\nabla u_i$, $\nabla u_j$ are orthogonal on the set of $N$ points $Q_l, l = 1, \dots, N$. (see Appendix)

$$ (u_i, u_j) = \sum_l \nabla u_i(Q_l) \cdot \nabla u_j(Q_l) = C_i \delta_{ij} \quad (1) $$

with $C_i$ being the corresponding norm of $\nabla u_i$.

## 2.2 Array of points close to an optimal array

Let us compute, as an example, a model based on night values of a given period of time. We first retain all measurement points $P_k$ satisfying this local time condition. Let $E(P_k) = E(P_k(r, \theta, \varphi))$ be this set of points. We pick up in $E(P_k)$ the point $P_i$ which is the closest to the point $Q_i$ of a given optimal array, in the sense that the angle $P_i O Q_i$ ($O$ being the sphere center) is the smallest (Fig. 1). The number $N$ of points $P_i$, equal to the number of points $Q_i$ of the optimal array, is much smaller than the number of points $P_k$. After this operation is made for each point $Q_i$, we are left with an $N$ array close to the optimal: $E_N(P_i(r_i, \theta_i, \varphi_i))$. Now, let the geomagnetic potential to be computed be written in the usual form:

$$ \begin{aligned} V(P) &= V(r, \theta, \varphi) \\ &= \sum_{j=1}^{K^{[\text{int}]}} g_j u_j^{[\text{int}]}(r, \theta, \varphi) + \sum_{j=1}^{K^{[\text{ext}]}} \gamma_j u_j^{[\text{ext}]}(r, \theta, \varphi) \end{aligned} \quad (2) $$

Let $\vec{B}(P_l) = (X(P_l), Y(P_l), Z(P_l))$ be the vectorial measurement of $\vec{B}$ at $P_l \in E_N(P)$. We want to compute the Gauss coefficients $g$ and $\gamma$ such that

$$ \sum_{l=1}^{N} \left( \vec{B}(P_l) - \nabla V(P_l) \right)^2 \quad (3) $$

$V$ contains $K^{[\text{tot}]} = K^{[\text{int}]}(K^{[\text{int}]} + 2) + K^{[\text{ext}]}(K^{[\text{ext}]} + 2)$ unknown coefficients, and there are $3N$ equations (3). In this paper we will address the core field and, briefly, a component of the so-called external field, with $K^{[\text{int}]} = 165$, and $K^{[\text{ext}]} = 2$. For such values, $N = 1000$ already provides a largely overdetermined set of equations. The system
<figure><img src="image_4.png" /><figcaption>Fig. 2. Examples of the evaluation the Gauss coefficients estimation as a function of the number of points N: left coefficient $g_1^1$, right coefficient $h_3^2$. Error bars are shown.</figcaption></figure>

(3) will be solved by the usual least squares technique, and the inversion of the resulting normal matrix by the Singular Value Decomposition algorithm. When computing this matrix, the following approximations of scalar products in Eq. (1) appear:

$$ \langle u_i^{[\text{int}]}, u_j^{[\text{int}]} \rangle, \quad \langle u_i^{[\text{ext}]}, u_j^{[\text{ext}]} \rangle, \quad \langle u_i^{[\text{ext}]}, u_j^{[\text{int}]} \rangle $$

with, for example,

$$ \langle u_i^{[\text{int}]}, u_j^{[\text{int}]} \rangle = \sum_{l=1}^{N} \nabla u_i^{[\text{int}]}(P_l) \cdot \nabla u_j^{[\text{int}]}(P_l) $$

which allows for computing the following angles between the $K^{[\text{tot}]}$ column vectors of the matrix:

$$ \beta_{ij} = \arccos \left( \frac{\langle u_i, u_j \rangle}{\sqrt{\langle u_i, u_i \rangle \langle u_j, u_j \rangle}} \right) \quad u_i = u_i^{[\text{int}]}, u_j^{[\text{ext}]} $$

Looking at the values of $\beta_{ij}$ is key to the method. Indeed, $u_j$ vectors are not expected to be strictly orthogonal on the $\{P_i\}$ array: their projections on the sphere ($r = a$) are not located exactly at points of the optimal array and, furthermore, radii of $P_i$ points may vary within 200 kilometers (see next paragraph).

**Remark** We could correct $\vec{B}(P_i)$ for the difference $\vec{B}(P_i) - \vec{B}(Q_i)$ using a priori model. However, here we prefer to present a self-contained algorithm.

## 2.3 A synthetic example

First, we pick up a model $\{g_j, \gamma_j\}$ from the literature (Langlais et al., 2003; $K^{[\text{int}]} = 16$, $K^{[\text{ext}]} = 2$) and compute the values X, Y, Z of the model field at M points which are “real” points in the sense that there are points of the orbits of Champ where real measurements were made (specifically the data points corresponding to 130 days centered at 2003.0). We select N points $P_l$ from the whole set, the one closest to an optimal array of 1000 or 3000 points $\{Q_l\}$, as explained in Fig. 1, and compute back the Gauss coefficients $\{g_j, \gamma_j\}$. The differences between the initial Gauss coefficients and the recovered ones, $g_j$ and $\gamma_j$, $j = 1, \dots, 288$ $k = 1, \dots, 8$, for N = 1000 and N = 3000, are of the order of $10^{-5}$ nT.

Let us take the opportunity with this synthetic example—but with real orbit points—to be more specific about the orthogonality of $u_j$ vectors, taking the case N = 1000

(1000 points in array $\{Q_i\}$). For $K^{[\text{int}]} = 16$ and $K^{[\text{ext}]} = 2$, we find that 99% of couples $u_i, u_j$ vectors make angles $\frac{\pi}{2} \pm \epsilon$ with $\epsilon < 0.01$ radian. The value of $\epsilon$ does not change much with the degree: from 0.005 for the low harmonic couples to 0.03 for the high harmonics. Simulations show that for such values of the departure $\epsilon$ from orthogonality, the coefficients $g_j$ and $\gamma_j$ can be recovered independently of one another with the required accuracy. Such will be the case for all the computations in the paper.

## 2.4 The time sampling of the modeling

A huge advantage of satellite data is their high density (~10⁸ per year!), which allows for a massive number of computations which could not be dreamed of in the pre-satellite times; this density makes it possible to split the data set into many subsets for various applications. In particular, it is possible to compute models at a succession of close time moments; in this paper, we will compute a model every 10 days, at times $t_k = t_0 + k \times 10$ days. Nevertheless, those models are not genuine instantaneous models at $t_k$. Indeed, to get a uniform enough distribution of local times when computing a model of the main field Eq. (2), it is necessary to consider data spanning a time interval $\tau$ around $t_k$; the duration $\tau$ is different for Champ and Oersted. So, our model at time $t_k$ is a model computed from data in a time interval of length $\tau$ centered on $t_k$. But we keep $g(t), h(t)$ (with a daily sampling) in the following temporal series, and not the decimated series with a sampling of $\tau$ days ($\tau = 130$ in case of Champ; see below); a surprising amount of information on short term ($\ll \tau$) features is preserved in the process of computing the model from data in a window of length $\tau$ (e.g., Blanter et al., 2005).

## 2.5 The choice of data amount

We have addressed the question of the number N of data points required to compute a model (an expansion) to degree d with an empirical point of view; indeed, not only the distribution of data points, but also the noise intervene in a broad sense, on the data. We compute a large number (~100) of models $\{g_i, h_j\}$ from suitable random simulations of data; we infer a mean value m and error bar $\sigma$ on each coefficient. Each random simulation at N points is obtained by biasing the initial vector field (calculated from a 16-degree model) by a Gaussian random vector whose distribution parameters were estimated at each point using real data records. And we repeated those computations varying N from 700 to 6000 in 25 steps. For the large majority of coefficients, m is approximately stabilized at
<figure><img src="image_5.png" /><figcaption>Fig. 3. Examples of Gauss coefficients computed as described in the main text. Crosses correspond to estimates obtained with no selection of measurements in the function of magnetic activity, dots to $a_m \le 20$. (a) $g_1^1$, (b) $h_3^2$, (c) $h_3^3$, (d) $g_4^2$.</figcaption></figure>

$N = 2400$ (Fig. 2). The value of $\sigma$ decreases from 0.75 nT for $N = 700$ to 0.2–0.4 nT (depending on the coefficient) for $N = 2400$ and keeps decreasing slowly thereafter. From these results we retained $N = 2400$ in most of the computations of the present paper, for models of degree $\le 16$.

## 2.6 Selecting the data as a function of magnetic activity

We now switch to real data. In most studies, measurements are selected according to the magnetic situation at the time of the measurement, as characterized by the value of the planetary $K_p$ index. We made a series of experiments to evaluate the influence of magnetic activity on our modeling. We take advantage of redundant observations. Let us first retain all of the data, without any selection, versus magnetic activity. For each date $t_k$ (multiple of 10 days) we build a number of quasi-optimal arrays $\{P_i\}$ close to the same optimal array $\{Q_i\}$ (with, in general, $N = 2400$), from different disjoint sets of data; we then compute as many co-efficient sets or models. Some of these appear to be perturbed by big magnetic storms; but it is always possible to find some which are not. The trend of the representative curves allows us to easily discard perturbed values due to the tightness of time sampling. Figure 3(a) has been chosen to illustrate the situation; it represents the evolution of the $g_1^1$ estimate. Crosses are for estimates computed from 130 days of Champ measurements without any selection in the function of magnetic activity. A segment of the curve, in the second half of 2002, is shifted 12 nT below the general trend. And a blank is observed in the first months of 2004. Both segments (S. Maus, personal communication) are characterized by a relatively low number of usable measurement points, which makes it harder to find points close to optimal data sets in the 130 days of data, especially if high magnetic activity is present during those intervals; in

this situation, steps in (g, h) estimate may occur. In fact, due to the abundance of data, it is always possible—except in the case of long gaps—to find close to optimal subsets from which g, h estimates fill the gaps of the graphs of Fig. 3, and are on the general trend. Those conclusions hold for all low degree coefficients of the main field. We kept those gaps here for illustration. Nevertheless, it would be awkward not to take advantage of the large redundancy of data to avoid computations from data corresponding to high activity, e.g., $a_m > 20$ nT (Mayaud, 1980). An extra verification is easily obtained by changing the threshold for $a_m$. Coefficients computed with this condition are represented by dots in Fig. 3 together with those estimated without selection. Along the same lines, we keep the vectorial measurements in high latitudes (as Maus *et al.*, 2006).

# 3. Fitting Oersted and Champ Data

## 3.1 The data

We use Champ data provided by the German team in the form of a list $\{t, X(t), Y(t), Z(t), \vec{r}(t), \vec{r}(t)\}$ being the current point and $t$ being the time of the measurement on the orbit, counted in seconds from June 2001 to December 2004, and Oersted data provided by the Danish team in the same form, $t$ being then a multiple of 1.3 sec (see Stolle *et al.*, 2006), running from March 1999 to June 2003. Note that these data are transformed data. To obtain geocentric components X, Y, Z, Euler angles of the sensor attitude have been determined by the teams in charge. We do not discuss this determination.

## 3.2 Fitting the data

A first opportunity to check the efficiency and accuracy of the algorithm is to look at how coefficients issued from Oersted and Champ fit together, paying special attention to
<figure><img src="image_6.png" /><figcaption>Fig. 4. Gauss coefficients derived from Oersted data (blue points) and Champ data (red points) with $a_m < 20$. (a) $g_1^0$ (b) $h_1^1$, (c) $h_5^3$, (d) $h_5^4$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. (a) Evolution of the moment of the centered dipole $\frac{4}{\mu}\pi a^3 [(g_1^0)^2 + (g_1^1)^2 + (h_1^1)^2]^{1/2}$ computed from Oersted data (blue points) and Champ data (red points). (b) Evolution of the angle $\theta$ between the rotation axis and the axis of the centered dipole.</figcaption></figure>

the overlapping period (June 2001–June 2003). We treat the data from both satellites exactly the same way, except for the time span $\tau$ required to build a model from data reasonably uniformly distributed in local time; indeed, 130 days are required in the case of Champ, while 2 years are needed in the case of Oersted to get the same performance; actually we use a 90-day time span for Oersted to the cost of a less strict condition on the uniformity of local time distribution.

Results are illustrated, for a few Gauss coefficients in Fig. 4. Let us look, for example, at $h_1^1$ (Fig. 4(b)); at each time $t_k$ there are eight estimates for Champ, five for Oersted. We will systematically use this control in all our computations. Clearly, $h_1^1$ values derived from Oersted and Champ agree within 2 nT most of the time, without any further averaging. Other examples are for $h_5^3$ and $h_5^4$ (mind the enlarged scale): the mean values from Oersted and Champ coincide within a nT. In the computations of this section, we

retained data subsets built in such a way that local times are reasonably uniformly distributed. The agreement between Oersted- and Champ-derived coefficients is pretty good but not perfect. There is a physical limitation to this agreement: the ionospheric field is not sampled in exactly the same way by the two satellites. Resulting departures are small and quite variable with the coefficient. For example, for coefficient $g_1^0$, a tiny drift of Oersted estimates with respect to Champ ones is observed Fig. 4(a) (the differences remaining smaller than 2 nT).

We also computed the dipole moment and the angle between the geographical axis and the dipole axis over the 1999–2004 time span (Fig. 5(a)). Except for some trouble at the beginning of Oersted life, the linear decreasing trend appears almost perfect.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Evolution of Gauss coefficients computed from Champ data. Black—all local times, red—day local times (6:00 to 18:00 LT), blue—night local times (18:00 to 6:00 LT). (a) $g_1^0$, (b) $g_2^0$.</figcaption></figure>

## 4. Internal and External Field

The satellite sees as a field of internal origin the sum of the main field generated by the dynamo, the lithospheric field, and the field generated by electric currents flowing in the ionosphere below the perigee of the satellite, especially according to the classical views, in the *E* layer at 110 km altitude (Ratcliffe, 1972). The satellites also encounter electric currents since they are flying in the upper part of the *F* layer. In particular, the fields associated with field-aligned electric currents in the polar regions can be large. But, as discussed in Section 2.6, despite their large magnitude, they do not severely affect the aimed modeling (the computation of internal low ($\le 13$) degree Gauss coefficients from data covering a given time interval). So, we may reasonably assume that the essential part of the field which is external to the solid Earth but internal to the satellites orbits is generated by currents in the *E* layer, i.e., currents driven by the atmospheric dynamo and currents driven in high latitude regions by forces originating high in the ionosphere (the polar current system) (Ratcliffe, 1972; Encrenaz et al., 2004).

To obtain a model of the main field, we have to get rid of this ionospheric field. As is well known, this is not an easy task (Olsen, 1996; Thomson, 2000). It is nevertheless possible to estimate the magnitude of the ionospheric contribution to the internal Gauss coefficients. We will not develop this question at length in the present paper. A classical method is to select data on the basis of local time, e.g., all local times, only day times or only night times (e.g., 6:00–18:00 or 18:00–6:00 LT). Our results will be illustrated by a few graphs. Figure 6(a) shows estimates of $g_1^0$, over a time span of 4 years, derived from Champ data, computed respectively from data at all local times, day times, and night times. The difference between night time estimates and the all times estimates is almost everywhere smaller than 3 nT (2 nT in the second half of 2004). Results for $g_2^0$ are shown in Fig. 6(b) (mind the scale). Generally, the amplitudes of the differences between the different estimates as well as their evolution in time depend on the considered coefficient in a way which is not straightforward to understand. For all of the coefficients, differences between different estimates do not exceed 2 or 3 nT. A last graph Fig. 7 shows two estimates of $g_1^0$ using all local times; for the first one, no selection is made in the function of activity; for the second one, only measurements corresponding to $a_m \le 20$ nT are

<figure><img src="image_9.png" /><figcaption>Fig. 7. Two estimates of $g_1^0$. All local times. Red points—all data, blue points—$a_m \le 20$.</figcaption></figure>

retained. A significant difference of ~2 nT shows up, but the trend is the same.

### 4.1 The internal field

To obtain the main (dynamo) field model, we first take the average of models of the internal field (sources within the sphere $r = r_{perigee}$) computing on all the universal times (practically at 0:00 UT, 1:00 UT to 23:00 UT). That comes down to compute a model I from points uniformly distributed in longitude; to check it we choose points $Q_i$ (see Section 2.1) whose measurements correspond to local times drawn randomly in the interval 0:00 LT–24:00 LT and compute a model II derived from this set of points. Figure 8(c) shows that the two models are indeed identical. Of course, adopting such a model for the main field means that the 24-h averaged ionospheric field is supposed to be zero—when averaged over a full day. It is such a model that we compute at days $t_k = t_0 + k \times 10$ days, from data in a 130-day time span $\tau$ centered on $t_k$.

We compare our two models, relative to 2003.0, with the POMME-3 model (Maus *et al.*, 2006) relative to the same epoch, in two ways. First, we compare the coefficients of our models and those of the POMME-3 model by computing their differences as well as the mean and standard deviation of these differences (Fig. 8(a) and 8(b)). Except for $g_1^0$ and $g_3^0$, the models can be said to be very similar. Our computation of the main field supposes we eliminate the ionospheric field by averaging in longitude, i.e. that this latter field has no zonal component, or, for approximate sym-
<figure><img src="image_10.png" /><figcaption>Fig. 8. (a) Differences between the coefficients of POMME-3 model and our model I computed from data uniformly distributed in local times. Standard deviation 0.35 nT, 0.17 nT when removing $g_1^0, g_3^0$. (b) Differences between coefficients of POMME-3 model and our model II computed as the average of 24 UT models. Standard deviation 0.32 nT, 0.16 nT. (c) Differences between coefficients of models I and II. Standard deviation 0.16 nT. Lexicographic ranking of Gauss coefficients is used.</figcaption></figure>

metry reasons, no component $g_{2k+1}^0$. Choosing only night values leads to slightly different estimates (Fig. 6(a)). We then draw the maps of the vertical components of models I and POMME-3 at the core-mantle boundary (CMB) more precisely on a sphere of radius 3480 km. We also map their differences (Fig. 9). We note that, despite the smallness of difference in coefficients illustrated by Fig. 8, a few small scale anomalies may reach a notable amplitude. The geometrical factor $(\frac{6350}{3480})^{14}$ is indeed equal to 4500; the downward continuation to the CMB of a core field contaminated by crustal anomalies is known to require some precautions. As expected, the largest discrepancies are observed in high latitudes, due to field aligned currents.

## 4.2 The external (ring current) field

We keep here a rather formal point of view, without addressing the physical nature of the field. This is the way the question has been treated for decades—the external field was globally called the ring current field up to recent modelings of satellite data (Olsen *et al.*, 2000). Only external field coefficients of the first degree were considered to be safely determined. So we compute, with the same sampling interval (every 10 days), from the same sets of data, coefficients (we note them $\gamma$, $\eta$ for simplicity, $\gamma$ for the cosine and $\eta$ for sine term) of the external field, i.e., external to the sphere containing the satellites’ orbits ($r > r_{apogee max}$). It is useful to recall again that we are not computing an instantaneous
<figure><img src="image_11.png" /></figure>

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

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

Fig. 9. Map of the radial field at the CMB ($r = 3480$ km), top: our model computed from 130 days of Champ data centered on 2003.0; middle: POMME-3, Potsdam Magnetic Model of the Earth, centered on 2003.0; bottom: Difference of models. Model degrees 1–13 are used in both cases.

field, but, at day $t_k$, a field based on a data set extracted from 130 days of data centered at $t_k$ (for Champ). For each data set we check that the orthogonality conditions required in Section 2 are verified: the $\nabla u_i^{[ext]}$ are orthogonal, with a high accuracy, to one another, and orthogonal to all the $\nabla u_i^{[int]}$ (Section 3). There is no contamination of the external coefficients by the internal ones, i.e., no contamination of the external field by the internal field. The results depend on what is being looked for. For example, $\gamma_1^0$ value is much larger when computed only from night (18:00–06:00 LT)

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

Fig. 10. The evolution of $\gamma_1^0$ coefficient. Blue points: Oersted data. Red points: Champ data. Black curve: $D_{st}$ index 90 days averaged.

values than when computed from day (06:00–18:00) values. This asymmetry was pointed out as early as in 1970 by Olson (1970).

We retain all the local times—which comes down to averaging in longitude—and focus on the axisymmetric coefficients $\gamma_k^0$. The graph of Fig. 10 represents the variation from April 1999 to December 2004 of $\gamma_1^0$. Data from both satellites, Champ and Oersted, are used. Again, the fit of the models derived from the two data sets is excellent, within 1 or 2 nT. $\gamma_1^0$ displays variations with time constants of a few months and amplitudes of some 15 nT. At the bottom of Fig. 10 the evolution of $D_{st}$ index is presented, averaged over a running window of 90 days to make the two graphs comparable (such a comparison is not a new idea (e.g. Cain *et al.*, 1967). Let us recall that $D_{st}$ (Sugiura, 1964) is “the disturbance field which is axially symmetric with respect to the dipole axis, and which is regarded as a function of storm time”. The correlation between $D_{st}$ and $\gamma_1^0$ is good up to the end of 2002; amplitudes of $\gamma_1^0$ variations are smaller. An interesting observation concerns the base level; it is zero by construction for $D_{st}$ index, while $\gamma_1^0$ evolves between 15 and 40 nT.

Maus and Lühr recently performed a study of the magnetospheric field during magnetically quet times (Maus and Lühr, 2005) using Oersted and Champ data from the years 1999–2004. The field is decomposed into contributions from sources in the solar-magnetic frame, and those in the geocentric-solar-magnetospheric frame. Such a separation is probably necessary for a coherent study of the external field. We pointed out at the beginning of this section the limited scope of our study of the “external field”.

## 5. Conclusion

The general objective of the analysis we presented here is to model different ingredients of the field altogether. The method reported here is the realization. This paper also has also a methodological character; we computed Gauss coefficients of the main (dynamo) field, and computed coefficients of the so-called ring current field. The characteristics of the analysis are the following. To compute a given model of the main field, for example, it is possible to use for each $\tau$-interval centered on day $t$ a number of disjoint data sets, each including a rather small (1000, 2400) number of
points. In this way, spurious values of coefficient estimates are made conspicuous; this provides a control of the model (g, h) which allows us to release data selection. An advantage of a close time-spacing of the model is to provide time series whose trends can be studied in the usual way.

**Acknowledgments.** We acknowledge fruitful discussions with S. Maus. We deeply thank the two referees, S. Macmillan and H. Utada, for the attention they paid to the first version of the paper and the resulting numerous and helpful comments and recommendations.

## Appendix A. The Numerical Integration Over Spheres

The question of distributing uniformly N points on a surface—we consider here the case of a sphere—is by no means a trivial one; it has been of interest to mathematicians since Antiquity, and is still the object of research (Hardin and Saff, 2004).

We want an algorithm which distributes a set of N points on a sphere in such a way that the distribution of these points converges to an uniform distribution when N gets large. But uniformity makes sense only for infinite sets. For N finite, we look for a configuration which is optimal with respect to some property, but which may be poor with respect to another property.

For example, we may want to approximate the integral of a function over the sphere $S^2$ by an arithmetic sum (without weights) of values $f(\vec{r}_k)$ at some N well-chosen points $\vec{r}_k$ on $S^2$, in other words we want the difference

$$ \left| \frac{1}{4\pi} \int_{S^2} f(\vec{r}) d\sigma(\vec{r}) - \frac{1}{N} \sum_{k=1}^{N} f(\vec{r}_k) \right| \quad (A.1) $$

to be small for a large class of functions. Such a configuration of points will be called optimal for the evaluation of the integral. The problem has not yet received a general precise solution, but explicit particular solutions have been found. For example Delsarte et al. (1977) considered so-called {$N, t$} spherical designs: configurations of N points $\vec{r}_k$ such that for all polynomials $P_m(\vec{r})$ (i.e., polynomials in three variables $(x, y, z) = \vec{r}$) of degree $m \le t$, the difference in Eq. (A.1) is equal to zero. By means of a computer search, Hardin and Sloane (1996) found spherical designs for all $t \le 13$ with a minimal number of points; for instance, they produced a (94, 13)-design. A variety of algorithms have also been proposed for explicitly constructing asymptotically uniform distributions of points on $S^2$. The most recent approach is to look for configurations corresponding to the minimum potential energy of N repelling points (repelling force needs to be redefined). The problem with the minimum energy configurations algorithm is that it is long and cumbersome. That is why we consider in this paper much more straightforward “spiral sets” which are almost as good as the “polynomial adapted” or “energy sets” referred to above. They are also good with respect to the property of orthogonality of spherical harmonics; this essential property will be simply checked.

Let us now describe the spiral set devised by Rakhmanov, Saff, and Zhou (Rakhmanov et al., 1994), used in the present study. In spherical coordinates $(\theta, \varphi)$, $0 \le \theta \le \pi$,

0 $\le \varphi \le 2\pi$, we take the following coordinates of the N points as:

$$ \theta_k = \arccos(h_k), \quad h_k = -1 + \frac{2k-2}{N-1}, \quad 1 \le k \le N \\ \varphi_k = \left( \varphi_{k-1} + \frac{\text{const}}{\sqrt{N}} \cdot \frac{1}{\sqrt{1-h_k^2}} \right) (\text{mod } 2\pi), \quad 1 < k < N; \\ \varphi_1 = \varphi_N = 0 \quad (A.2) $$

The estimate of the maximum diameter of N non-overlapping disks on the spherical surface (the so-called “best packing argument” (Habicht and van der Waerden, 1951)) suggests the constant in Eq. (A.2) to be chosen such that

$$ (\varphi_k - \varphi_{k-1}) \sqrt{1 - h_k^2} \approx \sqrt{\frac{8\pi}{N\sqrt{3}}} \approx \frac{3.8}{\sqrt{N}} $$

The following version of the latter construction with a good choice of the constant is easy to remember: we first generate N points $(x, y)$ belonging to the unit square $0 \le x, y \le 1$ and then use the cylindrical equal-area projection (i.e., $\theta = \arcsin(2y-1)$, $\varphi = 2\pi x$) onto the spherical surface. The generation of the initial sequence in the square is also simple: $x_k = \{(k-1) \cdot \phi\}$ (here braces denote the non-integer part of the real value) and $y_k = (k-1)/(N-1)$. The good value of the constant in Eq. (A.2) corresponds to $\phi = \frac{1}{2}(\sqrt{5}-1)$ (the reciprocal of the so-called “golden number”). This spiral set is illustrated in the main text for $N=1000$ and $N=3000$ (Fig. 1).

Let the surface harmonic functions $P_n^m(\cos\theta) \cos m\varphi$, $P_n^m(\cos\theta) \sin m\varphi$ be ranked in the usual lexicographic order, let $u_j$ be the jth surface harmonic in the corresponding series. Their gradients $\nabla u_i$, $\nabla u_j$ are orthogonal to one another on the set of N points $\vec{r}_l, l = 1, \dots, N$:

$$ C_i \delta_{ij} = (u_i, u_j) = \sum_l \nabla u_i(\vec{r}_l) \cdot \nabla u_j(\vec{r}_l) \quad (A.3) $$

This ensures that computing Gauss coefficients $g_n^m, h_n^m$ in the main text through non-weighted averages taken on the points $\vec{r}_l$ of the spiral set ($Q_l$ in the main text) is valid. To keep a practical point of view we checked that this orthogonality property is satisfied with a sufficient accuracy for our needs in the present study (see main text).

As we do not know of any other algorithm as simple to operate as this spiral set one, while producing better results for the problem at hand, we did not extend the analysis to other constructions.

## References

* Blanter, E., M. Shnirman, and J. L. Le Mouël, Solar variability, evaluation of correlation properties, *J. Atmos. Terr. Phys.*, **67**, 521–534, 2005.

* Cain, J. C., S. J. Hendricks, R. A. Langel, and W. V. Hudson, A proposed model for the International Geomagnetic Reference Field-1965, *J. Geomag. Geoelectr.*, **19**, 335–355, 1967.

* Delsarte, P., J. Goethais, and J. Seidel, Spherical codes and designs, *Geom. Dedicata*, **6**, 363–388, 1977.

* Encrenaz, T., J. Bibring, M. Blanc, M. Barucci, P. Zarka, and F. Roques, Le système solaire, *CNRS editions*, 2004.
Habicht, W. and B. L. van der Waerden, Lagerung von Punkten auf der Kugel, *Math. Ann.*, **123**, 223–234, 1951.

Hardin, D. P. and E. B. Saff, Discretizing manifolds via minimum energy points, *Notices of AMS*, **51**(N10), 1186–1194, 2004.

Hardin, R. and N. Sloane, McLaren’s improved snub cube and other new spherical designs in three dimensions, *avmalgin. livejournal.com/1203330.html*, *Disc. Comp. Geom.*, **15**, 429–442, 1996.

Langlais, B., M. Mandea, and P. Ultre-Guerard, High-resolution magnetic field modeling: application to Magsat and Ørsted data, *Phys. Earth Planet. Inter.*, **135**, 77–91, 2003.

Maus, S. and H. Lühr, Signature of the quiet-time magnetospheric magnetic field and its electromagnetic induction in the rotating Earth, *Geo-phys. J. Int.*, **162**, 755–763, 2005.

Maus, S., M. Rother, C. Stolle, W. Mai, S. Choi, H. Lühr, D. Cooke, and C. Roth, Third generation of the Potsdam Magnetic Model of the Earth (POMME), *G³*, **7**, N7, 2006.

Mayaud, P. N., Derivation, meaning and use of geomagnetic indices, *Geo-physical Monograph*, **22**, AGU, Washington D.C., 1980.

Olson, W. P., Variations in the Earth’s surface magnetic field from the magnetopause current system, *Planet. Space Sci.*, **18**, 1471–1484, 1970.

Olsen, N., A new tool for determining ionospheric currents from magnetic satellite data, *Geophys. Res. Lett.*, **23**, 3635–3638, 1996.

Olsen, N., T. J. Sabaka, and L. Toffner-Clausen, Determination of the IGRF 2000 model, *Earth Planets Space*, **52**, 1175–1182, 2000.

Olsen, N., H. Lühr, T. J. Sabaka, M. Mandea, M. Rother, L. Toffner-

Clausen, and S. Choi, *CHAOS*—a model of the Earths magnetic field derived from CHAMP, Oersted, and SAC-C magnetic satellite data, *Geophys. J. Int.*, **166**, 67–75, 2006.

Rakhmanov, E., E. Saff, and Y. Zhou, Minimal discrete energy on the sphere, *Math. Res. Lett.*, **1**, 647–662, 1994.

Ratcliffe, J. A., *An Introduction to the Ionosphere and Magnetosphere*, Cambridge University Press, 1972.

Sabaka, T. J., N. Olsen, and R. A. Langel, A comprehensive model of the quiet-time near-Earth magnetic field: Phase 3, *Geophys. J. Int.*, **151**, 32–68, 2002.

Sabaka, T. J., N. Olsen, and M. Purucker, Extending comprehensive models of the Earth’s magnetic field with Oersted and CHAMP data, *Geo-phys. J. Int.*, **159**, 521–547, 2004.

Stolle, C., H. Lühr, M. Rother, and G. Balasis, Magnetic signatures of equatorial spread F as observed by the CHAMP satellite, *J. Geophys. Res.*, **111**, A02304, 2006.

Sugiura, M., Hourly values of equatorial $D_{st}$ for IGY, in *Annals of the International Geophysical Year*, **35**, 945–948, Pergamon Press, Oxford, 1964.

Thomson, A. W. P., Improving the modelling of the geomagnetic main-field: Isolating the average ionospheric field in satellite data, *Earth Planets Space*, **52**, 1199–1206, 2000.

J. L. Le Mouël, P. Shebalin, and A. Khokhlov (e-mail: khokhlov@mitp.ru)