LABEL:firstpage–LABEL:lastpage

# Secular Variation Signals in Magnetic Field Gradient Tensor Elements derived
from satellite-based Geomagnetic Virtual Observatories

Magnus D. Hammer  Christopher C. Finlay and Nils Olsen
Division of Geomagnetism  Technical University of Denmark - DTU Space  Kgs.
Lyngby  Denmark
E-mail<EMAIL_ADDRESS>

(Received XXX; in original form XXX; xxxx)

###### keywords:

Satellite magnetics, Rapid time variations, Magnetic field variations through
time

††volume: XXX

We present local time series of the magnetic field gradient tensor elements at
satellite altitude derived using a Geomagnetic Virtual Observatory (GVO)
approach. Gradient element time series are computed in four-monthly bins on an
approximately equal-area distributed worldwide network. This enables global
investigations of spatio-temporal variations in the gradient tensor elements.
Series are derived from data collected by the Swarm and CHAMP satellite
missions, using vector field measurements and their along-track and East-West
differences, when available. We find evidence for a regional Secular Variation
impulse (jerk) event in 2017 in the first time derivative of the gradient
tensor elements. This event is located at low latitudes in the Pacific region.
It has a similar profile and amplitude regardless of the adopted data
selection criteria and is well fit by an internal potential field. Spherical
harmonic models of the internal magnetic field built from the GVO gradient
series show lower scatter in near-zonal harmonics compared with models built
using standard GVO vector field series. The GVO gradient element series are an
effective means of compressing the spatio-temporal information gathered by
low-Earth orbit satellites on geomagnetic field variations, which may prove
useful for core flow inversions and in geodynamo data assimilation studies.

##  1 Introduction

The main part of the geomagnetic field is generated in the Earth’s fluid outer
core by a process known as the geodynamo. Knowledge of how this core field
varies with space and time provides information on the fluid flow dynamics in
the liquid metal outer core. Although the temporal behaviour of the
geomagnetic field is well characterized in time series from ground
observatories, a global spatio-temporal study is hampered by the uneven
distribution of these observatories. Even though low-Earth-orbit (LEO)
satellites do provide good global coverage on timescales of weeks and longer,
the direct study of the first time derivative of the core field, the secular
variation (SV), from satellite measurements is not straightforward. This is
because LEO satellites are not geostationary (i.e. do not have the same
orbital period as the Earth’s rotation, such that their position does not stay
fixed as seen from ground), which leads to an ambiguity, since it is not
possible to establish whether an observed field variation is caused by a
temporal or spatial change of the field (Olsen & Stolle, 2012) . Spherical
harmonic (SH) field models derived from satellite measurements provide an
established way of studying the SV field and its time derivative, the secular
acceleration (SA), globally. However, such harmonic functions have global
support, which means that an estimate of the SV at a specific position may be
affected by noise from remote locations.

These issues lead Mandea & Olsen (2006) to introduce the concept of
Geomagnetic Virtual Observatories (GVOs) in space, in which satellite magnetic
measurements from within a selected region, collected during one month time
windows, were used to derive a local monthly mean vector field at the
satellite mean altitude. The resulting GVO time series resemble monthly mean
series computed using ground observatory magnetic measurements, by providing
the magnetic vector field elements at fixed locations. However, since they are
based on satellite data, regular sampling in both space and time is possible.
Thus, the GVO method provides a tool for compressing satellite magnetic field
measurements into a dataset that contains time series distributed in a global
grid, and which may also be supplied with error estimates (Hammer et al.,
2021c). Olsen & Mandea (2007) used CHAMP measurements to derive GVO vector
component time series, and carried out a global investigation of SV that
identified a regional geomagnetic jerk event in 2003.

In the original GVO approach of Mandea & Olsen (2006), processing of the
satellite measurements followed that of simple monthly field means at ground
observatories, taking measurements from all local times and with all levels of
geomagnetic activity, and relied on the assumption that short period external
fields would have zero mean over the course of one month. However, later
studies revealed that external fields, especially those due to the
magnetospheric ring current and ionospheric current systems, cause
contamination of the retrieved internal GVO field signal (Olsen & Mandea,
2007; Beggan et al., 2009; Domingos et al., 2019). In addition, insufficient
local time sampling from within one month of polar orbiting satellites
resulted in a bias due to the local time dependence of ionospheric and
magnetosphere-ionosphere coupling currents (Shore, 2013). Recently, the GVO
processing algorithm has been further developed in an effort to reduce
contamination from magnetospheric and ionospheric sources, and the local time
sampling bias, with the aim of better isolating the field signal generated by
the Earth’s outer core (Hammer et al., 2021a). These new GVO vector field
series have been used to study global patterns of field changes (Hammer et
al., 2021a, c), for inferring fluid flows close to the core surface (Kloss &
Finlay, 2019; Rogers et al., 2019) and for data assimilation studies (Barrois
et al., 2018; Huder et al., 2020).

In parallel to the development of these GVO-based techniques there has also
been recent progress in the theory of space-based magnetic gradiometry,
inspired by advances in satellite gravimetry. Initial studies have
demonstrated that knowledge of the second-order $3\times 3$ magnetic gradient
tensor may be beneficial when seeking to retrieve small scale features of the
field (both the lithospheric field and the time-dependent core field). This is
possible because gradient elements effectively give more weight to shorter
wavelengths, while at the same time some noise sources (e.g. unmodeled
magnetospheric fields) are predominantly of long wavelength, which can result
in a higher signal-to-noise ratio for short wavelengths compared to using the
vector field components (Kotsiaros & Olsen, 2012, 2014). Information on these
shorter spatial length scales is crucial for core dynamic studies, yet their
robust determination is a major obstacle in core field modelling. This issue
has motivated us to derive time series of field gradient elements using the
GVO method, and to investigate whether such series allow, via for example
spherical harmonic analysis, an improved retrieval of field changes. Gradient
GVO series have the potential to also be used directly in core flow inversions
and data assimilation studies of the geodynamo.

Assuming a potential field due to internal sources and no in-situ electrical
currents, the field becomes a solenoidal irrotational vector field and the
gradient tensor has the special property of being symmetric with a trace of
zero. The assumption of a symmetric gradient tensor reduces the number of
independent gradient tensor elements from nine to six, while a trace of zero
further reduces this number to five. Each element of the magnetic gradient
tensor may be considered as a directional filter providing specific
information on the magnetic field structures. Thereby, certain gradient tensor
elements better constrain specific spherical harmonics (Olsen & Kotsiaros,
2011). According to the studies of Kotsiaros & Olsen (2012) and Kotsiaros &
Olsen (2014), knowledge of the radial gradient of the radial field, written as
$\left[\nabla B\right]_{rr}$, is particularly suitable for resolving the
higher degree parts and zonal harmonics. The East-West gradient of the
azimuthal field, $\left[\nabla B\right]_{\phi\phi}$, and radial field,
$\left[\nabla B\right]_{r\phi}$, are especially sensitive towards sectorial
harmonics, while the North-South gradient of the radial, $\left[\nabla
B\right]_{r\theta}$, and meridional, $\left[\nabla B\right]_{\theta\theta}$,
fields are especially useful for determining near-zonal harmonics. The East-
West gradient of the meridional field, $\left[\nabla B\right]_{\theta\phi}$
does not provide significant additional information. In addition, knowledge of
how certain external fields may influence certain gradient tensor elements is
important to consider, for instance the magnetospheric ring current is
expected to affect zonal terms constrained by the $\left[\nabla B\right]_{rr}$
element but not the $\left[\nabla B\right]_{r\phi}$ element (Kotsiaros &
Olsen, 2014). Although it is not yet possible to directly measure the full
magnetic gradient tensor in space (Nogueira et al., 2015), it is nonetheless
possible to compute the tensor elements from a magnetic potential determined
using satellite magnetic measurements.

In this paper, our aim is to estimate local time series of the magnetic field
gradient tensor elements using the GVO method. Our primary motivation for
deriving such gradient field series, is to improve the recovery of the small
length-scales of the field’s secular variation, compared with the traditional
use of vector field data. In particular, the gradient elements are expected to
enable a higher signal-to-noise ratio as compared to vector elements
(Kotsiaros & Olsen, 2014). This is because they are more sensitive to the
small length scales of the field, and less affected by large scale external
fields. We use satellite magnetic measurements to derive the GVO gradient
series, and follow Hammer et al. (2021a) in implementing dark and quiet time
data selection criteria and use 4-month time windows to minimize problems
related to local time sampling.

In Section 2 we provide a detailed description of the satellite magnetic
measurements and selection criteria used, and in Section 3 we describe the GVO
method and the computation of both GVO vector and GVO gradient element time
series. In Section 4.1 we present results of the GVO series for each of the SV
gradient elements, and visually inspect these. In order to investigate the
possible benefits of using GVO gradient data, we compare SH field models
derived epoch by epoch from the GVO vector data and GVO gradient tensor data
in Section 4.2. We study the detailed behaviour of the gradient tensor
elements going from 2015 to 2018 with a focus on the Pacific region. Finally,
Section 5 provides discussions and conclusions.

##  2 Data

To derive the GVO time series we select vector magnetic field measurements
from the CHAMP and Swarm satellite missions. We used CHAMP L3 magnetic data
between July 2000 and September 2010 and Swarm Level 1b MAG-L, version
0505/0506, from all three Swarm satellites Alpha, Bravo and Charlie between
January 2014 and April 2020, and sub-sample at 15s intervals (i.e. taking
measurements every 15 second), in the vector field magnetometer (VFM) frame.
Next, the magnetic data in the VFM frame are rotated into an Earth-Centered
Earth-Fixed (ECEF) local Cartesian North-East-Centre (NEC) coordinate frame
(for details see Olsen et al. (2006)) by using the Euler rotation angles from
the CHAOS-7.2 model (Finlay et al., 2020). Measurements from known problematic
days (e.g. where satellite manoeuvres took place) were removed and gross data
outliers for which the vector field components deviated more than 500nT from
CHAOS-7.2 field model [http://www.spacecenter.dk/files/magnetic-
models/CHAOS-7/], see also Finlay et al. (2020) values were rejected. The
measurements were then selected using a dark quiet time criteria defined here
as: a) the sun is at least $10^{\circ}$ below horizon, b) geomagnetic activity
index $K_{p}<3^{\circ}$, c) ring current index $|dRC/dt|<3\mathrm{nThr^{-1}}$
(Olsen et al., 2014), merging electric field (averaged over two-hours) at
magnetopause $E_{m}\leq 0.8\mathrm{mVm^{-1}}$ (Olsen et al., 2014), and
placing constraints on the interplanetary magnetic field (IMF) (averaged over
two-hours) requiring $B_{z}>0\mathrm{nT}$ and $|B_{y}|<10\mathrm{nT}$ (Ritter
et al., 2004). Using the definition in Olsen et al. (2014), the merging
electric field was derived using 1-min values of the solar wind, solar clock
angle and IMF extracted from the OMNI database, http://omniweb.gsfc.nasa.gov.

Previous studies have demonstrated the benefits of using along-track
differences of the satellite magnetic field measurements for retrieving higher
spatial resolution of the core and also the lithospheric fields, as such
differences filter out correlated noise caused by external sources (Olsen et
al., 2015; Kotsiaros et al., 2015; Kotsiaros, 2016; Finlay, 2019). However,
from data differences alone it is not possible to obtain robust information on
the longer wavelength part of the field. Therefore, we also use the
complementary data means in a similar manner as in previous studies of the
core signal (Sabaka et al., 2013; Hammer, 2018). From the satellite magnetic
field measurements, $B_{k}(\mathbf{r})$, where $k$ is the vector component of
a given coordinate system, we use measurement means, $\Sigma d_{k}$, and
differences, $\Delta d_{k}$ as data. The differences, $\Delta d_{k}=(\Delta
d_{k}^{\mathrm{AT}},\Delta d_{k}^{\mathrm{EW}})$, and the means, $\Sigma
d_{k}=(\Sigma d_{k}^{\mathrm{AT}},\Sigma d_{k}^{\mathrm{EW}})$, are taken
along-track (AT) for each satellite and East-West (EW) between the Swarm Alpha
(SWA) and Charlie (SWC) satellites. Here along-track differences are
calculated from the 15 s differences $\Delta
d_{k}^{\mathrm{AT}}=[B_{k}(\mathbf{r},t)-B_{k}(\mathbf{r}+\delta\mathbf{r},t+15s)]$
while the means are given by $\Sigma
d_{k}^{AT}=[B_{k}(\mathbf{r},t)+B_{k}(\mathbf{r}+\delta\mathbf{r},t+15s)]/2$.
The East-West differences were calculated as $\Delta
d_{k}^{\mathrm{EW}}=[B_{k}^{\mathrm{SWA}}(\mathbf{r}_{1},t_{1})-B_{k}^{\mathrm{SWC}}(\mathbf{r}_{2},t_{2})]$,
and the means as $\Sigma
d_{k}^{\mathrm{EW}}=[B_{k}^{\mathrm{SWA}}(\mathbf{r}_{1},t_{1})+B_{k}^{\mathrm{SWC}}(\mathbf{r}_{2},t_{2})]/2$.
Considering a given orbit of Swarm Alpha, the corresponding Swarm Charlie
measurement were chosen to be that closest in colatitude provided that
$|\Delta t|=|t_{1}-t_{2}|<50s$ (Olsen et al., 2015).

