Co-estimation of core and lithospheric magnetic fields by a maximum entropy method

Satellite observations of the geomagnetic field contain signals generated in Earth's interior by electrical currents in the core and by magnetized rocks in the lithosphere. At short wavelengths the lithospheric signal dominates, obscuring the signal from the core. Here we present details of a method to co-estimate separate models for the core and lithospheric fields, which are allowed to overlap in spherical harmonic degree, that makes use of prior information to aid the separation. Using a maximum entropy method we estimate probabilistic models for the time-dependent core field and the static lithospheric field that satisfy constraints provided by satellite observations while being consistent with prior knowledge of the spatial covariance and expected magnitude of each field at its source surface.


INTRODUCTION
Earth's magnetic field is a result of sources located both within the Earth and above its surface in the upper atmosphere (Chapman & Bartels 1940;Langel & Hinze 1998;Olsen & Stolle 2012).Spherical harmonic analysis indicates that internal sources are responsible for the majority of the field (Gauss 1839;Hulot et al. 2015).The spatial power spectrum of the internal field at Earth's surface is steep at low degree (up to approximately spherical harmonic degree 13) and essentially flat at higher degree (Lowes 1974), indicating two sources: one deep within the planet and one located near to the surface (Voorhies et al. 2002;Voorhies 2004).These deep and shallow sources are thought to correspond to the core dynamo and lithospheric magnetisation.
Measurements of Earth's magnetic field from space have provided an increasingly detailed picture of the magnetic field due to internal sources.The MAGSAT mission (Langel et al. 1982) delivered the first set of vector measurements with global coverage allowing the change in the slope of the spatial power spectra, between the wavelengths where core and respectively lithospheric sources dominate, to be definitively observed (Langel & Estes 1982).With more recent satellite missions, in particular the CHAMP (Reigber et al. 2002) and Swarm (Friis- Christensen et al. 2006;Olsen & Floberghagen 2018) missions, it is possible to determine the internal field spectrum out to beyond degree 130 (Maus 2010;Olsen et al. 2017).Knowledge of the small scale core field (which we define here as above spherical harmonic degree 13) has on the other hand advanced rather little since the time of MAGSAT due to it being obscured by the lithospheric field.The conventional approach is to estimate a single internal spherical harmonic field model and then to truncate at degree 13 to study the core field, for example Core and lithospheric fields 3 when plotting maps of the radial field at the core-mantle boundary (CMB) (e.g.Cain et al. 1989;Olsen et al. 2014;Sabaka et al. 2020).
Truncation at a fixed spherical harmonic degree has however limitations when seeking to isolate the core field.Abrupt truncation in spectral space may cause ringing in physical space (Whaler & Gubbins 1981;Gubbins 2007).Furthermore, the lithospheric field does, of course, not stop at degree 14, and there will be some contribution to the internal field below 13 from lithospheric sources.Most seriously all information on the small scale core field above degree 13 is lost.
In the 1980s it was suggested that a better way to estimate the core field would be to minimize suitable norms of the field complexity at the CMB (Shure et al. 1985;Gubbins & Bloxham 1985).This approach, known as spatial regularization of the field, has been widely adopted for studying the core field over historical (Bloxham et al. 1989;Jackson et al. 2000) and paleomagnetic timescales (Korte et al. 2011;Panovska et al. 2018) when data coverage is sparse; it makes use of prior information from seismology on the depth of the CMB, along with asking for a field that is simple in a specific way (as measured by a chosen regularization norm) at the source radius.A drawback is that traditional regularizations norms, such as the squared value of the radial field or the horizontal gradient of the radial field integrated over CMB, or Ohmic heating norms (Gubbins 1976;Jackson et al. 2000), strongly penalize small length scales and typically cause the spatial power spectrum to decay in an unphysical fashion above degree 13 (Backus 1988;Buffett & Christensen 2007).
Geodynamo simulations for which the magnetic Reynolds number is of order 1000, as expected in Earth's core (Christensen & Tilgner 2004;Lhuillier et al. 2011), involve localized, high amplitude, flux features and spatial spectra at the CMB that are rather flat, decreasing only very slowly at spherical harmonic degrees 10 to 30 (see, for example, Schaeffer et al. 2017;Aubert et al. 2017;Sheyko et al. 2018).Jackson (2003) and Jackson et al. (2007) showed that regularization norms based on the entropy of the radial field at the CMB allowed the estimation of core fields with flatter spatial spectra and localized, high amplitude flux features.The entropy regularization technique was adapted to timedependent spherical harmonic field models by Gillet et al. (2007) and applied to satellite observations from the Ørsted, SAC-C and CHAMP missions by Finlay et al. (2012).A drawback in these studies was the need to abritrarily pick a value for the so-called default parameter (the magnitude of the radial field expected in the absence of data constraints) that controlled the width of the entropy distribution and hence the sharpness of the field structures (Maisinger et al. 2004;Jackson et al. 2007).Jackson (2003) and Jackson et al. (2007) focused on default values around 10 µT for the core field, while Gillet et al. (2007) used 30 µT.Despite many desirable features, entropy-based field reconstruction techniques have been little exploited in subsequent years in part due to doubts as to how to pick the troublesome default parameter.
An important conceptual step forward in co-estimating core and lithospheric field sources was made by Holschneider et al. (2016).They suggested how various field sources (including the core and lithospheric fields) could be co-estimated within a Bayesian framework making use of prior information, for example on the expected source depth and its correlation (or covariance) structure.This approach has been used to develop temporal sequences of field models using a Kalman filter algorithm, being applied to the modelling of satellite and ground magnetic observations by Ropp et al. (2020), Ropp & Lesur (2023) and Baerenzung et al. (2020Baerenzung et al. ( , 2022)).Using a simple correlation function and treating the source depth as a free parameter Baerenzung et al. (2020) were able to construct stable maps of the field at the CMB up to spherical harmonic degree 20, although it was found to be difficult to reliably separate the large scale lithospheric field.More detailed prior information on the covariance between spherical harmonic coefficients in dynamo simulations has also been used in combination with observation-based internal field models up to degree 13 to infer the core field up to degree 30 (Aubert 2015(Aubert , 2020)).On the theoretical side Baratchart & Gerhards (2017) have shown that core and lithospheric fields can be formally separated if the lithospheric field sources are localized to a sub-region of the spherical surface.Non-Gaussian field distributions thus seem to aid the separation of fields from different sources, as is well known in other contexts such as independent component analysis (e.g.Hyvärinen & Oja 2000).
Here we build on the above studies and seek to estimate separate models for the core and lithospheric fields within a Bayesian framework using a maximum entropy method that accounts for spatial covariances found in first principles simulations of the core dynamo and the lithospheric magnetisation.Similar maximum entropy based techniques have previously been applied to signal separation problems in cosmology (Hobson et al. 2010).Section 2 sets out details of our Bayesian model estimation scheme and specifies the prior information used.Section 3 describes the satellite and ground magnetic observations employed.Section 4 presents our results, with Appendix A collecting findings from a synthetic test based on a similar data and modelling setup.We conclude in Section 5 with a discussion of what has been achieved and suggestions for future improvements of the method.