##  3 Theory and Method

The purpose of this paper is to derive time series of the magnetic field
gradient tensor elements using the GVO method. In section 3.1, we begin by
recalling the GVO method and how this is used to derive time series of the
magnetic field vector elements. Following this, in section 3.2, we describe
our new extension of the GVO method and how this can be used to derive time
series of the magnetic field gradient tensor elements.

###  3.1 Geomagnetic Virtual Observatory Method

The Geomagnetic Virtual Observatory method allows for epoch estimates of the
magnetic vector field components at a given target point (referred to as a GVO
target location) to be derived using available satellite measurements enclosed
within a cylinder of radius 700 km during the course of four months. A radius
of 700 km enables enough data for computing reliable and independent GVO
estimates every four months (Hammer, 2018). From these measurements, provided
in an ECEF coordinate frame given by the spherical polar components,
$\mathbf{B}^{obs}=(B_{r},B_{\theta},B_{\phi})$, magnetic field residuals are
calculated as (Hammer et al., 2021a)

$\delta\mathbf{B}=\mathbf{B}^{obs}-\mathbf{B}^{MF}-\mathbf{B}^{lit}-\mathbf{B}^{mag}-\mathbf{B}^{iono}\;,$
(1)

where model fields subtracted are: the main field (MF), $\mathbf{B}^{MF}$, for
SH degrees $n\in[1,13]$ determined using the CHAOS-7.2 model
[http://www.spacecenter.dk/files/magnetic-models/CHAOS-7/], see also Finlay et
al. (2020); the static lithospheric field, $\mathbf{B}^{lit}$, for SH degrees
$n\in[14,185]$ determined using the LCS-1 model (Olsen et al., 2017); the
large-scale magnetospheric and associated Earth induced fields,
$\mathbf{B}^{mag}$, as given by the CHAOS-7.2 model parameterized in time by
the RC index (Finlay et al., 2020); and the ionospheric and associated Earth
induced fields, $\mathbf{B}^{iono}$, as determined using the CIY4 model
parameterized by 90-day averages of solar flux F10.7 (Sabaka et al., 2018).

Note here that we remove values of the time-dependent main field and then at a
later stage, add back main field values at the GVO target position and target
epoch (Mandea & Olsen, 2006; Hammer et al., 2021a). However, the precise
choice of main field used in both steps is not crucial (Hammer, 2018; Hammer
et al., 2021c). Note that in removing estimates of the time-dependent main
field in eq.(1) estimates of the SV field (up to degree 13) are also
effectively removed for each data point during the 4 month time windows, which
aids a robust estimation. Adding back a main field estimate at the GVO epoch
time then synchronizes the data from the considered four month window to this
common epoch. This way of correcting the individual data includes spatial
gradients of the main field model used, so the resulting local potential does
no longer contain the tensor information from the main field model which has
been removed. It is thus necessary to reinstate the tensor information from
the main field model when computing the magnetic field gradient tensor in
section 3.2. Note that despite implementing a dark quiet time data selection
scheme and also having removed estimates of magnetospheric and ionospheric
fields together with their associated Earth-induced fields, contamination from
non-core electrical currents may persist in the residual GVO field eq.(1).
Such contamination could be related to ionospheric currents such as the polar
electrojets and F-region currents.

Next, for each data point, both the position of the data point and the
corresponding residual magnetic vector, eq.(1), are transformed from the
spherical system to a right-handed local topocentric Cartesian system
$(x,y,z)$ having its origin at the GVO target location, as detailed in
(Hammer, 2018, p. 64). At this specific GVO location (and only at this
location), $x$ points towards geographic south, $y$ points towards east and
$z$ points radially upwards (Hammer et al., 2021a).

At the GVO target point, the unit vectors of the local Cartesian frame
coincide with the spherical polar unit vectors, i.e.
$(\mathbf{\hat{e}_{z}},\mathbf{\hat{e}_{x}},\mathbf{\hat{e}_{y}})=(\mathbf{\hat{e}_{r}},\mathbf{\hat{e}_{\theta}},\mathbf{\hat{e}_{\phi}})$.
Assuming that the magnetic field measurements are made in a source free
region, the residual field, $\delta\mathbf{B}$, is a Laplacian potential field
that is approximately quasi-stationary (Sabaka et al., 2010). This means that
a magnetic scalar potential, $V$, is associated with the residual field, which
in the local Cartesian coordinate system can be expanded as a sum of
polynomials having the form $C_{abc}x^{a}y^{b}z^{c}$ (Backus et al., 1996). In
this application we expand to cubic terms following Hammer et al. (2021a)

$\displaystyle V(x,y,z)=C_{100}x+C_{010}y+C_{001}z+C_{200}x^{2}+C_{020}y^{2}$
$\displaystyle-(C_{200}+C_{020})z^{2}+C_{110}xy+C_{101}xz+C_{011}yz$
$\displaystyle-\frac{1}{3}(C_{102}+C_{120})x^{3}-\frac{1}{3}(C_{210}+C_{012})y^{3}$
$\displaystyle-\frac{1}{3}(C_{201}+C_{021})z^{3}+C_{210}x^{2}y+C_{201}x^{2}z$
$\displaystyle+C_{120}y^{2}x+C_{021}y^{2}z+C_{102}z^{2}x+C_{012}z^{2}y+C_{111}xyz\;.$
(2)

The means and differences of the residual magnetic field are linked to this
potential via appropriate design matrices constructed as described in Hammer
et al. (2021a). The coefficients of the potential are determined from a robust
least-squares solution, which includes a) an a prior data covariance matrix
derived from standard deviations of the residuals between the data (means and
differences) and estimates of an un-weighted least-squares solution, b) a
diagonal weight matrix consisting of robust (Huber) weights, using a scale
constant of 1.5 (e.g., Constable, 1988), and c) an additional down-weighting
factor of 1/2 when data comes from Swarm satellites Alpha and Charlie, taking
into account that these fly side-by-side and thus provide similar
measurements. From these potential coefficients, a mean residual magnetic
field for the given GVO target point position and epoch is computed as
$\delta\mathbf{B}_{GVO}(x=0,y=0,z=0)=-\nabla V=-(C_{100},C_{010},C_{001})$.
This mean residual field is then rotated back into the vector components in
spherical polar coordinates, $\delta B_{GVO,r}=\delta B_{GVO,z}$, $\delta
B_{GVO,\theta}=\delta B_{GVO,x}$, $\delta B_{GVO,\phi}=\delta B_{GVO,y}$ and
afterwards a main field estimate evaluated at the GVO epoch using SH degrees
$n\in[1,13]$ is added back to obtain the GVO field (Hammer et al., 2021a). We
recall here, that the CHAOS-7.2 main field up to degree 13 has to be
reinstated, as time-dependent main field estimates were removed by eq.(1).

Figure 1: Number (No.) of GVOs for each epoch, given the applied dark and
quiet-time selection criteria, during CHAMP (blue) and Swarm (green) times.

Following the same procedure as in Hammer et al. (2021a), we compute a global
grid of 300 GVO’s located on an approximately equal area grid based on the
sphere partitioning algorithm of Leopardi (2006). The global grid of 300 GVOs
are listed starting from a position of latitude $89.9^{\circ}$N at longitude
$0^{\circ}$ going to a position of latitude $89.9^{\circ}$S at longitude
$0^{\circ}$. The distance between the GVOs in this grid is $\approx 1400$ km,
and with a target cylinder radius of 700 km close to 80% of the measurements
are used. This choice of radius avoids overlap between the GVO data cylinders,
such that each satellite measurement is only used once in the global GVO grid
computation, which helps to minimize possible correlated errors and biases due
to the GVO binning method (Beggan et al., 2009; Shore, 2013). The GVO height
above ground is taken to be 370 km and 490 km (approximate mean orbital
altitude) during CHAMP and Swarm times, respectively. Figure 1 presents the
available number of GVOs for each epoch (the maximum possible number at each
epoch is 300). Especially when considering the early CHAMP period, there are
several epochs with very few available GVOs. Such strongly depleted epochs are
primarily caused by our restrictions to geomagnetically quiet and night-time
data, and to selection based on CHAMP data quality flags. Table 1 presents the
Huber-weighted mean and root-mean-square (rms) of residuals between the
satellite measurements (after applying the corrections of eq.(1)) and the
field estimates as computed from the determined potential coefficients
$(C_{100},C_{010},C_{001})$. These residuals are summed for each component and
split into regions of 78 polar and 222 non-polar GVO’s, defining polar to be
GVOs poleward of $\pm 54^{\circ}$ geographic latitude. The polar rms values
for both data sums and differences are higher than the non-polar, and the
CHAMP values are higher than the Swarm values. The non-polar rms values for
all components are below 2nT during both CHAMP and Swarm times.

Table 1: GVO model rms misfit statistics between contributing satellite data and GVO estimates from the total number (No.) of GVOs, as compiled in Figure 1, during CHAMP and Swarm. Here $\sum$ and $\Delta$ represent data means and data differences, respectively. | CHAMP | Swarm |  | CHAMP | Swarm
---|---|---|---|---|---
Component | No. | Mean | rms | No. | Mean | rms | Component | No. | Mean | rms | No. | Mean | rms
|  | [nT] | [nT] |  | [nT] | [nT] |  |  | [nT] | [nT] |  | [nT] | [nT]
Polar | 2574 |  |  | 1872 |  |  | Non-polar | 7326 |  |  | 5328 |  |
$\sum B_{x,NS}$ |  | -0.30 | 6.61 |  | -0.52 | 6.26 | $\sum B_{x,NS}$ |  | -0.80 | 1.75 |  | 0.01 | 1.69
$\sum B_{y,NS}$ |  | 0.00 | 6.52 |  | -0.02 | 6.79 | $\sum B_{y,NS}$ |  | 0.00 | 1.46 |  | 0.00 | 1.74
$\sum B_{z,NS}$ |  | 0.00 | 3.33 |  | 0.01 | 3.02 | $\sum B_{z,NS}$ |  | 0.00 | 1.30 |  | -0.00 | 0.95
$\sum B_{x,EW}$ |  |  |  |  | 0.05 | 5.92 | $\sum B_{x,EW}$ |  |  |  |  | -0.03 | 1.57
$\sum B_{y,EW}$ |  |  |  |  | -0.01 | 6.44 | $\sum B_{y,EW}$ |  |  |  |  | 0.01 | 1.48
$\sum B_{z,EW}$ |  |  |  |  | 0.01 | 2.89 | $\sum B_{z,EW}$ |  |  |  |  | -0.01 | 0.88
$\Delta B_{x,NS}$ |  | -0.01 | 4.35 |  | 0.01 | 3.80 | $\Delta B_{x,NS}$ |  | -0.01 | 0.50 |  | 0.00 | 0.26
$\Delta B_{y,NS}$ |  | -0.01 | 5.20 |  | -0.01 | 4.86 | $\Delta B_{y,NS}$ |  | 0.00 | 0.58 |  | 0.00 | 0.38
$\Delta B_{z,NS}$ |  | 0.01 | 1.61 |  | -0.00 | 1.36 | $\Delta B_{z,NS}$ |  | 0.00 | 0.53 |  | 0.00 | 0.27
$\Delta B_{x,EW}$ |  |  |  |  | 0.10 | 3.17 | $\Delta B_{x,EW}$ |  |  |  |  | 0.10 | 0.51
$\Delta B_{y,EW}$ |  |  |  |  | 0.00 | 3.17 | $\Delta B_{y,EW}$ |  |  |  |  | 0.02 | 0.70
$\Delta B_{z,EW}$ |  |  |  |  | -0.07 | 0.95 | $\Delta B_{z,EW}$ |  |  |  |  | -0.02 | 0.50

###  3.2 Computing the Magnetic Field Gradient Tensor within the GVO framework

We now proceed to formulate the magnetic field gradient tensor and describe
how it transforms from a spherical polar coordinate system to the local
topocentric Cartesian right-handed coordinate system used in the GVO method.
This transformation will allow us to compute GVO time series for the magnetic
field gradient tensor elements in analogy to the concept of GVO vector field
time series.

We begin by expressing the magnetic field gradient tensor in the local
Cartesian system of the GVO method described in Section 3.1. This is given by
(see Appendix A for full details)

$\displaystyle\nabla\mathbf{B}$
$\displaystyle=-\left(\begin{array}[]{lll}\frac{\partial^{2}V}{\partial
z^{2}}&\frac{\partial^{2}V}{\partial x\partial
z}&\frac{\partial^{2}V}{\partial y\partial z}\\\ \frac{\partial^{2}V}{\partial
z\partial x}&\frac{\partial^{2}V}{\partial
x^{2}}&\frac{\partial^{2}V}{\partial y\partial x}\\\
\frac{\partial^{2}V}{\partial z\partial y}&\frac{\partial^{2}V}{\partial
x\partial y}&\frac{\partial^{2}V}{\partial
y^{2}}\end{array}\right)=\left(\begin{array}[]{lll}\left[\nabla
B\right]_{zz}&\left[\nabla B\right]_{zx}&\left[\nabla B\right]_{zy}\\\
\left[\nabla B\right]_{xz}&\left[\nabla B\right]_{xx}&\left[\nabla
B\right]_{xy}\\\ \left[\nabla B\right]_{yz}&\left[\nabla
B\right]_{yx}&\left[\nabla B\right]_{yy}\end{array}\right)\;.$ (9)