Geomagnetic field model
We model Earth's magnetic field B as a potential field, representing it by the gradient of a scalar potential V such that Core and lithospheric fields 5 with V int the potential due to internal sources and V ext that due to external sources.Both core and lithospheric sources contribute to V int , we represent each by a separate spherical harmonic expansion where (r, θ, ϕ) are geocentric spherical polar coordinates, a is the Earth's mean spherical reference radius, n is the degree of the spherical harmonic, m the order of the spherical harmonic and P m n (cos θ) are associated Legendre functions.g L n,m and h L n,m are spherical harmonic coefficients describing the lithospheric field, assumed here to be static, and here considered up to a maximum degree of t) are time-dependent spherical harmonic coefficients for the core field, considered up to a maximum degree N C = 30.These are expanded in time using a B-spline basis, of order 6 and with a 0.5 year knot spacing, where B k is the kth basis function of the order 6 B-splines.We collect the coefficients describing the core field in a vector and the coefficients describing the lithospheric field in a vector m L = g L n,m , h L n,m .Note that the core and lithospheric field representations overlap between spherical harmonic degrees 1 and 30, additional prior information is therefore needed in order to separate them.
As in the CHAOS-7 model (Finlay et al. 2020) these internal field coefficients are supplemented by model coefficients m ext describing the external field, and coefficients m q describing the in-flight alignment of the vector magnetometers on each satellite, to give the full model vector m = m C , m L , m ext , m q T .