This is a second-order tensor where the minus sign comes from defining the
magnetic field as the negative gradient of the potential. The gradient tensor
elements are denoted here by $\left[\nabla B\right]_{jk}$, where the first
subscript, $j$ denotes the vector component under consideration and the second
subscript, $k$, denotes direction of the field derivative. Using the local
cubic potential eq.(2) estimated from the residual magnetic field as described
in Section 3.1, a second-order residual field gradient tensor at the GVO
target point can be derived using eq.(9) as

$\displaystyle\nabla\delta\mathbf{B}_{GVO}$
$\displaystyle=\left(\begin{array}[]{lll}2(C_{200}+C_{020})&-C_{101}&-C_{011}\\\
-C_{101}&-2C_{200}&-C_{110}\\\
-C_{011}&-C_{110}&-2C_{020}\end{array}\right)\;.$ (13)

Because the magnetic field is a solenoidal vector field, the divergence is
zero, such that the trace of the gradient tensor vanishes, i.e.
$tr(\nabla\delta\mathbf{B}_{GVO})=2(C_{200}+C_{020})-2C_{200}-2C_{020}=0$,
reducing the number of independent elements from 9 to 8. In addition to this,
because the field is a Laplacian potential field, the curl of the field
vanishes and the number of independent tensor elements reduces to 5; in other
words, the magnetic gradient tensor is diagonally symmetric and its trace is
zero (Kotsiaros & Olsen, 2012).

Following the GVO algorithm, gradient tensor estimates from a field model then
have to be added back at the GVO target point for harmonic degrees $n\leq 13$.
The time-dependent CHAOS-7.2 main field up to degree 13 has here to be
reinstated, as it was removed by eq.(1), in order to synchronize the data to
the particular GVO epoch. Removal of the main field estimates results in the
local potential in eq.(2) also lacking the main field tensor information; this
must therefore be added back, in order to obtain the full magnetic field
gradient tensor. To do this, we will need to consider how the gradient tensor
elements in spherical polar and Cartesian coordinate systems are related. The
magnetic gradient tensor elements as expressed in the spherical coordinate
system are given by (Olsen & Kotsiaros, 2011; Kotsiaros & Olsen, 2012), see
also Appendix A

$\displaystyle\nabla\mathbf{B}$
$\displaystyle=\left(\begin{array}[]{lll}-\frac{\partial^{2}V}{\partial
r^{2}}&-\frac{1}{r}\frac{\partial^{2}\partial V}{\partial\theta\partial
r}+\frac{1}{r^{2}}\frac{\partial
V}{\partial\theta}&-\frac{1}{r\mathrm{sin}\theta}\frac{\partial^{2}V}{\partial\phi\partial
r}+\frac{1}{r^{2}\mathrm{sin}\theta}\frac{\partial V}{\partial\phi}\\\
-\frac{1}{r}\frac{\partial^{2}V}{\partial
r\partial\theta}+\frac{1}{r^{2}}\frac{\partial
V}{\partial\theta}&-\frac{1}{r^{2}}\frac{\partial^{2}V}{\partial\theta^{2}}-\frac{1}{r}\frac{\partial
V}{\partial
r}&-\frac{1}{r^{2}\mathrm{sin}\theta}\frac{\partial^{2}V}{\partial\phi\partial\theta}+\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}^{2}\theta}\frac{\partial
V}{\partial\phi}\\\
-\frac{1}{r\mathrm{sin}\theta}\frac{\partial^{2}V}{\partial
r\partial\phi}+\frac{1}{r^{2}\mathrm{sin}\theta}\frac{\partial
V}{\partial\phi}&-\frac{1}{r^{2}\mathrm{sin}\theta}\frac{\partial^{2}V}{\partial\theta\partial\phi}+\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}^{2}\theta}\frac{\partial
V}{\partial\phi}&-\frac{1}{r^{2}\mathrm{sin}^{2}\theta}\frac{\partial^{2}V}{\partial\phi^{2}}-\frac{1}{r}\frac{\partial
V}{\partial
r}-\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}^{2}\theta}\frac{\partial
V}{\partial\theta}\end{array}\right)$ (17)
$\displaystyle=\left(\begin{array}[]{lll}\left[\nabla
B\right]_{rr}&\left[\nabla B\right]_{r\theta}&\left[\nabla B\right]_{r\phi}\\\
\left[\nabla B\right]_{\theta r}&\left[\nabla
B\right]_{\theta\theta}&\left[\nabla B\right]_{\theta\phi}\\\ \left[\nabla
B\right]_{\phi r}&\left[\nabla B\right]_{\phi\theta}&\left[\nabla
B\right]_{\phi\phi}\end{array}\right)\;.$ (21)

Here the first column of the tensor contains the derivatives of the magnetic
field components along the radial direction, the second column contains the
derivatives along the co-latitudinal direction and the third column contains
the derivatives along the longitudinal direction. The gradient element in the
first column and row contains one term only, the field derivative term e.g.
$\partial^{2}/\partial r^{2}$, while the other gradient tensor elements in
addition to this also have an additional field term i.e. $\partial/\partial
r$, $\partial/\partial\theta$ or $\partial/\partial\phi$. The transformations
relating the gradient tensor elements in the local Cartesian system to the
tensor elements of the spherical system, only at the GVO target location, are
in the end given by the following simple relations (see Appendix A for a full
derivation).

$\displaystyle\left[\nabla B\right]_{zz}$ $\displaystyle=\left[\nabla
B\right]_{rr}\qquad\left[\nabla B\right]_{zx}=\left[\nabla
B\right]_{r\theta}\qquad\left[\nabla B\right]_{zy}=\left[\nabla
B\right]_{r\phi}$ (22) $\displaystyle\left[\nabla B\right]_{xz}$
$\displaystyle=\left[\nabla B\right]_{\theta r}\qquad\left[\nabla
B\right]_{xx}=\left[\nabla B\right]_{\theta\theta}\qquad\left[\nabla
B\right]_{xy}=\left[\nabla B\right]_{\theta\phi}$ $\displaystyle\left[\nabla
B\right]_{yz}$ $\displaystyle=\left[\nabla B\right]_{\phi r}\qquad\left[\nabla
B\right]_{yx}=\left[\nabla B\right]_{\phi\theta}\qquad\left[\nabla
B\right]_{yy}=\left[\nabla B\right]_{\phi\phi}\;.$

Having determined the potential from the residual magnetic field eq.(1), we
can compute a residual field gradient tensor by eq.(13) and add back main
field gradient tensor estimates from the CHAOS-7.2 field model using eq.(21)
for SH degrees $n\leq 13$, using the above relations, in order to obtain the
required GVO field gradient estimates $\nabla\mathbf{B}_{GVO}$. Note that this
procedure is analogous to the procedure applied in deriving vector field GVOs
where the main vector field is added back. The above procedure is then
repeated at each GVO location and for each epoch to compute all the desired
GVO field gradient time series.

Field contamination from non-core sources dominates the instrumental errors of
the satellites. Thus, the primary limitation in obtaining core field GVOs are
contributions due to unmodelled fields, such as the polar electrojet, and not
measurement errors (e.g. Finlay et al., 2017). Error estimates for each tensor
element $jk$, and separately for CHAMP and Swarm, are computed using the
residuals $e_{jk}=d^{GVO}-d^{CHAOS}$ between the GVO gradient tensor data,
$d^{GVO}=\left[\nabla B_{GVO}\right]_{jk}$, and the gradient element values of
the CHAOS-7.2 for SH degree $n=1-16$, $d^{CHAOS}=\left[\nabla B\right]_{jk}$.
Considering all epochs for each GVO in the grid, the error estimates for
tensor element $jk$ are given by the total mean square error
$\sigma_{jk}=\sqrt{\sum_{i}(e_{jk,i}-{\mu_{jk}})^{2}/M+\mu_{jk}^{2}}$ (e.g.
Bendat & Piersol, 2010), where $e_{jk,i}$ is the residual of the $i$th data
element, $M$ is the number of data in a given series and $\mu_{jk}$ is the
residual mean for a given component. Hammer et al. (2021a) computed similar
uncertainty estimates for the vector components using the vector field
residuals with respect to the CHAOS-7.2 model. We note here, that using the
CHAOS model to quantify the variability of the GVO series, provides only a
crude indication of the data errors, which are assumed time independent.
Moreover, correlation in the data errors are presently not accounted for,
although this is expected especially at high latitudes. It is however
challenging to empirically estimate non-diagonal data error covariance
matrices with the short GVO time series presently available. Further work is
needed on this important topic.

Table 2: Mean of the rms differences (in $\mathrm{pT/km\>yr^{-1}}$) between GVO SV series and GCV cubic spline fits for six of the gradient tensor elements. Results are shown for GVO SV gradient series derived from Swarm and CHAMP data using CHAOS-7.2 (Finlay et al., 2020) as MF model in the GVO processing. Component | $rms_{[rr]}$ | $rms_{[\theta\theta]}$ | $rms_{[\phi\phi]}$ | $rms_{[r\theta]}$ | $rms_{[r\phi]}$ | $rms_{[\theta\phi]}$
---|---|---|---|---|---|---
| $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$
CHAMP |  |  |  |  |  |
Polar | 5.40 | 3.90 | 4.20 | 2.10 | 2.60 | 2.60
non-Polar | 2.10 | 0.70 | 2.00 | 0.70 | 1.60 | 0.60
All | 2.98 | 1.55 | 2.55 | 1.03 | 1.83 | 1.14
Swarm |  |  |  |  |  |
Polar | 4.20 | 3.20 | 3.90 | 1.60 | 1.70 | 2.10
non-Polar | 1.20 | 0.30 | 1.20 | 0.30 | 0.70 | 0.60
All | 1.98 | 1.01 | 1.86 | 0.61 | 0.95 | 0.97

As with the ordinary GVO vector field time series, we estimate GVO gradient
tensor time series on a global grid of 300 GVOs. We compute the SV as annual
differences at each GVO for each tensor element. In order to quantify the
scatter level in each series, we then fit cubic smoothing splines to the time
series, with a knot spacing of 4 months and a smoothing parameter determined
using a generalized cross-validation (GCV) approach (Green & Silverman, 1993).
Table 2 presents the mean rms differences between the GVO SV gradient tensor
elements and GCV spline fits, separated into polar and non-polar regions.
These rms numbers provide an indication of the scatter level in the GVO SV
gradient data derived from the CHAMP and Swarm measurements first using
CHAOS-7.2 (Finlay et al., 2020) as a main field model. Comparing the numbers
between CHAMP and Swarm, we see that overall the values are lower for Swarm,
i.e. Swarm gradient tensor element SV time series have a lower scatter than
similar series for CHAMP. We note that the misfits are considerably lower at
non-Polar GVOs, compared with Polar GVOs, and this applies to all of the
gradient elements, during both CHAMP and Swarm times. Furthermore, the off-
diagonal elements and the $d\left[\nabla B\right]_{\theta\theta}/dt$ element
show considerably lower misfit values being below $2\>\mathrm{pT/km\>yr^{-1}}$
and $1\>\mathrm{pT/km\>yr^{-1}}$ at non-Polar latitudes during CHAMP and Swarm
times, respectively.

Table 3: Mean of the rms differences (in $\mathrm{pT/km\>yr^{-1}}$) between GVO SV series and GCV cubic spline fits for six of the gradient tensor elements. Results are shown for GVO SV gradient data derived from Swarm and CHAMP data using COV-OBS.x2 model (Huder et al., 2020) as MF model in the GVO processing. Component | $rms_{[rr]}$ | $rms_{[\theta\theta]}$ | $rms_{[\phi\phi]}$ | $rms_{[r\theta]}$ | $rms_{[r\phi]}$ | $rms_{[\theta\phi]}$
---|---|---|---|---|---|---
| $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$ | $[\mathrm{pT}/\mathrm{km}\>\mathrm{yr}^{-1}]$
CHAMP |  |  |  |  |  |
Polar | 5.40 | 3.90 | 4.20 | 2.10 | 2.60 | 2.60
non-Polar | 2.10 | 0.70 | 2.00 | 0.70 | 1.60 | 0.60
All | 2.99 | 1.55 | 2.55 | 1.03 | 1.83 | 1.14
Swarm |  |  |  |  |  |
Polar | 4.20 | 3.20 | 3.80 | 1.60 | 1.70 | 2.10
non-Polar | 1.20 | 0.30 | 1.20 | 0.30 | 0.70 | 0.60
All | 1.97 | 1.01 | 1.86 | 0.61 | 0.96 | 0.97

We also tested how the choice of main field model (used for subtracting and
adding back main field estimates) would impact the results. We produced test
GVO tensor element series from both CHAMP and Swarm measurements using the
main field values for SH degrees $n\in[1,13]$ of the COV-OBS.x2 model (Huder
et al., 2020). Table 3 presents the mean rms differences using the COV-OBS.x2
model instead of CHAOS-7.2. This results in almost identical misfit levels to
the GCV splines (i.e. scatter), between the GVO gradient series during CHAMP
and Swarm times, regardless of whether CHAOS-7.2 or COV-OBS.x2 is chosen.

##  4 Results

###  4.1 Field Gradient Element SV time series

We begin by investigating the temporal behaviour of the annual differences of
each gradient tensor element at an example GVO location above Honolulu ground
observatory in Hawaii, from which there are well known vector field records.
To do this, we compute dedicated GVO gradient element series above the
Honolulu observatory using the method described in Section 3.2. Here we are
motivated by studies which have pointed out a change in secular acceleration
of the radial component in the Pacific occurring around 2017 (Sabaka et al.,
2018; Finlay et al., 2020). In particular, we are interested to see if it is
possible to identify this event in the GVO gradient tensor time series, and
how this will display in the various tensor elements. Figure 2 present plots
of the SV for each gradient tensor element above Honolulu, showing the GVOs
derived from CHAMP (in blue) and Swarm (in red) measurements. For comparison
purposes we have mapped the two GVO series to a common altitude of 500 km by
subtracting the SV gradient field differences between the GVO altitudes and
500 km altitude using the CHAOS-7.2 model. The uncertainty estimates
$\pm\sigma$ are computed as for the global grid GVOs, i.e. the total mean
square error from residuals between the Honolulu GVO gradient tensor data and
gradient estimates from the CHAOS-7.2 model up to degree 16.

Visual inspection clearly demonstrates that geophysical signals are captured
in all of the SV gradient tensor elements. Distinctive changes centred around
2017 can be observed having a ”$V$” shape in the $d\left[\nabla
B\right]_{rr}/dt$ and $d\left[\nabla B\right]_{r\theta}/dt$ elements, with a
corresponding ”$\Lambda$” shape in the $d\left[\nabla
B\right]_{\theta\theta}/dt$ and $d\left[\nabla B\right]_{\phi\phi}/dt$
elements. In addition to this, we note that during 2004-2010, especially the
$d\left[\nabla B\right]_{r\theta}/dt$ element displays a variation pattern
which resembles that found in the $\theta$-component of the annual differences
of monthly mean vector field series from Honolulu (not shown).

Figure 2: Annual differences of the GVO field gradient elements with $\pm
1\sigma$ uncertainties, during CHAMP (blue) and Swarm (red) times at altitude
500km for a case study above Honolulu, Hawaii. Units are
$\mathrm{pT}/\mathrm{km}\,\mathrm{yr}^{-1}$.

Next, we investigate the global behaviour of annual differences of the
gradient elements for GVOs derived from Swarm measurements during 2014-2020.
Here we have chosen to present global series for the $d\left[\nabla
B\right]_{rr}/dt$ element in Figure 3. By visual inspection, we find local
regions with similar temporal changes as those observed at the Honolulu SV
gradient series. In particular, a distinct ”$V$” shaped behaviour is found in
the eastern Pacific region in a band stretching from latitudes $20^{\circ}$S
to $20^{\circ}$N and longitudes $180^{\circ}$ to $220^{\circ}$ with a possibly
related opposite ”$\Lambda$” shaped behaviour in the western Pacific region
from latitudes $20^{\circ}$S to $20^{\circ}$N and longitudes $120^{\circ}$E to
$180^{\circ}$E. These regional changes occur over a time window of 6 years
reaching amplitudes of about $15\>\mathrm{pT}/\mathrm{km}\,\mathrm{yr}^{-1}$.
Note that a ”$V$”-shaped SV gradient time series means a strong positive
change in the SA, while a ”$\Lambda$”-shaped time series means a strong
negative change in the SA. Though more complex to interpret, the other SV
gradient tensor elements (not shown) also exhibit distinctive behaviour in the
Pacific region. These observed changes in the SV gradient elements indicate
regional jerk-type event happening in the Pacific centred on 2017. In the
ionosphere external fields tend to be organized according to the geometry of
Earth’s main field, and their signal in the GVO series may therefore be
grouped accordingly to magnetic latitude in quasi-dipole coordinates (Laundal
& Richmond, 2017). Here magnetic latitude $\pm 70^{\circ}$ (dark blue curve)
may be used to approximate the border between North/South Polar and Auroral
zones, while magnetic latitude $\pm 50^{\circ}$ (light blue curve) divides the
North/South Auroral and Low- to Mid-latitude zones (Hammer et al., 2021a). In
all of the SV gradient element maps, higher scatter are found at GVOs located
in the Polar and Auroral zones, which is consistent with noise (unmodeled
fields) from ionospheric and magnetosphere-ionosphere coupling currents. In
addition, some rapid variations can be seen, which decrease in amplitude when
going towards equatorial latitudes. An example of this is found along
longitude $150^{\circ}$W, stretching from latitudes $30^{\circ}$S to
$60^{\circ}$S, as highlighted in the side-panels of Figure 3 at three selected
GVO locations (marked in the global maps by the green ellipse)

Figure 3: Time series of SV field gradient element $d\left[\nabla
B\right]_{rr}/dt$ (map) during Swarm time from 2014-2020 at 490km altitude.
Magnetic latitudes $\pm 50^{\circ}$ and $\pm 70^{\circ}$ shown with blue
curves. GVO locations marked with a red cross. Highlighted are selected time
series (locations marked with green ellipse) after removing the mean trend in
order to ease comparison, at three GVO locations along the $150^{\circ}$W
meridian line at latitudes $30^{\circ}$S (top), $42^{\circ}$S (center) and
$54^{\circ}$S (bottom).

Figure 4: Annual differences of the field gradient elements from GVO’s derived
using different data selection criteria, as described in the text, during
Swarm time for a case study above Honolulu, Hawaii. Units are
$\mathrm{pT}/\mathrm{km}\,\mathrm{yr}^{-1}$.

An important question is whether the prominent change in SV observed centred
on 2017 is robust and of internal origin. To address this, we produced a set
of the GVO SV series above Honolulu during Swarm time, testing a range of
geomagnetic selection criteria. We considered five cases: Case A using a dark
quiet time data selection removing estimates of the magnetospheric and
ionospheric fields as described in Section 2, this is our ’preferred’ criteria
for studying core field variations, Case B using a dark quiet time data
selection removing estimates of the magnetospheric field but not removing
estimates of the ionospheric field, Case C using a dark time data selection
but with neither ionospheric nor magnetospheric corrections, Case D using a
quiet time data selection from both day and night (”all” local times) and
removal of the magnetospheric field, and Case E using data from ”all” local
times, without any quiet time data selection applied and without removal of
for magnetospheric or ionospheric fields. Here ”dark” means the sun is
required to be at least $10^{\circ}$ below horizon, and ”all” means that no
such selection is applied, i.e. sunlit data are also included. SV gradient
series for all five data selection cases are shown in Figure 4. The black dots
corresponding to Case A are for our ”preferred” criteria were also used to
derive the maps in Figure 3. All selection criteria results in the same
overall temporal ”$\Lambda/V$” shape behaviour with an amplitude of $\approx
15\>\mathrm{pT/km\>yr^{-1}}$. There is no increase in amplitude when including
more disturbed data. For example, comparing Case C (purple star) with Case
A/Case B (black/blue dots) should expose a signal from a magnetospheric
source; however, the same ”$\Lambda/V$” shape behaviour appears in all three
cases. Comparing instead Case E (yellow star) with Case C and Case D should
expose an ionospheric signal, which is expected to be larger during sunlit
conditions; even though more scatter is seen in Case E, the same overall
”$\Lambda/V$” shape is clearly visible and of similar amplitude. These results
are consistent with an internal origin of the 2017 SV impulse event.

###  4.2 Example Spherical Harmonic Models Derived From Gradient Data

In this section we demonstrate that spherical harmonic (SH) field models with
high temporal resolution (4 months) can be built from the global network of
GVO gradient tensor time series. We then use these models to investigate
global changes in SA during Swarm time, and in particular, we analyse the
possible benefits of the GVO field gradient tensor series over more standard
GVO vector field series. The magnetic field can be described as the negative
gradient of internal, $V^{int}$, and external, $V^{ext}$, potentials (e.g.
Sabaka et al., 2010)

$\mathbf{B}=-\nabla V^{int}-\nabla V^{ext}\;.$ (23)

Here the internal and external potentials at a given epoch can be expressed by
truncated spherical harmonic expansions up to degree $N^{int}$ and $N^{ext}$,
respectively

$\displaystyle V^{int}(r,\theta,\phi,t)$
$\displaystyle=r_{a}\sum_{n=1}^{N^{int}}\sum_{m=0}^{n}\left[g_{n}^{m}(t)\mathrm{cos}(m\phi)+h_{n}^{m}(t)\mathrm{sin}(m\phi)\right]\left(\frac{r_{a}}{r}\right)^{n+1}P_{n}^{m}(\mathrm{cos}\,\theta)$
(24) $\displaystyle V^{ext}(r,\theta,\phi,t)$
$\displaystyle=r_{a}\sum_{n=1}^{N^{ext}}\sum_{m=0}^{n}\left[q_{n}^{m}(t)\mathrm{cos}(m\phi)+s_{n}^{m}(t)\mathrm{sin}(m\phi)\right]\left(\frac{r}{r_{a}}\right)^{n}P_{n}^{m}(\mathrm{cos}\,\theta)\;,$
(25)

where $r_{a}=6371.2$km is the Earth’s mean spherical radius, $n$ and $m$ are
the SH degree and order, respectively, and $P_{n}^{m}$ are the associated
Schmidt semi-normalized Legendre functions. $\\{g_{n}^{m},h_{n}^{m}\\}$ and
$\\{q_{n}^{m},s_{n}^{m}\\}$ are the internal and external expansion
coefficients, respectively. Assuming an internal field source, the linear
forward problem of determining the SH expansion coefficients can be written

$\mathbf{d}=\underline{\underline{\mathbf{G}}}\mathbf{m},$ (26)

where $\mathbf{d}$ is a data vector containing the GVO epoch data (i.e. field
vector components or field gradient tensor elements),
$\underline{\underline{\mathbf{G}}}$ is the design matrix for an internal
potential relating each model coefficient to either the GVO epoch vector or
gradient data (for explicit expressions of the gradient tensor elements, we
refer to Appendix A of Kotsiaros (2012)), and vector $\mathbf{m}$ contains the
parameters of the potential, i.e. the internal SH coefficients here denoted as
$g_{n}^{m}$ and $h_{n}^{m}$ for order $m$ and degree $n$. For each GVO epoch
we estimate a single epoch SH model using an Iteratively Reweighted Least
Squares method which is applied in order to mitigate the impact of non-
Gaussian distributed data residuals (e.g. Olsen et al., 2006)

$\mathbf{m}=(\underline{\underline{\mathbf{G}}}^{T}\underline{\underline{\mathbf{W}}}\
\underline{\underline{\mathbf{G}}})^{-1}\underline{\underline{\mathbf{G}}}^{T}\
\underline{\underline{\mathbf{W}}}\mathbf{d}\;,$ (27)

where we assign weights according to the diagonal weight matrix
$\underline{\underline{\mathbf{W}}}=w/\sigma^{2}$ consisting of Huber weights,
$w$, having a Huber tuning constant of 1.5 (e.g., Constable, 1988) and error
estimates, $\sigma^{2}$, of either the field vector components or field
gradient tensor elements, depending on which type of data is solved for. These
error estimates are derived as the total mean squared error (Bendat & Piersol,
2010) of the residuals between the GVO vector/gradient tensor series and
vector/gradient estimates of the CHAOS-7.2 model up to degree 16 (see Section
3.2). That is, the weight matrix consists of the inverse data covariance
matrix elements multiplied by Huber weights applied in order to account for
non-Gaussian data residuals. The weights are then iteratively updated until
convergence. We derive models up to SH degree $N^{int}=14$ for each 4 month
interval.

To investigate the 2017 region jerk event we next compute the secular
acceleration change for each gradient tensor element between 2015.5 and 2018.5
at the Earth’s surface

$\Delta d^{2}\left[\nabla B\right]_{jk}/dt^{2}=d^{2}\left[\nabla
B\right]_{jk}/dt^{2}|_{2018.5}-d^{2}\left[\nabla
B\right]_{jk}/dt^{2}|_{2015.5}\;.$ (28)

Plotting global maps of this change in Figure 5 for each of the gradient
elements for degrees $n\leq 9$ at the Earth’s surface, distinct patterns of SA
change are seen to have occurred in the Pacific region during 2015.5-2018.5.
Only results for the upper right part of the gradient tensor are shown as the
tensor is symmetric. The $\Delta d^{2}\left[\nabla B\right]_{rr}/dt^{2}$ map
identifies two strong localized patches of opposite sign in SA change reaching
amplitudes of $40\>\mathrm{pT/km\,yr^{-2}}$ in a region confined to latitudes
$25^{\circ}$S to $25^{\circ}$N and longitudes $140^{\circ}$ to $220^{\circ}$.
Associated strong negative and positive patches are seen in the $\Delta
d^{2}\left[\nabla B\right]_{\phi\phi}/dt^{2}$ map in the same region. In
addition, the $\Delta d^{2}\left[\nabla B\right]_{r\theta}/dt^{2}$ and $\Delta
d^{2}\left[\nabla B\right]_{\theta\phi}/dt^{2}$ elements show a tiling pattern
of positive and negative field patches from latitudes $25^{\circ}$S to
$25^{\circ}$N and longitudes $120^{\circ}$ to $240^{\circ}$. Similar changes,
but in the radial field SA between 2014 to 2020, involving nearby features
with opposite sign in the Pacific region, have been found in the CHAOS-7.2
field model Finlay et al. (2020) and using the technique of Subtractive
Optimized Local Averages (SOLA) applied to Swarm data Hammer & Finlay (2019);
Hammer et al. (2021c).