Entropic priors for the core and lithospheric fields
We make use of prior information regarding the radial component of the core and lithospheric fields at their respective source surfaces, the CMB and Earth's surface.The radial field is evaluated at each The latent variables x C and x L therefore describe the core and lithospheric radial fields at their source surfaces in a space where their elements are normalized and decorrelated, as is appropriate for the application of maximum entropy methods (Maisinger et al. 2004).
It is assumed that x C and x L are each described by an entropic probability density function where S is the information entropy for variables that can take both positive and negative values (Gull & Skilling 1990;Hobson & Lasenby 1998) where M is in our case the number of grid points on the spherical surface, ψ i = x 2 i + 4ω 2 , and ω is a so-called 'default' parameter that defines the width of the entropy function.The information entropy function S is a measure of the amount of uncertainty inherent in the distribution of values x (Shannon 1948;Jaynes 2003), the form we use follows from requirements of subset independence, coordinate invariance, system independence and scaling (Skilling 1988).In the geomagnetic context it can be thought of as measuring the number of ways a given distribution of radial field on the source surface can be arranged from elementary flux bundles (Jackson 2003;Jackson et al. 2007); fields with larger entropy are simpler in the sense that they can be arranged in more ways.
Assignment of an entropic prior is argued to be an appropriate choice in the absence of precise information as to the form of a prior pdf (Skilling 1989;Hobson et al. 1998), and it is compatible with possibly non-Gaussian distributions of x.Maximizing the entropy essentially broadens the distribution of x as much as possible without violating the available observational constraints.The resulting distribution therefore agrees with what is known but expresses maximum uncertainty with respect to all other matters (Jaynes 1968).Maximizing the entropy does not introduce additional correlations amongst the latent variables (Gull & Skilling 1984).
The factors λ C and λ L appearing in the entropic pdfs are scaling factors.We are able to set these equal to 1 because x C and x L have already been normalized via the transform of b C and b L Core and lithospheric fields 7 to the latent space (Hobson et al. 2010).The transform to latent space using the a-priori covariance functions also ensures the entropy is computed from uncorrelated variables, an important condition for the maximum entropy method.
To put into practice the above scheme we require prior information concerning the covariance structure of the radial fields on the source surfaces and the expected widths of the distributions of x for each source.We obtain these from first principles simulations of the core dynamo and the lithospheric magnetisation.Full details are provided in Otzen (2022) only a short summary is given here.
For the core field prior, we use an ensemble of realizations of the core field produced by versions of the coupled-Earth dynamo of Aubert et al. (2013).This numerical dynamo is known for producing field structures and patterns of secular variation similar to those observed over the past centuries.To start with we used a collection of radial fields realizations, well separated in time, generated by the original version of the coupled-Earth dynamo (Aubert et al. 2013) that has been used in previous field modelling and data assimilation studies (Barrois et al. 2017;Ropp et al. 2020).To this we added radial field realizations from a long run of an updated version (71% of path) of the coupled-Earth dynamo (Aubert & Gillet 2021).Although these two cases involve different control parameters they lie on a path through control parameter space along which the field morphology is essentially invariant (Aubert et al. 2017).We finally augmented our set of realizations by carrying out rotations of the simulated core fields by an arbitrary amount in longitude, this was possible because the covariance functions we use do not depend on longitude and this enabled us to work with a larger ensemble.In all this resulted in an ensemble of 5688 core field realizations up to spherical harmonic degree 30.To be more consistent with the observed field we also adjusted the dipole fields from the dynamo simulations replacing the n = 1 coefficients with random samples from normal distributions with mean values of ḡC 1,0 = −29000 nT, ḡC 1,1 = 0 nT, and hC 1,1 = 0 nT and standard deviations of 5000nT, 3000nT, and 3000nT respectively, the latter being similar to those seen in the dynamo realizations.The spread of the power spectra from the resulting ensemble encompasses the observed internal field (e.g.Finlay et al. 2020) up to spherical harmonic degree 13 (Otzen 2022).
Our prior for lithospheric field comes from simulations of the lithospheric magnetization based on the forward modelling scheme developed by Hemant & Maus (2005), with revised oceanic magnetisation according to Masterton et al. (2012) and subduction zone magnetisations following Williams & Gubbins (2019).We produced an ensemble of 6000 realizations of the lithospheric field by (i) varying the crustal thickness within a range given by published crustal thickness models (Nataf & Ricard 1996;Reguzzoni & Sampietro 2015), (ii) varying parameters of the remanent vertically integrated magnetisation model (Masterton et al. 2012;Williams & Gubbins 2019), and (iii) using stochastic perturbations generated using a Gaussian random field approach.The observed power spectra for the Core and lithospheric fields 9 lithospheric field, for example from the LCS-1 model (Olsen et al. 2017), lies within the spread of the power spectra of this ensemble (Otzen 2022).The simulated fields were generated up to spherical harmonic degree 255, but here we only used them up to degree 120.
From each ensemble member, we evaluated the radial field on an approximately equal area grid on the source surface.We used HEALPix (Górski et al. 2005) grids with 3072 points for the core field, and 49152 points for the lithospheric field, which are suitable for representing fields up to the spherical harmonic truncation level chosen for each source.Based on these gridded values we computed empirical semi-variograms as a function of angular distance on the spherical surface and fit covariance functions to these.We used a multi-quadratic covariance function for the lithospheric field (Gneiting 2013) and a combination of a multi-quadratic covariance function and a spline function for large distances for the core field.The resulting covariance functions, along with the empirical covariances of the ensemble members, are shown in the top row of Fig. 1 along with their corresponding power spectra at the source surfaces generated using these covariance functions are shown in the middle row.
These covariance functions provide us with the a-priori expected covariance structure for our core and lithospheric fields.Fig. 2 shows example realizations of the core and lithospheric field generated using these covariance functions.
We also make use of the distribution of the radial field at the sources surfaces from our core and lithospheric field ensembles, after transformation to the latent space (see equation ( 6) and the related discussion).These distributions and relevant statistics are presented in the bottom row of Fig. 1 . In particular we use expected absolute values, < |x| >= |x| p(x) dx, calculated using the mean empirical pdfs shown in Fig. 1, to define the latent space default parameters i.e. we set ω =< |x| >, separately for the core and lithospheric fields.
The time-dependence of the core field is represented by a 6th order B-spline representation smoothed by third time derivative regularization, a standard choice in time-dependent geomagnetic field models when one wishes to study field accelerations (see e.g.Olsen et al. 2014).This is formally equivalent within the Bayesian framework to assuming a-priori that the spherical harmonic coefficients are realizations of a continuous process of the form (Wahba 1990) where f n,m (t) is a zero mean Gaussian process (see e.g.Rasmussen & Williams 2006) specified by the covariance function where t a and t b are arbitrary times, u is a dummy integration variable, and we have used the notation (z) + = z for z ≥ 0 and (z) + = 0 otherwise.K C,0 n,m , K C,1 n,m and K C,2 n,m are constants associated with constant, linear and quadratic time-dependences that can be different for each spherical harmonic coefficient.σ 2 n is the variance of the process that depends on the spherical harmonic degree n and is related to the choice of regularization parameter in the standard field modelling framework.

Likelihood of geomagnetic observations
Turning to the observations, we assume the geomagnetic measurements are contaminated by unmodelled signals that follow a long-tailed Huber error distribution.The appropriate likelihood function is then where χ 2 (m) = e T W e is a robust (Huber weighted) misfit norm.e = d − g(m) are the residuals between the ground and satellite magnetic observations d and the associated model predictions g(m).
where C e is the a-priori data error covariance matrix and W h a diagonal matrix that implements robust (Huber) weights and has elements is the a-priori expected error for the ith datum and c = 1.5 is a constant (Constable 1988;Olsen 2002;Sabaka et al. 2004).
Core and lithospheric fields 11