Figure 5: Change in SA gradient tensor elements between 2015.5 and 2018.5 for
SH degrees $n\leq 9$ at the Earth’s surface.

Figure 6: MF (dashed-dotted), SV (solid) and SA (dashed) CMB mean power
spectra of epoch models truncated at SH degree $n=14$. Left plot show spectra
from vector field (blue) and field gradient (red) GVOs, derived using all
Swarm satellites. Right plot show curves based on gradient GVOs derived using:
all Swarm satellites (red, SAT-ABC), Swarm-Alpha (dark blue, SAT-A), Swarm-
Alpha and Bravo (green, SAT-AB) and Swarm-Alpha and Bravo (light blue, SAT-
AC).

Next, we seek to further inspect and compare the SH models obtained from GVO
gradient series with similar models obtained using more traditional GVO vector
component series. The left plot of Figure 6 presents the mean of the MF
(dashed-dotted), SV (solid) and SA (dashed) spatial Lowes-Mauersberger power
spectra (Hulot et al., 2015), at the Core-mantle boundary (CMB) obtained by
computing the time-averaged spectrum taken over all the epoch spectral lines.
Spectral curves in blue and red are derived from GVO gradients and GVO vector
series, respectively. The SV and SA power spectra derived from GVO gradient
tensor series are seen to diverge less rapidly as compared to models derived
from GVO vector series, and the SV/SA intersection happens at a slightly
higher degree (11 compared to 9). This behaviour is consistent with the
analyzes of Kotsiaros & Olsen (2014), who found that gradient data better
constrain SV to higher SH degrees than vector data.

In order to characterize the influence of satellite configuration on the GVO
gradient series, we consider case studies consisting of: 1) a single
satellite, 2) two satellites at different altitudes which would provide a
better local-time coverage as compared to the first case, 3) East-West
capability from side-by-side flying satellites. For each of these cases, we
derived GVO gradient series using: 1) Swarm Alpha only (SAT-A), 2) both Swarm
Alpha and Bravo (SAT-AB) covering altitude ranges of approx. 430-515km and
490-530km during 2014-2020, respectively, 3) Swarm Alpha and Charlie (SAT-AC)
taking advantage of their East-West longitudinal separation of approx.
$1.4^{\circ}$ (corresponding to 155km at the equator) (Olsen et al., 2015).
The right plot of Figure 6 presents the time-averaged MF (dashed-dotted), SV
(solid) and SA (dashed) spatial power spectra, derived from each of these
three case. Note here, that all of the MF spectra overlap and are thus hidden
by the red SAT-ABC spectrum. In addition, note that this red SAT-ABC spectrum
in the right panel refers to the same model as the red MF-GVO Gradient
spectrum in the left panel of Figure 6. We discuss these results and their
implications in more detail in the next section.

Investigating further these SH models, Figure 7 shows the first time
derivative of the internal expansion coefficients, computed based on simple
first differences, derived from internal models using the GVO vector (blue,
and also including an external expansion, in green) and GVO gradient (red)
series derived using all three Swarm satellites. Example coefficients are
shown for zonal, $m=0$, terms (top row), tesseral, $m\neq n$, terms (middle
row) and sectorial, $m=n$, terms (bottom row). To quantify the scatter level
in the epoch coefficient series, standard deviations between the coefficient
series and GCV smoothing spline fits (solid curves) are given in each case.
Although robust estimation has been applied when deriving these models,
outliers can be seen in the both series. A change in the sign of the trend in
the SV signal is evident, especially in the sectorial coefficients
$\dot{h}_{3}^{3},\dot{h}_{4}^{4},\dot{h}_{6}^{6}$, but also in the
$\dot{h}_{5}^{3}$ coefficient centered on 2017. The scatter level in the near
zonal coefficients of the GVO vector series are reduced if an external SH
expansion (in green dots and curve) is included as well. This points to the
presence of external contamination in the vector series, which seems to be
reduced in the gradient series.

Figure 8 collects such standard deviations for SH degrees and order up to 12,
from models derived using the GVO vector (left plot) and GVO gradient (right
plot) series. We generally find less scatter in the GVO gradient series, and
especially in the zonal and near-zonal (where $m$ is close to zero)
coefficients. For the sectorial terms the scatter levels are low for both
series. Use of the vector series results in higher scatter levels for the
near-zonal terms for degrees $n>2$ and orders $m\leq 2$. When also including
an external SH expansion, i.e. using both internal and external expansion
terms up to degree 14, for the GVO vector series (middle plot), we are able to
reduce the scatter level in the near-zonal coefficients (middle plot of Figure
8), illustrating that using the GVO gradient elements in SH modelling helps in
excluding external field signals. Note here that we focus on comparing SH
models derived from GVO vector data with those derived solely using GVO
gradient data (including an external SH expansion for the GVO gradient series
would require vector information as well to obtaining a robust estimation).
Global maps (not shown) show that much of the enhanced scatter is related to
signals in polar regions being spuriously mapped into the internal field.

Finally, we show in Figure 9, the standard deviations from SH models based on
the case study GVO gradient series, which were derived using Swarm Alpha only
(case 1, left plot), Swarm Alpha and Bravo (case 2, center plot), and Swarm
Alpha and Charlie (case 3, right plot).

Figure 7: Series of first time derivatives for some internal coefficients
$dg_{n}^{m}/dt$ and $dh_{n}^{m}/dt$ derived from GVO vector (blue dots), GVO
vector including an external SH expansion (green dots) and GVO gradient (red
dots) data. Standard deviations of differences between the series and a GCV
smoothing spline fit (solid curves) to the coefficients are given. Units are
nT/yr.

Figure 8: Standard deviations of differences between the first time derivative
of internal SH model coefficient series and spline-fitted curves for each
series derived from GVO vector data (left plot) and GVO gradient data (right
plot), and GVO vector data including an external SH expansion (middle plot).
Positive orders $m$ refer to the coefficients $dg_{n}^{m}/dt$, while negative
orders refer to $dh_{n}^{m}/dt$ coefficients. Units are nT/yr.

Figure 9: Standard deviations of differences between the first time derivative
of internal SH model coefficient series and spline-fitted curves for each
series based on GVO gradient series which has been derived using: Swarm-Alpha
data (left plot, SAT-A), Swarm Alpha and Bravo data (center plot SAT-AB), and
Swarm Alpha and Charlie data (right, SAT-AC). Positive orders $m$ refer to the
coefficients $dg_{n}^{m}/dt$, while negative orders refer to $dh_{n}^{m}/dt$
coefficients. Units are nT/yr.

##  5 Discussion and Conclusions

In this study we have extended the existing GVO concept and derived time
series of the second-order gradient tensor elements of the geomagnetic field
at a global network of 300 locations. We have computed such GVO gradient time
series from the mean and differences of vector magnetic field measurements,
along track and in the east-west direction, from the low Earth orbiting CHAMP
and Swarm satellites.

Inspecting SV gradient tensor elements for a GVO located above the Honolulu
ground observatory we found evidence in the gradient series for a regional
jerk-type event centered on 2017, observed as a characteristic ”$V$” shaped
change in the $d\left[\nabla B\right]_{rr}/dt$ and $d\left[\nabla
B\right]_{r\theta}/dt$ elements, and as a ”$\Lambda$” shape in the
$d\left[\nabla B\right]_{\theta\theta}/dt$ and $d\left[\nabla
B\right]_{\phi\phi}/dt$ elements. In the global GVO SV gradient element
records, spanning the years from 2014 to 2020, we find evidence of robust time
variations in many of the tensor elements. In particular, intense fluctuations
in the Pacific region confined in longitude, suggest a regionally localized
geomagnetic impulse event taking place around 2017. This is consistent with
ground observatory measurements of the SV of the radial magnetic field
component at the Honolulu observatory (e.g. Finlay et al., 2020; Sabaka et
al., 2018). By changing the geomagnetic quiet-time and local time selection
criteria, we see little change in the amplitude of this jerk signal,
supporting the hypothesis that the 2017 event is of internal origin. At the
Earth’s surface nearby patches of intense change in the SA gradient field,
with opposite signs, occurs between 2015.5 and 2018.5. These are found to be
limited to latitudes between $25^{\circ}$S to $25^{\circ}$N and to longitudes
between $140^{\circ}$ to $220^{\circ}$E. In particular, two strong patches of
change in the radial gradient of the radial field, with opposite signs, locate
the centre of the 2017 jerk event to approximately $0^{\circ}$N and
$170^{\circ}$E in the central Pacific. We also find (not shown) that we can
robustly map the change in SA at the CMB up to degree 7 using the 4-monthly
gradient tensor element data. Although maps at the CMB exhibit enhanced noise
due to downward continuation of the field they do display coherency regarding
the distinctive SA changes observed in the Pacific region which are clearly
related to those found at the Earth’s surface in Figure 5, further supporting
an internal origin of the 2017 SA impulse.

Various geophysical explanations of geomagnetic jerk events, similar to those
we have highlighted here in the Pacific region, have been proposed. The
possibilities still under discussion include equatorially trapped MAC waves in
a possible stratified layer close to the core surface (Buffett & Matsui, 2019;
Chi-Durán et al., 2020) and equatorial focusing of hydrodynamic waves
originating from turbulent convection deep within the core (Aubert & Finlay,
2019; Gerick et al., 2021). In that connection the new concept of GVO gradient
tensor time series may aid future studies of the appearance and dynamics of
geomagnetic jerks, related changes in core flows and core dynamics via e.g.
data assimilation. In a follow-up study with K. Whaler (in prep), we present
computations and investigations of core surface flows derived from GVO
gradient tensor elements series, paying particular attention to the jerk in
the Pacific region in 2017.

In addition to the signature of geomagnetic jerks, global maps of the GVO SV
tensor element series, also show some example rapid smaller amplitude SV
fluctuations within a few years, especially noticeable in the $d\left[\nabla
B\right]_{rr}/dt$ element (see left panel of Figure 3). Some of these small
amplitude variations in the SV gradient series might be indicative of external
field leakage. The most obvious sources of external field leakage are 1) rapid
variations in polar latitude GVO series that persist to lower latitudes along
the same meridional, which could be an indication of contamination by polar
electrojet or field aligned currents, 2) rapid variations at mid or low
latitude GVO series seen at all longitudes that could be caused by a signal
having magnetospheric origin. We find no signs of distinct temporal variations
along all longitudes at mid/low latitudes, suggesting that remaining
magnetospheric disturbances are likely small. We do find some examples of
rapid variations that seem to diminish towards lower latitudes, which could
signify contamination from polar electrojet or field aligned currents systems.
Longer time series are needed in order to study the origin of such variability
in the GVO gradient series.

In order to test for possible improvements in retrieving the SV signal using
GVO gradient tensor data as compared to GVO vector data, we produced simple
unregularized SH field models built from the GVO gradient and vector data
derived using Swarm measurements. Comparing the power spectra of these models
shown in Figure 6 supports the findings of Kotsiaros & Olsen (2014), that
harmonics of the SV above degree 6 can be better resolved when using gradient
tensor data than using vector data. In particular, analysis of the first time
derivatives of the SH coefficients as presented in Figure 8, shows that
especially zonal and near-zonal harmonics of models derived from GVO gradients
have less scatter compared to similar models derived from GVO vector data.

We find that the scatter level in the SV gradient tensor elements is higher
during CHAMP times compared with Swarm times. The orbital configuration of the
three Swarm satellites is advantageous for computing the gradient tensor as
more data are available and Swarm Bravo, having a slightly higher altitude
than Alpha and Charlie, provides information on the radial gradient which
enables better potential determination. The right panel of Figure 6 shows that
having two satellites, one of which is at a higher orbit, slightly lowers the
spatial power (green) as compared with having one satellite (dark blue). In
addition, from Figure 9, having two satellites lowers the scatter of the
tesseral coefficients for $m>2$. Considering instead two satellites in a side-
by-side orbit, the associated spectral line (light blue in the right panel of
Figure 6) is only marginally lower, but differences are most obvious from
Figure 9, where the scatter level of the tesseral coefficients is smaller for
orders $m>2$. The zonal terms for $n<7$ are found to have higher scatter, but
this is expected as East-West information alone does not allow determination
of zonal terms (Kotsiaros & Olsen, 2012). Comparing the three case studies
with the original case of using all three Swarm satellites, the associated
spatial power (red in the right panel of Figure 6) is lower, while the scatter
level across all SH coefficients as shown in Figure 8 is the lowest, except
for lowest zonal terms which seems to be related to including the East-West
gradient information as shown in the right plot of Figure 9. In order to
mitigate such behavior of the low degree zonal terms, in future studies it
might be worth exploring a Selective Infinite-Variance Weighting (SIVW)
approach (Kotsiaros & Olsen, 2014; Sabaka et al., 2013), wherein weights are
applied to data subsets that are more sensitive to certain parameter subsets.
By this approach, the East-West data could in the future be given less weight
with regard to parts of the GVO potential that provide information of the
zonal terms. The larger misfits during CHAMP times, may also result from a
closer proximity of the lower flying CHAMP satellite to ionospheric current
systems. From our tests we therefore conclude that the Swarm satellite trio is
advantageous for deriving GVO gradient series. Having satellites at different
altitudes better fills the 3D space of the GVO data cylinder, improving the
recovery of the gradient quantities, for example improving the determination
of the tesseral harmonics. Using all three satellites, enhances the recovery
of all harmonic coefficients, and is clearly superior to having a single
satellite.

## Availability of datasets and material

The GVO gradient tensor data underlying this article and their associated
uncertainty estimates are available from https://data.dtu.dk/, at (Hammer et
al., 2021b). The datasets used in this article are available in the following
repositories: Swarm data are available from
https://earth.esa.int/web/guest/swarm/data-access; CHAMP data are available
from https://isdc.gfz-potsdam.de/champ-isdc; Ground observatory data are
available from ftp://ftp.nerc-murchison.ac.uk/geomag/Swarm/AUX_OBS/hour/; The
RC-index is available from http://www.spacecenter.dk/files/magnetic-
models/RC/; The CHAOS-7.2 model and its updates are available at
http://www.spacecenter.dk/files/magnetic-models/CHAOS-7/; solar wind speed,
interplanetary magnetic field, and Kp-index are available from
https://omniweb.gsfc.nasa.gov/ow.html.

###### Acknowledgements.

We wish to thank Gauthier Hulot and an anonymous reviewer for constructive
comments that helped us improve the manuscript. We thank the GFZ German
Research Centre for Geoscience for providing access to the CHAMP MAG-L3 data
and the European Space Agency (ESA) for providing access to the Swarm L1b
data. The high resolution 1-min OMNI data were provided by the Space Physics
Data Facility (SPDF), NASA Goddard Space Flight Centre. We thank the staff of
the geomagnetic observatories and the INTERMAGNET for providing high-quality
observatory data. This project was partly supported by the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 772561). In addition, this study was also
partly funded by ESA through the Swarm DISC activities, contract no.
4000109587.

## References

  * Aubert & Finlay (2019) Aubert, J. & Finlay, C. C., 2019.  Geomagnetic jerks and rapid hydromagnetic waves focusing at Earth’s core surface, Nature Geoscience, 12(5), 393–398.
  * Backus et al. (1996) Backus, G., Parker, R., & Constable, C., 1996.  Foundations of Geomagnetism, Cambridge Univ. Press, New York.
  * Barrois et al. (2018) Barrois, O., Hammer, M. D., Finlay, C. C., Martin, Y., & Gillet, N., 2018.  Assimilation of ground and satellite magnetic measurements: inference of core surface magnetic and velocity field changes, Geophys. J. Int., 215, 695–712.
  * Beggan et al. (2009) Beggan, C. D., Whaler, K. A., & Macmillan, S., 2009.  Biased residuals of core flow models from satellite-derived virtual observatories, Geophys. J. Int., 177(2), 463–475.
  * Bendat & Piersol (2010) Bendat, J. & Piersol, A., 2010.  Random Data, Analysis and Measurement Procedures, Wiley, New Jersey.
  * Buffett & Matsui (2019) Buffett, B. & Matsui, H., 2019.  Equatorially trapped waves in Earth’s core, Geophys. J. Int., 218(2), 1210–1225.
  * Casotto & Fantino (2009) Casotto, S. & Fantino, E., 2009.  Gravitational gradients by tensor analysis with application to spherical coordinates, Journal of Geodesy, 83(7), 621–634.
  * Chi-Durán et al. (2020) Chi-Durán, R., Avery, M. S., Knezek, N., & Buffett, B. A., 2020.  Decomposition of geomagnetic secular acceleration into traveling waves using complex empirical orthogonal functions, Geophys. Res. Lett., p. e2020GL087940.
  * Constable (1988) Constable, C. G., 1988.  Parameter estimation in non-Gaussian noise, Geophys. J. Int., 94(1), 131–142.
  * Domingos et al. (2019) Domingos, J., Pais, M. A., Jault, D., & Mandea, M., 2019.  Temporal resolution of internal magnetic field modes from satellite data, Earth, Planets and Space, 71(1), 1–17.
  * Finlay (2019) Finlay, C. C., 2019.  Models of the main geomagnetic field based on multi-satellite magnetic data, Ionospheric Multi-Spacecraft Analysis Tools: Approaches for Deriving Ionospheric Parameters, 17, 255.
  * Finlay et al. (2017) Finlay, C. C., Lesur, V., Thébault, E., Vervelidou, F., Morschhauser, A., & Shore, R., 2017.  Challenges handling magnetospheric and ionospheric signals in internal geomagnetic field modelling, Space Science Reviews, 206(1-4), 157–189.
  * Finlay et al. (2020) Finlay, C. C., Kloss, C., Olsen, N., Hammer, M. D., Tøffner-Clausen, L., Grayver, A., & Kuvshinov, A., 2020.  The CHAOS-7 geomagnetic field model and observed changes in the South Atlantic Anomaly, Earth, Planets and Space, 72(1), 1–31.
  * Gerick et al. (2021) Gerick, F., Jault, D., & Noir, J., 2021.  Fast Quasi-Geostrophic Magneto-Coriolis Modes in the Earth’s core, Geophys. Res. Lett., 48(4), e2020GL090803, doi:10.1029/2020GL090803.
  * Green & Silverman (1993) Green, P. J. & Silverman, B. W., 1993.  Nonparametric regression and generalized linear models: a roughness penalty approach, Chapman and Hall.
  * Hammer (2018) Hammer, M. D., 2018.  Local Estimation of the Earth’s Core Magnetic Field, Ph.D. thesis, Technical University of Denmark.
  * Hammer & Finlay (2019) Hammer, M. D. & Finlay, C. C., 2019.  Local averages of the core–mantle boundary magnetic field from satellite observations, Geophys. J. Int., 216(3), 1901–1918.
  * Hammer et al. (2021a) Hammer, M. D., Cox, G., Brown, W., Beggan, C. D., & Finlay, C. C., 2021a.  Geomagnetic Virtual Observatories: monitoring geomagnetic secular variation with the Swarm satellites, Earth, Planets and Space, 73(1), 1–22.
  * Hammer et al. (2021b) Hammer, M. D., Finlay, C. C., & Olsen, N., 2021b.  Secular variation signals in magnetic field gradient tensor elements derived from satellite-based geomagnetic virtual observatories., Technical University of Denmark. Dataset. doi:10.11583/DTU.14695590.
  * Hammer et al. (2021c) Hammer, M. D., Finlay, C. C., & Olsen, N., 2021c.  Applications of CryoSat-2 satellite magnetic data in studies of the Earth’s core field variations, Earth, Planets and Space, 73(1), 1–22.
  * Huder et al. (2020) Huder, L., Gillet, N., Finlay, C. C., Hammer, M. D., & Tchoungui, H., 2020.  COV-OBS.x2: 180 years of geomagnetic field evolution from ground-based and satellite observations, Earth, Planets and Space, 72(1), 1–18.
  * Hulot et al. (2015) Hulot, G., Sabaka, T. J., Olsen, N., & Fournier, A., 2015.  The present field, in Treatise on Geophysics, vol. 5, chapter 5.02, ed. Kono, M., Elsevier Ltd.
  * Kloss & Finlay (2019) Kloss, C. & Finlay, C. C., 2019.  Time-dependent low latitude core flow and geomagnetic field acceleration pulses, Geophys. J. Int., 217, 140–168.
  * Koop (1993) Koop, R., 1993.  Global gravity field modelling using satellite gravity gradiometry, Netherlands Geodetic Comission, Publications on Geodesy, New Series, Number 38, Delft, The Netherlands.
  * Kotsiaros (2012) Kotsiaros, S., 2012.  Determination of Earth’s magnetic field from satellite constellation magnetic field observations, Ph.D. thesis, DTU Space.
  * Kotsiaros (2016) Kotsiaros, S., 2016.  Toward more complete magnetic gradiometry with the Swarm mission, Earth, Planets and Space, 68(1), 1–13.
  * Kotsiaros & Olsen (2012) Kotsiaros, S. & Olsen, N., 2012.  The geomagnetic field gradient tensor, GEM-International Journal on Geomathematics, 3(2), 297–314.
  * Kotsiaros & Olsen (2014) Kotsiaros, S. & Olsen, N., 2014.  End-to-end simulation study of a full magnetic gradiometry mission, Geophys. J. Int., 196(1), 100–110.
  * Kotsiaros et al. (2015) Kotsiaros, S., Finlay, C., & Olsen, N., 2015.  Use of along-track magnetic field differences in lithospheric field modelling, Geophys. J. Int., 200(2), 878–887.
  * Laundal & Richmond (2017) Laundal, K. M. & Richmond, A. D., 2017.  Magnetic coordinate systems, Space Science Reviews, 206(1-4), 27–59.
  * Leopardi (2006) Leopardi, P., 2006.  A partition of the unit sphere into regions of equal area and small diameter, Electronic Transactions on Numerical Analysis, 25(12), 309–327.
  * Mandea & Olsen (2006) Mandea, M. & Olsen, N., 2006.  A new approach to directly determine the secular variation from magnetic satellite observations, Geophys. Res. Lett., 33(15).
  * Nogueira et al. (2015) Nogueira, T., Scharnagl, J., Kotsiaros, S., & Schilling, K., 2015.  NetSat-4G A four nano-satellite formation for global geomagnetic gradiometry, in Proceedings of 10th IAA Symposium on Small Satellites for Earth Observation.
  * Olsen & Kotsiaros (2011) Olsen, N. & Kotsiaros, S., 2011.  Magnetic satellite missions and data, in Geomagnetic Observations and Models, pp. 27–44, Springer.
  * Olsen & Mandea (2007) Olsen, N. & Mandea, M., 2007.  Investigation of a secular variation impulse using satellite data: The 2003 geomagnetic jerk, Earth Planet. Sci. Lett., 255(1), 94–105.
  * Olsen & Stolle (2012) Olsen, N. & Stolle, C., 2012.  Satellite geomagnetism, Annu. Rev. Earth Planet. Sci., 40, 441–465.
  * Olsen et al. (2006) Olsen, N., Lühr, H., Sabaka, T. J., Mandea, M., Rother, M., Tøffner-Clausen, L., & Choi, S., 2006.  CHAOS—a model of the Earth’s magnetic field derived from CHAMP, Ørsted, and SAC-C magnetic satellite data, Geophys. J. Int., 166(1), 67–75.
  * Olsen et al. (2014) Olsen, N., Lühr, H., Finlay, C. C., Sabaka, T. J., Michaelis, I., Rauberg, J., & Tøffner-Clausen, L., 2014.  The CHAOS-4 geomagnetic field model, Geophys. J. Int., 197(2), 815–827.
  * Olsen et al. (2015) Olsen, N., Hulot, G., Lesur, V., Finlay, C. C., Beggan, C., Chulliat, A., Sabaka, T. J., Floberghagen, R., Friis-Christensen, E., Haagmans, R., et al., 2015\.  The Swarm Initial Field Model for the 2014 geomagnetic field, Geophys. Res. Lett., 42(4), 1092–1098.
  * Olsen et al. (2017) Olsen, N., Ravat, D., Finlay, C. C., & Kother, L. K., 2017.  LCS-1: a high-resolution global model of the lithospheric magnetic field derived from CHAMP and Swarm satellite observations, Geophys. J. Int., 211(3), 1461–1477.
  * Reed (1973) Reed, G., 1973.  Application of kinematical geodesy for determining the shorts wavelength component of the gravity field by satellite gradiometry, Ph.D. thesis, The Ohio State University, Dept. of Geod Science, Rep. No. 201, Columbus, Ohio.
  * Riley et al. (2004) Riley, K. F., Hobson, M., & Bence, S., 2004.  Mathematical Methods for Physics and Engineering, Cambridge University Press, Cambridge.
  * Ritter et al. (2004) Ritter, P., Lühr, H., Maus, S., & Viljanen, A., 2004.  High-latitude ionospheric currents during very quiet times: their characteristics and predictability, Annales Geophysicae, 22, 2001–2014.
  * Rogers et al. (2019) Rogers, H. F., Beggan, C. D., & Whaler, K. A., 2019.  Investigation of regional variation in core flow models using spherical slepian functions, Earth, Planets and Space, 71(1), 19\.
  * Sabaka et al. (2010) Sabaka, T. J., Hulot, G., & Olsen, N., 2010.  Mathematical properties relevant to geomagnetic field modeling, in Handbook of Geomathematics, pp. 503–538, Springer.
  * Sabaka et al. (2013) Sabaka, T. J., Tøffner-Clausen, L., & Olsen, N., 2013.  Use of the Comprehensive Inversion method for Swarm satellite data analysis, Earth, Planets and Space, 65(11), 1201–1222.
  * Sabaka et al. (2018) Sabaka, T. J., Tøffner-Clausen, L., Olsen, N., & Finlay, C. C., 2018.  A Comprehensive Model of the Earth’s Magnetic Field Determined From 4 Years of Swarm Satellite Observations, Earth, Planets and Space, 70(1), 1–26.
  * Shore (2013) Shore, R. M., 2013.  An improved description of Earth’s external magnetic fields and their source regions using satellite data, Ph.D. thesis, The University of Edinburgh.
  * Tscherning (1976) Tscherning, C., 1976.  Computation of the second-order derivatives of the normal potential based on the representation by a Legendre series, Manuscripta geodaetica, 1, 71–92.