Estimation of the model posterior probability density function
Applying Bayes theorem the posterior probability function is This can be maximized by minimizing the loss function Here S tav (x C ) = 1 approximates the information entropy of the (spatially decorrelated) CMB radial field averaged over time (Gillet et al. 2007) by summing the entropy at each discrete L is assumed to be static.C T is the a-prior model (temporal) covariance matrix, which includes the temporal prior information from Eqn. 10 for the core part of the model.
The minimization is achieved iteratively using a Newton-type descent algorithm.The (k + 1)th estimate of the posterior mean model is obtained based on the model at the previous kth iteration by where is the Jacobian matrix of partial derivatives of the forward model for each datum with respect to the model parameters, evaluated using the model parameters m k .We have used a notation similar to that of Stockmann et al. (2009) to define matrices, related to the Hessian matrix of the second-order partial derivatives of the entropy function, of the form and vectors, related to the first order partial derivatives of the entropy function, of the form 14) the superscripts C and L refer to the core and lithospheric fields respectively, while subscript k denotes that the computation is carried out using model parameters from the previous k th iteration.In α C and β C the expressions for α and β at each epoch t p must be averaged over time in the same way as S tav is defined above.
Minimizing measures of the data misfit and the temporal complexity while maximizing the entropy of the (decorrelated) core and lithospheric radial fields at their source surfaces results in internal fields that fit the observations and are compatible with the temporal prior but allow high dynamic ranges of x at the source surfaces while satisfying the spatial covariance properties of the core and lithospheric priors.
After convergence of the above scheme we describe the dispersion of the posterior distribution using an approximate model covariance matrix defined about the maximum of the posterior pdf, computed from the Hessian of the loss function Φ(m) by (e.g.Tarantola 2005;Hobson et al. 1998) where A, W, α C and α L are the values of A k W k , α C k and α L k from the final iteration, when the scheme is considered to have converged.

OBSERVATIONS
The models reported here are built from a dataset of geomagnetic observations similar to that used to construct the CHAOS-7 geomagnetic field model (Finlay et al. 2020), but restricted to the period between 2005 and 2020.
Observations from the CHAMP, Cryosat-2 and Swarm A and B satellites are used, three-component vector measurement for quasi-dipole latitudes up to 55 degrees and scalar intensity data at higher latitudes.Level 3 CHAMP magnetic field data, Cryosat-2 L1b magnetic field data (FGM 1, the dataset used in the CHAOS-7 model, here pre-calibrated using CHAOS-7) and Swarm L1b magnetic field data (version 0601) are used with 1 minute sampling for CHAMP and Cryosat-2 and 2 minute sampling for each Swarm satellite.We also used along-track gradients from CHAMP, Swarm A and Swarm B; gradients are particularly useful for constraining the high degree lithospheric field.Geomagnetic quiet-time selection criteria were employed such that Kp ≤ 2 0 , d|RC|/dt ≤ 2 nT/hr (Olsen et al. 2014), averaging over the previous 2 hours the merging electric field at the magnetopause E m ≤ 0.8 mV/m, the interplanetary magnetic field (IMF) B Z > 0 , and IMF B Y is less than 3 nT in the northern hemisphere and greater than -3 nT in the southern hemisphere.Only data from dark conditions (sun at least 10 degrees below the horizon) were used.A more detailed description of these data selection criteria is found in Finlay et al. (2020).In addition to satellite observations, as in CHAOS-7 we used annual differences of revised monthly means from ground observatories are used, as in the CHAOS-7 model, based on hourly mean values from the BGS AUX OBS database, version 0129.A stacked histogram of the number of vector field observations used versus time is presented in Fig. 3.

Implementation and model diagnostics
We now briefly document here some details regarding the practical implementation before moving on to the results.In order to make our modelling setup as close as possible to the CHAOS field modelling scheme, in C −1 T we consider sub-matrices associated with calculating quadratic norms of (i) the 3rd time derivative of the internal radial field integrated over the core surface and and throughout the model time span (formally equivalent to the the prior defined in Eqn. 9) (ii) the acceleration of the core field at the model endpoints and (iii) temporal first differences of the estimated offsets of the external dipole in solar-magnetic coordinates (related to imperfections in the RC index).The related hyperparameters denoted by λ i3 , λ i2e , and λ sm are fixed throughout this study and implicitly included within C −1 T .They were chosen so as to produce time-variations similar to the CHAOS-7 field model although for simplicity we did not use degree-dependent tapering or treat zonal terms differently.
The relevant hyperparameters related to the entropy default parameters for the dimensionless latent variables, ω C and ω L , were set to values 0.412 and 0.422 based on the expected absolute values of x C and x L from the distributions of the dimensionless latent variables found in the prior ensembles (see Fig. 1,bottom).The adopted hyperparameters are collected in table 1.
The full estimated model, including the time-dependent core field to degree 30, static lithospheric field to degree 120, external field parameters and alignment parameters for each satellite, consists of 49495 parameters in all.We started the iterative model estimation scheme with model parameters for the core taken from CHAOS-6.9 up to degree n = 12, and for the lithosphere from the LCS-1 model at degree n = 16 and above.The small-scale core field and large-scale lithospheric field were otherwise initialized with zeros.After 24 iterations the largest change in a model parameter relative to its amplitude was 0.0037% and no further change was seen in CMB maps of the posterior mean core field.Below we refer to the resulting model, including co-estimated core and lithospheric parts, as model CL.
For comparison, we also built a more traditional CHAOS-type field model, with a single timedependent internal field up to spherical harmonic degree 20 and a static internal field for degrees 21 to 120, using the same external field parameterization and hyperparameters, covering the same period, and based on the same dataset as used to build the model CL.For this model we considered the estimated core field to be the time-dependent internal field up to degree 13, as has been the standard practise when interpreting the CHAOS model (e.g.Olsen et al. 2014).
Table 2 collects Huber-weighted means and RMS residuals between the vector field data and the model predictions (in nT), comparing model CL and our CHAOS-type reference model.Model CL fits the satellite and ground data overall to a similar level as the CHAOS-type model, while simultaneously minimizing the information entropy of the (spatially-decorrelated) time-dependent core and lithospheric fields.The two models are found to have very similar temporal regularization norms, which is not surprising at they were built using the same temporal hyperparameters.The non-dimensional information entropy norms for the decorrelated core and lithospheric fields, S C tav and S L , after 24 iterations were respectively 0.94 and 5.31 for the CL model.