##  Appendix A The Magnetic Field Gradient Tensor and General Coordinate
Transformations

Here we provide details on the magnetic gradient tensor and how its elements
transform between different coordinate systems. In particular, we are
interested in the transformation between the local topocentric Cartesian
coordinate system described in Section 3.1 and the spherical polar coordinate
system. Formulations from gravimetry of the gravitational gradient tensor
(also referred to as the Marussi tensor) can be found in Reed (1973); Koop
(1993); Casotto & Fantino (2009); Tscherning (1976). Here we follow the
notation of Casotto & Fantino (2009), which is inspired by common usage in
general relativity. The reader should however take care concerning the
differences between the magnetic and gravity cases, and in particular, of the
coordinate systems adopted, i.e. their orientation and whether they are left-
or right-handed systems.

Referring to a point $P$ (which would denote a given GVO target point), the
usual geocentric system is given by the Cartesian coordinates as
$\tilde{x}^{p}=(\tilde{x},\tilde{y},\tilde{z})$ and by the spherical polar
coordinates as $(r,\theta,\phi)$, where $\theta$ is the colatitude. The
geocentric system can be described by the Cartesian unit vectors
$(\mathbf{\hat{i}_{1}},\mathbf{\hat{i}_{2}},\mathbf{\hat{i}_{3}})$ denoting
the basis $i_{p}$. At $P$ a local Cartesian coordinate system $(z,x,y)$ is
defined by the basis $\mathbf{e}_{p}$ where $p=1,2,3$, which is same one as
used in GVO method, see Section 3.1. This covariant right-handed orthogonal
basis is determined by the components of the partial derivatives of the
position vector $\mathbf{r}$ as: $\mathbf{e}_{1}=\partial\mathbf{r}/\partial
r$ pointing radially outwards,
$\mathbf{e}_{2}=\partial\mathbf{r}/\partial\theta$ pointing to the south and
$\mathbf{e}_{3}=\partial\mathbf{r}/\partial\phi$ pointing to the east, i.e.
similar to the spherical polar basis vector at the target point $P$. The
normalized basis vectors for the local system are given as
$\hat{\mathbf{e}}_{1}=\mathbf{e}_{1}$,
$\hat{\mathbf{e}}_{2}=(1/r)\mathbf{e}_{2}$,
$\hat{\mathbf{e}}_{3}=1/(r\mathrm{sin}\,\theta)\mathbf{e}_{1}$, such that the
position vector can be written as
$\mathbf{r}=z\hat{\mathbf{e}_{1}}+x\hat{\mathbf{e}_{2}}+y\hat{\mathbf{e}_{3}}$.
Notice that while the basis vectors $\mathbf{i}_{p}$ are constant in magnitude
and direction, the basis vectors $\hat{\mathbf{e}}_{p}$ have constant
magnitude but their directions vary (the same goes for the spherical basis
vectors). Thus when computing the spatial derivatives of a vector, the basis
vectors also need to be differentiated as they depend on position. The
position vector from origin $O$ to the point $P$ can be written in geocentric
Cartesian coordinates as (Riley et al., 2004)

$\mathbf{r}=\tilde{x}\mathbf{\hat{i}_{1}}+\tilde{y}\mathbf{\hat{i}_{2}}+\tilde{z}\mathbf{\hat{i}_{3}}=\tilde{x}^{p}i_{p}\;,$
(29)

where the summation convention has been used. The Cartesian coordinates are
related to the spherical coordinates via (Riley et al., 2004, p. 363)

$\mathbf{r}=\begin{bmatrix}\tilde{x}\\\ \tilde{y}\\\
\tilde{z}\end{bmatrix}=r\begin{bmatrix}\mathrm{sin}\theta\mathrm{cos}\phi\\\
\mathrm{sin}\theta\mathrm{sin}\phi\\\ \mathrm{cos}\theta\end{bmatrix}\;.$ (30)

The magnetic scalar potential, $V$, can be considered as a tensor of zero-
order (or rank). The gradient operator in the generalized coordinates $u^{p}$,
where $p=1,2,3$, having the covariant basis
$\mathbf{e}_{p}=\partial\mathbf{r}/\partial u^{p}$ and contravariant basis
$\mathbf{e}^{p}$ can be defined as (Riley et al., 2004; Casotto & Fantino,
2009)

$\nabla=\mathbf{e}^{p}\frac{\partial}{\partial u^{p}}\;.$ (31)

Applying the gradient operator to the potential generates a new tensor of one
order higher, which is the first-order tensor (vector) describing the magnetic
field

$\nabla V=V_{,p}\mathbf{e}^{p}\;,$ (32)

where we use the comma notation $V_{,p}=\partial V/\partial u^{p}$ to denote
partial differentiation. The components that arise by applying an operator
such as eq.(31) which is represented in the covariant basis, are here referred
to as the “natural” components. In contrast, describing the gradient in terms
of the normalized basis vectors, the ”physical” components of the first order
gradient can be expressed as

$V_{,p}^{*}=\frac{1}{h_{p}}\frac{\partial V}{\partial
u^{p}}=\frac{1}{h_{p}}V_{,p}\;.$ (33)

That is, we use the ”*” notation to denote the physical components represented
in the normalized basis, and distinguish it from the natural components. The
metric scale factor $h_{p}$ is determined by the elements of the metric tensor
$g_{pq}=\mathbf{e}_{p}\cdot\mathbf{e}_{q}$ (which completely characterize any
curvilinear coordinate system) as $h_{p}=\sqrt{g_{pp}}$. Note that the metric
tensor also facilitates the conversion between covariant and contravariant
bases (Riley et al., 2004; Casotto & Fantino, 2009). Applying the gradient
operator to eq.(31), produces the second-order gradient operator

$\displaystyle\nabla\nabla$
$\displaystyle=\mathbf{e}^{q}\frac{\partial}{\partial
u^{q}}\left(\mathbf{e}^{p}\frac{\partial}{\partial u^{p}}\right)$
$\displaystyle=\mathbf{e}^{p}\mathbf{e}^{q}\left(\frac{\partial^{2}}{\partial
u^{p}\partial u^{q}}-\Gamma_{pq}^{s}\frac{\partial}{\partial u^{s}}\right),$
(34)

where $\Gamma_{pq}^{s}$ denotes the Christoffel’s symbols of the second kind
(an array of numbers describing the derivatives of the covariant basis vector
along that same basis), which can be expressed in terms of the metric tensor
as (Riley et al., 2004, p. 814)

$\Gamma_{pq}^{s}=\frac{1}{2}g^{st}\left(\frac{\partial g_{qt}}{\partial
u^{p}}+\frac{\partial g_{pt}}{\partial u^{q}}-\frac{\partial g_{pq}}{\partial
u^{t}}\right)\;.$ (35)

Applying the operator eq.(34) to the magnetic potential $V$ generates the
second-order magnetic gradient tensor elements which can be written using the
Christoffel’s symbols

$V_{;pq}=\left(\frac{\partial^{2}V}{\partial u^{p}\partial
u^{q}}-\Gamma_{pq}^{s}\frac{\partial V}{\partial u^{s}}\right)\;.$ (36)

Here the semicolon notation is used to denote covariant differentiation. In
addition, indices which occurs after a semicolon, means differentiation. In
physical components eq.(36) is written as

$V_{;pq}^{*}=\frac{1}{h_{p}h_{q}}V_{;pq}\;.$ (37)

Recall here, that the ”*” denotes the physical components. An essential aspect
of the first- and second-order tensors is how their physical elements
$V_{,p^{\prime}}^{*}$ or $V^{*}_{;p^{\prime}q^{\prime}}$ in one coordinate
system $u^{p^{\prime}}$ transforms to a new coordinate system $u^{p}$ (Casotto
& Fantino, 2009; Riley et al., 2004, p. 811)

$\displaystyle V_{,p}^{*}$
$\displaystyle=\frac{h_{p^{\prime}}}{h_{p}}\frac{\partial
u^{p^{\prime}}}{\partial u^{p}}V_{,p^{\prime}}^{*}$ (38) $\displaystyle
V^{*}_{;pq}$
$\displaystyle=\frac{h_{p^{\prime}}}{h_{p}}\frac{h_{q^{\prime}}}{h_{q}}\frac{\partial
u^{p^{\prime}}}{\partial u^{p}}\frac{\partial u^{p^{\prime}}}{\partial
u^{q}}V^{*}_{;p^{\prime}q^{\prime}}\;,$ (39)

where the partial derivatives $\partial u^{p^{\prime}}/\partial u^{p}$ are
expressed by the Jacobian matrix. The Jacobian matrix times the metric scale
factor term, can be regarded as a rotation matrix such that we may re-write
eqs.(38) and (39)

$\displaystyle V^{*}_{,p}$ $\displaystyle=V^{*}_{,p^{\prime}}R$ (40)
$\displaystyle V^{*}_{;pq}$
$\displaystyle=RV^{*}_{;p^{\prime}q^{\prime}}R^{T}\;,$ (41)

having the transformation matrix determined as

$R=\frac{\partial u^{p^{\prime}}}{\partial u^{p}}D\;,$ (42)

where
$D=diag(h_{p^{\prime}}/h_{p})=diag(h_{1^{\prime}}/h_{1},h_{2^{\prime}}/h_{2},h_{3^{\prime}}/h_{3})$
is a diagonal $3\times 3$ matrix of the scale factor ratios between the two
coordinate systems using the physical components. Thus equations (38) and (39)
(equivalently eqs.(40) and (41)) allow us to transform the tensors in one
coordinate system, for instance the global $(\tilde{x},\tilde{y},\tilde{z})$,
to another, for instance $(r,\theta,\phi)$. Note here, that for the
equivalence of eqs.(40) and (41) in the natural components, the transformation
matrix is defined by eq.(42) but without the matrix $D$, i.e. without the
metric correction. Let us now consider the two transformations:

  * a)

Transformation from the global Cartesian $(\tilde{x},\tilde{y},\tilde{z})$ to
the spherical system $(r,\theta,\phi)$

  * b)

Transformation from the spherical system $(r,\theta,\phi)$ to the local system
$(z,x,y)$

First, we specify the inner products of the basis vectors, the covariant
metric tensors for the Cartesian system , which is valid for both the global
and local Cartesian coordinate systems having basis vectors
$(\hat{\mathbf{i}}_{1},\hat{\mathbf{i}}_{2},\hat{\mathbf{i}}_{3})$ and
$(\hat{\mathbf{e}}_{1},\hat{\mathbf{e}}_{2},\hat{\mathbf{e}}_{3})$,
respectively

$g_{pq}=\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&1\end{array}\right]\;.$ (43)

and the spherical system, having the un-normalized basis vectors
$\mathbf{e}_{p}$, i.e.
$(\mathbf{e}_{r},\mathbf{e}_{\theta},\mathbf{e}_{\phi})$

$g_{pq}=\left[\begin{array}[]{ccc}1&0&0\\\ 0&r^{2}&0\\\
0&0&r^{2}\mathrm{sin}^{2}\theta\end{array}\right]\;.$ (44)

Note, that the un-normalized vectors
$(\mathbf{e}_{r},\mathbf{e}_{\theta},\mathbf{e}_{\phi})$ are the same as the
vectors of the basis $\mathbf{e}_{1}$, $\mathbf{e}_{2}$ and $\mathbf{e}_{3}$
mentioned in the introduction of this appendix. Thus the metric scale factors
of the Cartesian system become

$h_{\tilde{x}}=1,\quad h_{\tilde{y}}=1,\quad h_{\tilde{z}}=1\;,$ (45)

and for the spherical system

$h_{r}=1,\quad h_{\theta}=r,\quad h_{\phi}=r\mathrm{sin}\theta\;.$ (46)

The Christoffel’s symbols for the spherical coordinate system determined by
eq.(35) yields 27 values of which 9 are non-zero

$\displaystyle\Gamma_{pq}^{1}$
$\displaystyle=\left[\begin{array}[]{ccc}0&0&0\\\ 0&-r&0\\\
0&0&-r\mathrm{sin}^{2}\theta\end{array}\right]$ (50)
$\displaystyle\Gamma_{pq}^{2}$
$\displaystyle=\left[\begin{array}[]{ccc}0&\frac{1}{r}&0\\\ \frac{1}{r}&0&0\\\
0&0&-\mathrm{cos}\theta\mathrm{sin}\theta\end{array}\right]$ (54)
$\displaystyle\Gamma_{pq}^{3}$
$\displaystyle=\left[\begin{array}[]{ccc}0&0&\frac{1}{r}\\\
0&0&\frac{\mathrm{cos}\theta}{\mathrm{sin}\theta}\\\
\frac{1}{r}&\frac{\mathrm{cos}\theta}{\mathrm{sin}\theta}&0\end{array}\right]\;.$
(58)

Note that for both the global and local Cartesian systems, the Christoffel’s
symbols are zero as the metric tensor is the identity matrix. In case a) the
Jacobian matrix between the spherical coordinates $u^{p}=(r,\theta,\phi)$ and
the Cartesian coordinates
$u^{p^{\prime}}=\tilde{x}^{p}=(\tilde{x},\tilde{y},\tilde{z})$ is derived
considering the position vector in eq.(30)

$\left(\frac{\partial u^{p^{\prime}}}{\partial
u^{p}}\right)=\frac{\partial(\tilde{x},\tilde{y},\tilde{z})}{\partial(r,\theta,\phi)}=\begin{pmatrix}\frac{\partial\tilde{x}}{\partial
r}&\frac{\partial\tilde{x}}{\partial\theta}&\frac{\partial\tilde{x}}{\partial\phi}\\\
\frac{\partial\tilde{y}}{\partial
r}&\frac{\partial\tilde{y}}{\partial\theta}&\frac{\partial\tilde{y}}{\partial\phi}\\\
\frac{\partial\tilde{z}}{\partial
r}&\frac{\partial\tilde{z}}{\partial\theta}&\frac{\partial\tilde{z}}{\partial\phi}\end{pmatrix}=\begin{pmatrix}\mathrm{sin}\theta\mathrm{cos}\phi&r\mathrm{cos}\theta\mathrm{cos}\phi&-r\mathrm{sin}\theta\mathrm{sin}\phi\\\
\mathrm{sin}\theta\mathrm{sin}\phi&r\mathrm{cos}\theta\mathrm{sin}\phi&r\mathrm{sin}\theta\mathrm{cos}\phi\\\
\mathrm{cos}\theta&-r\mathrm{sin}\theta&0\end{pmatrix}\;.$ (59)

In case b) the Jacobian matrix between the local Cartesian coordinates
$u^{p}=(z,x,y)$ and the spherical coordinates $u^{p^{\prime}}=(r,\theta,\phi)$
is found by considering the position vector in the local Cartesian basis

$\mathbf{r}=z\hat{\mathbf{e}}_{1}+x\hat{\mathbf{e}}_{2}+y\hat{\mathbf{e}}_{3}=z\frac{\partial\mathbf{r}}{\partial
r}+\frac{x}{r}\frac{\partial\mathbf{r}}{\partial\theta}+\frac{y}{r\mathrm{sin}\,\theta}\frac{\partial\mathbf{r}}{\partial\phi}\;.$
(60)

From this, we can compute the Jacobian matrix

$\left(\frac{\partial u^{p^{\prime}}}{\partial
u^{p}}\right)=\frac{\partial(r,\theta,\phi)}{\partial(z,x,y)}=\begin{pmatrix}\frac{\partial
r}{\partial z}&\frac{\partial r}{\partial x}&\frac{\partial r}{\partial y}\\\
\frac{\partial\theta}{\partial z}&\frac{\partial\theta}{\partial
x}&\frac{\partial\theta}{\partial y}\\\ \frac{\partial\phi}{\partial
z}&\frac{\partial\phi}{\partial x}&\frac{\partial\phi}{\partial
y}\end{pmatrix}=\begin{pmatrix}1&0&0\\\ 0&\frac{1}{r}&0\\\
0&0&\frac{1}{r\mathrm{sin}\,\theta}\end{pmatrix}\;.$ (61)

The first-order tensor (i.e. the magnetic field vector) in the spherical polar
coordinates is given by eq.(33) using the metric scale factors eq.(46)

$\nabla
V=V_{,r}\mathbf{\hat{e}_{r}}+\frac{1}{r}V_{,\theta}\mathbf{\hat{e}_{\theta}}+\frac{1}{r\mathrm{sin}\theta}V_{,\phi}\mathbf{\hat{e}_{\phi}}\;.$
(62)

The physical components of first-order magnetic tensor elements in the local
Cartesian system, can be derived from eqs.(40) and (42), using the Jacobian
matrix eq.(61). Here the transformation matrix $R$ becomes the identity matrix
because the Jacobian matrix, ($\partial u^{p^{\prime}}/\partial
u^{p}=\partial(r,\theta,\phi)/\partial(z,x,y)$), is the inverse of the $D$
matrix. Therefore the physical elements in the local Cartesian system are
identical to those in the spherical polar coordinates given by the relations

$V^{*}_{;z}=V^{*}_{,r}\quad,\quad V^{*}_{;x}=V^{*}_{,\theta}\quad,\quad
V^{*}_{;y}=V^{*}_{,\phi}\;.$ (63)

In addition, from eqs.(33) and (46), and using the fact that eq.(45) also
holds for the local Cartesian system having ($h_{z}=h_{x}=h_{y}=1$), we obtain
the relations

$V_{,z}=V_{,r}\;,\quad V_{,x}=\frac{1}{r}V_{,\theta}\;,\quad
V_{,y}=\frac{1}{r\mathrm{sin}\,\theta}V_{,\phi}\;.$ (64)

The second-order tensor in the spherical polar coordinates is given by eq.(36)
using the Christoffel’s symbols from eq.(58) and metric scale factors eq.(46)

$\displaystyle\nabla\nabla V$
$\displaystyle=V_{,rr}\mathbf{\hat{e}_{r}}\mathbf{\hat{e}_{r}}+\left(\frac{1}{r}V_{,\theta
r}-\frac{1}{r^{2}}V_{,\theta}\right)\mathbf{\hat{e}_{r}}\mathbf{\hat{e}_{\theta}}+\left(\frac{1}{r\mathrm{sin}\theta}V_{,\phi
r}-\frac{1}{r^{2}\mathrm{sin}\theta}V_{,\phi}\right)\mathbf{\hat{e}_{r}}\mathbf{\hat{e}_{\phi}}$
$\displaystyle+\left(\frac{1}{r}V_{,r\theta}-\frac{1}{r^{2}}V_{,\theta}\right)\mathbf{\hat{e}_{\theta}}\mathbf{\hat{e}_{r}}+\left(\frac{1}{r^{2}}V_{,\theta\theta}+\frac{1}{r}V_{,r}\right)\mathbf{\hat{e}_{\theta}}\mathbf{\hat{e}_{\theta}}+\left(\frac{1}{r^{2}\mathrm{sin}\theta}V_{,\phi\theta}-\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi}\right)\mathbf{\hat{e}_{\theta}}\mathbf{\hat{e}_{\phi}}$
$\displaystyle+\left(\frac{1}{r\mathrm{sin}\theta}V_{,r\phi}-\frac{1}{r^{2}\mathrm{sin}\theta}V_{,\phi}\right)\mathbf{\hat{e}_{\phi}}\mathbf{\hat{e}_{r}}+\left(\frac{1}{r^{2}\mathrm{sin}\theta}V_{,\theta\phi}-\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi}\right)\mathbf{\hat{e}_{\phi}}\mathbf{\hat{e}_{\theta}}+\cdots$
$\displaystyle+\left(\frac{1}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi\phi}+\frac{1}{r}V_{,r}+\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}\theta}V_{,\theta}\right)\mathbf{\hat{e}_{\phi}}\mathbf{\hat{e}_{\phi}}\;.$
(65)

At the position $P$ (being the GVO target point), we can relate the physical
components of eq.(65) to the local Cartesian system by using eq.(41) having
the transformation matrix equal to the identity matrix, i.e. $R=I$. This leads
to the following relations

$\displaystyle V^{*}_{;zz}$ $\displaystyle=V^{*}_{;rr}\quad
V^{*}_{;zx}=V^{*}_{;r\theta}\quad V^{*}_{;zy}=V^{*}_{;r\phi}$ (66)
$\displaystyle V^{*}_{;xz}$ $\displaystyle=V^{*}_{;\theta r}\quad
V^{*}_{;xx}=V^{*}_{;\theta\theta}\quad V^{*}_{;xy}=V^{*}_{;\theta\phi}$ (67)
$\displaystyle V^{*}_{;yz}$ $\displaystyle=V^{*}_{;\phi r}\quad
V^{*}_{;yx}=V^{*}_{;\phi\theta}\quad V^{*}_{;yy}=V^{*}_{;\phi\phi}\;.$ (68)

Since the natural and physical component are identical in the local Cartesian
system (in fact in any Cartesian system), we have the relations

$\displaystyle V_{,zz}$ $\displaystyle=V^{*}_{;rr}\quad
V_{,zx}=V^{*}_{;r\theta}\quad V_{,zy}=V^{*}_{;r\phi}$ (69) $\displaystyle
V_{,xz}$ $\displaystyle=V^{*}_{;\theta r}\quad
V_{,xx}=V^{*}_{;\theta\theta}\quad V{,xy}=V^{*}_{;\theta\phi}$ (70)
$\displaystyle V_{,yz}$ $\displaystyle=V^{*}_{;\phi r}\quad
V_{,yx}=V^{*}_{;\phi\theta}\quad V_{,yy}=V^{*}_{;\phi\phi}\;.$ (71)

By using eqs.(69) to (71) together with eq.(65), we can infer the relations
between the gradient tensor elements in the local Cartesian system and the
gradient tensor described in the spherical system. We can do this by
recognizing, that the components in front of the elementary unit tensor (e.g.
$\hat{\mathbf{e}}_{r}\hat{\mathbf{e}}_{r}$) of eq.(65), are the physical
components of the $\nabla\nabla V$ tensor, so that:

$\displaystyle V_{,zz}$ $\displaystyle=V_{,rr}$ $\displaystyle V_{,xz}$
$\displaystyle=\frac{1}{r}V_{,\theta r}-\frac{1}{r^{2}}V_{,\theta}$
$\displaystyle V_{,yz}$ $\displaystyle=\frac{1}{r\mathrm{sin}\theta}V_{,\phi
r}-\frac{1}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi}$ $\displaystyle V_{,zx}$
$\displaystyle=\frac{1}{r}V_{,r\theta}-\frac{1}{r^{2}}V_{,\theta}$
$\displaystyle V_{,xx}$
$\displaystyle=\frac{1}{r^{2}}V_{,\theta\theta}+\frac{1}{r}V_{,r}$
$\displaystyle V_{,yx}$
$\displaystyle=\frac{1}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi\theta}-\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi}$
$\displaystyle V_{,zy}$
$\displaystyle=\frac{1}{r\mathrm{sin}\theta}V_{,r\phi}-\frac{1}{r^{2}\mathrm{sin}\theta}V_{,\phi}$
$\displaystyle V_{,xy}$
$\displaystyle=\frac{1}{r^{2}\mathrm{sin}\theta}V_{,\theta\phi}-\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi}$
$\displaystyle V_{,yy}$
$\displaystyle=\frac{1}{r^{2}\mathrm{sin}^{2}\theta}V_{,\phi\phi}+\frac{1}{r}V_{,r}+\frac{\mathrm{cos}\theta}{r^{2}\mathrm{sin}\theta}V_{,\theta}\;.$
(72)

In order to express the gradient tensor in Cartesian coordinates, we note that
the metric tensor becomes the identity matrix meaning that the metric scale
factors $h_{p}$ becomes unity, and all of the Christoffel’s symbols becomes
zero such that the gradient tensor is given by eq.(36)

$\displaystyle\nabla\mathbf{B}$
$\displaystyle=-\left(\begin{array}[]{lll}\frac{\partial^{2}V}{\partial
z^{2}}&\frac{\partial^{2}V}{\partial x\partial
z}&\frac{\partial^{2}V}{\partial y\partial z}\\\ \frac{\partial^{2}V}{\partial
z\partial x}&\frac{\partial^{2}V}{\partial
x^{2}}&\frac{\partial^{2}V}{\partial y\partial x}\\\
\frac{\partial^{2}V}{\partial z\partial y}&\frac{\partial^{2}V}{\partial
x\partial y}&\frac{\partial^{2}V}{\partial
y^{2}}\end{array}\right)=-\left(\begin{array}[]{lll}V_{,zz}&V_{,zx}&V_{,zy}\\\
V_{,xz}&V_{,xx}&V_{,xy}\\\
V_{,yz}&V_{,yx}&V_{,yy}\end{array}\right)=\left(\begin{array}[]{lll}\left[\nabla
B\right]_{zz}&\left[\nabla B\right]_{zx}&\left[\nabla B\right]_{zy}\\\
\left[\nabla B\right]_{xz}&\left[\nabla B\right]_{xx}&\left[\nabla
B\right]_{xy}\\\ \left[\nabla B\right]_{yz}&\left[\nabla
B\right]_{yx}&\left[\nabla B\right]_{yy}\end{array}\right)\;,$ (82)

where the minus sign comes from defining the field as the negative gradient of
the potential. We recall that the semicolon notation denotes tensor elements
following Casotto & Fantino (2009), and not the second order spatial
derivatives. However, in the case of the gradient tensor in Cartesian
coordinates, these two are equivalent cf. eq.(82).