The core-mantle boundary field
The spatial power spectra (Lowes 1966;Mauersberger 1956) of model CL at the CMB in 2020 is presented in Fig. 4. The obtained posterior mean model closely follows the internal field from CHAOS-7 up to degree 11 but contains slightly less power at degree 12 and 13.Here and below comparisons to CHAOS-7 used version CHAOS-7.9.Internal field models (such as the CHAOS model) are usually truncated at degree 13 when carrying out interpretations at the CMB since above this degree their CMB spectra diverge as they also include signals from the lithospheric field.The posterior realizations from the model CL core field have power spectra that are approximately flat out to degree 30 and do not diverge.The posterior mean model shows a gradual drop in power for degrees 15 to 19 and a slight increase again for degrees 20 to 22. Above degree 22 the power in the posterior mean model drops to much lower levels indicating that field realizations essentially average to zero at higher degrees where almost all the observed signal comes from the lithospheric field.
The power in the lithospheric field realizations and mean models mapped down to the CMB is also shown in Fig. 4, these do diverge.The core and lithospheric field spectra cross between degree 14 and 16 for realizations of model CL, and at degree 15 for the mean models.Note that the estimated lithospheric field is presented only for degree 2 and above, at degree 1 it is not well separated from the core field, which we believe is a consequence of the entropy function of the decorrelated latent variables not being greatly affected by changes in the dipole field.
Maps of the radial magnetic field at the CMB in 2020.0 from model CL, for the posterior mean core field model and four example posterior realizations all truncated at degree 22, are presented in Siberia in 2020 is found to be more localized and stronger than in models truncated at degree 13.This feature has moved north-westwards between 2005 and 2020, as seen in Fig. 6 which shows the posterior mean map up to spherical harmonic degree 22 at a sequence of times.
Concerning reversed flux patches in the South Atlantic, we find evidence for two reversed patches under South Africa adjacent to strong norm flux patches under central Africa (see also their time evolution in Fig. 6).Regarding the reverse flux region under the Southern Atlantic ocean, there are several distinct reversed flux concentrations visible within this region, which are observed to evolve separately.

The lithospheric field at Earth's surface
Although not our main focus here, we present for reference details of our co-estimated lithospheric field, which includes an estimate of the large scale of the lithospheric field which is usually neglected.The cross-over between the mean core and lithospheric field models also occurs at degree 15 at Earth's surface.In Fig. 9 we present a map of the posterior mean lithospheric field from model CL at Earth's surface for degree 2 to 120.

Modelled secular variation and comparisons with ground observatories
To document the time dependence in model CL and show that this is also reasonable, Fig. 10 presents example time-series of the first time derivative (or secular variation, SV) of spherical harmonic coefficients from model CL's core field model.It shows the posterior mean, 5000 posterior realizations, and the CHAOS-7 model for reference.The time-dependent SV of g 0 1 is slightly smoother in model CL than in CHAOS-7, agrees well at intermediate degrees, following the same trends and showing some (expected) differences close to the endpoints.Model CL generally shows lower amplitude changes in SV at high degree.This behaviour is expected because the CHAOS-7 model tapered its regularization to lower values at high degrees.The dispersion of the posterior realizations also increases with spherical harmonic degree.We see no evidence for unrealistic features in the SV coefficients of model CL.
This conclusion is supported by comparisons of SV from model CL with observed SV data at ground observatories, for example as shown in Fig. 11.Model CL predictions agree well with the annual differences of monthly means for stations from Africa (M'Bour), Europe (Niemegk) and in the Pacific (Honolulu), with the fit being similar to that of CHAOS-7 although with somewhat smoother time variations.The posterior realizations give an indication of the formal uncertainties in the SV predictions of model CL.

Insights from a synthetic test
A key question is the extent to which our method is able to retrieve the core field above degree 13.
To investigate this issue we carried out a synthetic test where the true core and lithospheric fields were known.We took as input a time-dependent core field up to spherical harmonic degree 30 from a dynamo simulation, along with a synthetic lithospheric field based on a simulation of the induced and remanent lithospheric magnetisation up to degree 120.Magnetic data were synthesized at the same locations, and for the same field components, as in the observed dataset and this was inverted using the maximum entropy co-estimation scheme described in 2. Full details and the results from this synthetic test are collected in Appendix A. The synthetic data used in the test is consistent with our prior, however the utilized prior information is rather weak (involving only the source radii and isotropic spatial covariance functions); the main purpose of the test is to investigate what level of separation of the core and lithospheric fields is possible using our approach.Fig. A1 presents the resulting spatial spectra at the CMB and at Earth's surface while Fig. A2 compares side-by-side the radial field at the CMB from the estimated posterior mean core field model and the input dynamo field for increasing truncation degrees of 13, 16, 19 and 22.The posterior mean model obviously has less power than the dynamo synthetic truth from degree 16 to 22 with the missing power at small scales particularly obvious at low latitudes where the dynamo solution is most complex.
Despite under-estimation of the power at degrees 16 to 22, Fig. A2 shows the estimated posterior mean model does contain useful information on CMB field structures above spherical harmonic degree 13.Field structures remain coherent as power is added from degree 13 to 22, with some important details recovered.For example, in the southern polar region, under the Australian-Antarctic basin (near latitude 60 • S, longitude 135 • E) there is a localized intense flux feature present in both the dynamo model and in the posterior mean; this is weaker and smeared when the CMB field is truncated at degree 13.Similarly in the northern polar region under Siberia and Alaska, both intense normal polarity features and reversed flux features are better retrieved in the posterior mean model to degree 22 compared to a more conventional model truncated at degree 13.
At low latitudes a number of intense features are better retrieved in the estimated posterior mean model to degree 22 than in the model truncated at degree 13, for example the flux concentration under central Africa and the positive-negative pair of flux patches arranged north-south across the equator under India and under the equator south of Mexico.However some smaller scale features, for example a positive flux concentration in the dynamo model under the equator at 15 degrees west, are still poorly retrieved.Isolated small scale flux features with no imprint on larger scales, and which change rapidly in time, are not well recovered.
Flux features recovered in the posterior mean core field model up to degree 22 are however always related to features in the true dynamo field.We do not find any evidence for artefacts due to leakage of the lithospheric field.The reconstructed fields are by construction as simple as possible (in terms of maximizing the information entropy) while satisfying the observational constraints.The adopted temporal prior (spline smoothing) also involves strong time-averaging over small length scales that is absent in the dynamo, this will contribute to the loss of detail on small scales.We note that interpretations in terms of the model resolution matrix (e.g.Bloxham et al. 1989) are not straightforward here.With the available satellite data coverage we are able to well resolve the internal field up to high spherical harmonic degree.It is however constraints from prior information that allow us to partially separate the core and lithospheric fields.
Variations are found in the synthetic test results regarding the amount of small scale power in the posterior mean model as the data constraints change.For example there is less power in the small scale field in the gap between the CHAMP and Swarm missions when only high altitude and lower quality satellite data from Cryosat-2 was available.This is an expected consequence of the maximum entropy method -at times when the data constraints are weaker the posterior mean model becomes simpler and the model uncertainty at high degree larger.

Interpretations of the inferred CMB field
Returning to the model CL derived from the real observations, in Fig. 12 we present maps of the radial magnetic field at the CMB from the posterior mean model, for increasing truncation degrees of 13, 16, 19 and 22.As in the synthetic test, power adds coherently with increasing degree with some structures becoming more localized and intense.This is also the case for the SV at the CMB, which we have examined in manner similar to Holme et al. (2011)  Reversed flux features are found in the northern polar region under the Canadian Arctic, centred under the Queen Elizabeth Islands, and also under the Nansen Basin (between Spitzbergen and the new Siberian Islands).In the Southern polar region there is an extended but weak reverse flux feature under Eastern Antarctica southwards from Africa, this feature, also present in models truncated at degree 13, persists to higher degree and is interesting as it crosses the tangent cylinder.Reverse flux features inside the tangent cylinder could be related to a strong omega effect driven by strong azimuthal flows inside the tangent cylinder and associated flux expulsion, similar features have been seen in turbulent dynamos (Schaeffer et al. 2017;Sheyko et al. 2018).
At low latitudes the strong flux feature found under the equator between Africa and south America in models truncated at degree 13 is split into two features, as is a normal feature under central Africa.
Such splitting of low latitude flux features has been suggested in previous attempts to retrieve the core field above degree 13 (Baerenzung et al. 2020;Otzen et al. 2022), and it is consistent with the patterns of core surface SV retrieved above degree 13 (Finlay et al. 2020).Many of the field concentrations at low latitudes seem to occur in oppositely signed pairs and they are observed to move westwards.A possible explanation could be that they are the signature of toroidal flux being expelled from the core at low latitudes and subsequently propagating as a wave (e.g.Aubert et al. 2013Aubert et al. , 2022)).Such features are ubiquitous at low latitudes in strongly forced geodynamo simulations when the viscosity is sufficiently low (Sheyko 2014;Schaeffer et al. 2017;Aubert & Finlay 2019).

Features of the inferred large-scale lithospheric field
In the synthetic test reported in Appendix A we were also able to test the retrieval of the lithospheric field.The lithospheric field at the Earth's surface generally compares well with the input synthetic truth lithospheric field, albeit with less power below degree 15 and above degree 90 (Fig. A1 and Fig. A3).
Of particular interest is whether or not any details of the synthetic truth large scale lithospheric field were retrieved.Fig. A4 shows that, somewhat surprisingly, some details of the large scale lithospheric field can be recovered, albeit with reduced amplitude.The largest anomalies in the recovered large scale lithospheric field, for example between the North Pole and the Bering strait, near Australia and in north-Eastern Europe are also present in the synthetic truth model.On the other hand some prominent structures are missing or incomplete, for example in the Atlantic-Indian-Antarctic basin or under North America, and the recovered amplitudes are too low.Fig. 13 shows a similar map of the large scale lithospheric field from the model CL posterior mean derived from the real data.It shows strong anomalies in the northern part of Eastern Europe, Australia, around eastern Antarctica and under eastern North America, which are known locations of strong continental magnetic anomalies.

Limitations and future prospects
The aim of this study was to better separate core and lithospheric magnetic fields; it has only been partially successful.The recovered small scale core field and large scale lithospheric field lack power.
Use of more informative priors, if these can be justified, would certainly bring improvements.For example, we used a single spatial covariance function for each source, which assumes field structures that are statistically the same for all locations on the sphere.The ensembles of prior fields from the dynamo and magnetisation models could instead be used to build full covariance matrices characterising the statistical covariances between all locations on the sphere for each source.Such dense spatial covariance matrices have already been used by other authors (Gillet et al. 2019;Ropp et al. 2020;Istas et al. 2023), in the context of lower resolution core field and flow modelling.Once sufficiently large prior ensembles are available this will be a relatively simple extension of the method presented above.
A concern with this approach is that imperfect aspects of the dynamo and magnetisation simulations might be mapped into the estimated field models, it was for this reason we started here by using rather simple information from the prior ensembles.
Regarding the small scale core field, the spline-based temporal prior employed here prevents the recovery of rapid changes on small length scales.This limitation could be particularly serious at low latitudes.This situation could be remedied by adopting temporal priors that better reflect the expected Otzen et al.
A variant of the approach presented here is to take the spherical harmonic (Gauss) coefficients of the internal potential as input data for the separation into core and lithospheric fields rather than satellite data.This has some computational advantages and may prove useful in future applications, further details and an example are presented in Appendix B. Our scheme could also be applied to dedicated studies of the lithospheric field.This would require improving the lithospheric prior to allow for more power at small length scales, use of higher data sampling rates, and use of Swarm east-west gradients that were not included in this study.
Much is still to be learnt regarding the small scale core field.The maps of the posterior mean core field presented here show how information on the core dynamo is lost when CMB fields are truncated at degree 13.On the other hand our ability to retrieve the small scale core field, and avoid lithospheric field contamination, depends on correctly formulating and utilizing prior information regarding the sources.Further effort is needed on how best to extract reliable prior information from a variety of simulations of the core dynamo and the lithospheric magnetisation.Improved Bayesian field modelling requires prior ensembles that are both informative and broadly representative of the diversity of possible core and lithospheric fields.Here we report results from a synthetic test designed to test the extent to which we can retrieve the small scale core field above spherical harmonic degree 13 and the large scale lithospheric field below degree 13 using the scheme described in Section 2.
We use as input a time-dependent core field taken from a preliminary version of the dynamo model assimilation runs described by (Aubert 2023).This simulation contained realistic core field structures and time dependence up to degree 30 and was not contained in the prior ensemble.For the synthetic lithospheric field we used one realizatation from our set of prior magnetisation models, based on the forward models by Hemant & Maus (2005), Masterton et al. (2012) and Williams & Gubbins (2019) with the perturbations described by Otzen (2022), but generated separately and not included in the ensemble used to construct the prior statistics.From these models we synthesized data at the same times and positions, and with the same measured components, as the real satellite and ground data described in Section 3. Gaussian noise was added using a standard deviation of 3.0nT for CHAMP, 5.0nT for CryoSat-2, 2.0nT for Swarm and 3.0nT/yr for ground observatory data.
Inversions were then carried out in exactly the same way as for the real data.The same priors, and same fixed regularization settings the same starting models and the same number of iterations (24) were performed.
In Fig. A1 we present the Lowes-Mauersberger spectrum at the core surface (top) and at the Earth's surface (bottom), with the synthetic truth marked in the dashed line.In Fig. A2 we present how the estimated posterior mean and the synthetic truth for the core field changes with the truncation degree.Although some small scale feature are lost and the amplitude is reduced, the posterior mean model above degree 22 retrieves more details of real features.

6
source surface on an approximately equal area grid, and values are collected into vectors b C and b L for the core and lithospheric fields respectively.Such knowledge of the radial field at the source surface completely defines the potential due to an internal source.These are related to the spherical harmonic model coefficients discussed in the previous section by b C (t p ) = G C,tp m C and b L = G L m L (5) Otzen et al. with G C,tp and G L being matrices that synthesize the radial field from the relevant spherical harmonic model coefficients, for the core field at some epoch t p .Knowledge regarding the spatial covariance of each field at its source surface is provided in the form of a-priori model covariance matrices C C and C L , with lower triangular Cholesky factors L C and L L that can be used to transform b C and b L to latent variables x C and x L such that b

Figure 1 .
Figure 1.Statistical properties of the prior ensemble of core and lithospheric fields.Top row: Empirical covariance functions (a) for the CMB radial field, (b) for the lithospheric radial field at Earth's surface showing ensemble members (grey), and estimated covariance functions (black).For computational reasons empirical covariances are shown only out to 18 degrees for the lithospheric field.Middle row: power spectra for 5000 realizations generated from the estimated covariance functions (c) at the CMB, and (d) at Earth's surface with the CHAOS-7 model (up to degree 13) and from the LCS-1 model (above degree 16) for reference.Bottom row: empirical probability density functions for radial magnetic field values on approximately equal area grids at (e) the CMB and (f) Earth's surface, derived from the prior ensembles of core and lithospheric magnetic fields respectively, after transformation to an uncorrelated and normalized latent space.
Figure 2. Example realizations of radial magnetic fields generated using the estimated covariance functions for the core and lithospheric fields.The core fields are shown at the core surface (top row) and the lithospheric fields are shown at Earth's surface (bottom row).These illustrate the a-prior correlation structures assumed for the core and lithospheric fields.Note the covariance models used to generate these prior fields are rotationally invariant.

Figure 3 .
Figure 3. Stacked histogram showing number of vector field observations employed per month between 2005 and 2020.Colours indicate the data source.

Figure 4 .
Figure 4. Lowes-Mauersberger spectra showing power as a function of spherical harmonic degree at the coremantle boundary in 2020.0 for model CL.Solid black lines with circles show the estimated posterior mean core and lithospheric field models which cross at degree 15.Grey lines show realizations based on the adopted prior covariance model, cyan lines show posterior realizations of the core field and purple lines posterior realizationsof the lithospheric field.Degree 1 is not shown for the lithospheric field as it is not well separated from the core field.CHAOS-7 (up to degree 13) and LCS-1 (at degree 16 and above) are shown for reference.

Fig. 5 .
Fig. 5.A similar map from the CHAOS-type reference model, truncated in the conventional fashion at degree 13, is shown for reference.The posterior realizations contain more power at small length scales, but all realizations agree on the larger-scale structure as represented by the mean model.More details of the CMB field structures are evident in the model CL posterior mean compared to the traditional CHAOS-type model truncated at degree 13.Some low latitude flux concentrations are split, see for example the two strong positive radial field features near the equator between Africa and South America which are usually interpreted as single feature (as in the CHAOS-type model).A strong high latitude flux feature under Siberia, located under the Taymyr peninsula in central northern

Figure 5 .
Figure 5. Radial magnetic field maps at the core-mantle boundary in 2020 for model CL.Estimated posterior mean model (top left), a CHAOS-type model constructed from the same dataset (top right), and four example posterior realizations.All models are truncated at degree 22 except for the CHAOS-type model which is truncated at degree 13.The common features across all posterior models, as captured in the mean model, provide information on the core field to beyond degree 13.

Figure 6 .Figure 7 .
Figure 6.Time sequence of maps of radial field at the core-mantle boundary in 2006, 2010, 2014 and 2019 from the model CL posterior mean truncated at degree 22. Changes are due both to the true evolution of the core field and changes in observational constraints across the epochs.

Fig. 8
Fig.8shows the spatial power spectra of the co-estimated core and lithospheric fields from model CL at the Earth's mean spherical reference surface.The spread in the posterior realizations is very small (almost invisible in the plot) at degrees 17 up to 70, indicating that the lithospheric field is very well constrained by the observations for these degrees.At lower degrees the posterior spread increases, with the mean model containing lower power than most of the individual realizations.The posterior spread also increases above degree 70, becoming as large as the spread in the prior above degree 110 by which point the mean model contains less power than any of the posterior realizations.

Figure 8 .
Figure 8. Lowes-Mauersberger spectra showing power as a function of spherical harmonic degree at Earth's mean spherical reference radius in 2020.0 for model CL.Solid black lines with circles show the estimated posterior mean core and lithospheric field models which cross at degree 15.Grey lines show prior realizations based on the adopted spatial covariance model, cyan lines show posterior realizations of the core field and purple lines posterior realizations of the lithospheric field.CHAOS-7 (up to degree 13) and LCS-1 (at degree 16 and above) are shown for reference.

Figure 9 .
Figure 9. Map of estimated posterior mean lithospheric field at Earth's spherical reference radius, for degrees 2 to 120.

Figure 10 .
Figure 10.First time derivative (secular variation) of selected spherical harmonic coefficients of the core field over the model timespan from 2005 to 2020.Top left: dg 0 1 /dt, top right: dh 2 3 /dt, bottom left: dg 1 7 /dt and bottom right: dh 14 14 /dt.The estimated posterior mean model is shown as the black solid line, posterior realizations in grey and for reference CHAOS-7 in red.

Figure 11 .
Figure 11.Comparison of secular variation predicted by model CL for the posterior mean model (black solid line) and posterior realizations (blue lines) with annual differences of revised monthly means at selected ground observatories (thin black line with dots).CHAOS-7 is shown for reference in red.Left column shows the time derivative of the radial field, middle column the time derivative of the southward field component and right column the time derivative of the eastward field component.Top row: M'bour observatory (MBO) from low latitudes in west Africa, middle row: Niemegk observatory (NGK) from mid-latitudes in Europe and bottom row: Honolulu observatory (HON) from mid-latitudes in the Pacific.

Figure 12 .
Figure 12.Maps of radial magnetic field at the core-mantle boundary in 2020.0 from the model CL posterior mean, as a function of the truncation degree, up to degree 13 (top left), to degree 16 (top right), to degree 19 (bottom left) and to degree 22 (bottom right).Note how the power adds coherently to existing features already evident at lower truncation degrees.
(not shown); SV features seen at lower degree remain coherent although more small scale SV features begin to appear by degree 22 indicating this is at the edge of what we are able to reliably study.At high latitudes strong normal flux features are found close to the where the inner core tangent cylinder intersects the CMB (at latitudes 69.6 degrees north and south).There is a particularly intense flux concentration under the Taymyr Peninsula in Siberia in 2020, and a number of high amplitude features arranged near latitude 60 degrees north under Greenland, Canada and eastern Siberia.In the Southern polar region there are strong normal flux features under Wilkes land in eastern Antarctica and under western Antartica near the Antarctic peninsula.Normal flux patches localized close to the tangent cylinder are consistent with the poloidal dipole field being produced by an alpha-effect in energetic eddies originating in vigorous convection close to the inner core boundary.Another possibility is that such eddies are spawned by powerful azimuthal flows inside the tangent cylinder occasionally ejected across the tangent cylinder (e.g.Schaeffer et al. 2017).The different locations of the flux concentrations in the northern and southern hemisphere seem difficult to explain in terms of purely columnar flows.

Figure 13 .
Figure 13.Estimated posterior mean large-scale lithospheric field (degrees 2 to 14) plotted on a spherical surface at Earth's mean reference radius.

Fig
Fig.A3presents a comparison of the estimated posterior mean lithospheric field and Fig.A4compares the recovered (posterior mean) large scale lithospheric field, and the input synthetic truth for the large scale lithospheric field at Earth's spherical reference radius, for degrees 2-14.

Figure A1 .
Figure A1.Lowes-Mauersberger spectra showing power as a function of spherical harmonic degree of magnetic fields at the CMB (top) and at Earth's surface spherical reference radius, (bottom) in 2020.0.The input synthetic truth models for the core and lithospheric fields are shown in the black dashed lines.Solid black lines with circles show the estimated posterior mean core and lithospheric field models which cross at degree 15.Grey lines show prior realizations based on the adopted spatial covariance model, cyan lines show posterior realizations of the core field and purple lines posterior realizations of the lithospheric field.CHAOS-7 (up to degree 13) and LCS-1 (at degree 16 and above) are shown for reference.

Figure A2 .
Figure A2.Maps of radial magnetic field at the core-mantle boundary from the estimated posterior mean model derived from the synthetic test dataset (left column) and the input truth dynamo model (right column), as a function of the truncation degree, up to degree 13 (top row), to degree 16 (second row), to degree 19 (third row) and to degree 22 (bottom row).

Figure A3 .
Figure A3.Map of the radial component of the estimated posterior mean lithospheric field on a spherical surface at Earth's reference radius, for degrees 2 to 120 (top) and the same quantity from the synthetic truth input model based on lithospheric magnetisation models (bottom).

Figure A4 .
Figure A4.Estimated posterior mean large-scale lithospheric field (degrees 2 to 14) derived from the synthetic dataset plotted at Earth's spherical reference surface (top) and the same quantity from the synthetic truth model based on a prior realization of lithospheric magnetisation (bottom).Note the change in the scale used for the two panels.

Table 2 .
Misfit statistics for vector field data in the non-polar region and scalar data in the polar region.