CM6: a comprehensive geomagnetic field model derived from both CHAMP and Swarm satellite observations

From the launch of the Ørsted satellite in 1999, through the CHAMP mission from 2000 to 2010, and now with the Swarm constellation mission starting in 2013, satellite magnetometry has provided excellent monitoring of the near-Earth magnetic field regime. The advanced Comprehensive Inversion scheme has been applied to data before Swarm and to the Swarm data itself, but now for the first time to all the satellite data in this new era, culminating in the CM6 model. The highlights of this model include not only a continuous core magnetic field description over the entire time period 1999 to 2019.5 in good agreement with the CHAOS model series, but the addition of two new oceanic tidal magnetic sources: the larger lunar elliptic semi-diurnal constituent N2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_2$$\end{document} and the lunar diurnal constituent O1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O_1$$\end{document}. CM6 is also the parent model of the NASA/GSFC candidates for the DGRF2015 and IGRF2020 in response to the IGRF-13 call. This paper provides a full report on the development of CM6.


Introduction
The "modern era" of satellite magnetometry can be said to have been initiated with the launch of the Danish satellite Ørsted on 23 February, 1999 (still in orbit, but not collecting magnetic measurements since summer 2014), followed by the "CHAllenging Minisatellite Payload" (CHAMP) launched on 15 July, 2000 and flying for 10 years before its demise in September, 2010. Other satellites have flown in the interim, such as SAC-C (data available for 2000 to 2004), leading to the European Space Agency (ESA) Swarm mission launched on 22 November, 2013. Swarm is a constellation of three satellites, two side-by-side low-fliers, referred to as "Alpha" and "Charlie", whose polar orbits allow for a unique East-West cross-track sampling of field gradients and a high-flier, referred to as "Bravo", whose polar orbit precesses so as to sweep out different local times to improve separation between internal and external current systems. When supplemented with measurements from ground-based observatories, in particular, observatory hourly mean values (OHMs), this fleet of magnetic monitors provides an excellent survey of Earth's magnetic field.
For more than two decades the magnetic modeling group at the NASA/Goddard Space Flight Center (GSFC), in collaboration with the group at the Technical University of Denmark (DTU), developed a modeling approach known as the "Comprehensive Inversion" (CI) (see Sabaka et al. 2013) in which the majority of dominant magnetic signals detected in the near-Earth regime are parameterized and subsequently co-estimated to obtain optimal signal separation, while taking into account both random and systematic errors. This has resulted in the production of a number of "Comprehensive Models" (CMs) using data from the pre-Swarm era (e.g., Sabaka et al. 2002Sabaka et al. , 2004Sabaka et al. , 2015 and related models based only on Swarm data (see Sabaka et al. 2018), which have been delivering a consistent set of Swarm Level-2 (L2) magnetic data products. While the named CMs have never included Swarm data and while the only satellite data included in the Swarm L2 models are from Swarm, the model reported on here, denoted "CM6", is the first CM to analyze data from the Ørsted, SAC-C, CHAMP and Swarm satellites, thus covering the entire "modern era".
The magnetic sources originally included in CI parameterizations are current systems in the core, lithosphere, ionospheric E and F-regions, magnetosphere, and associated induced currents. However, with the advent of high-quality satellite vector magnetic measurements, particularly from CHAMP, additional current systems were added such as the oceanic M 2 tidal field (see Sabaka et al. 2015Sabaka et al. , 2016Sabaka et al. , 2018. With the recent work of Grayver and Olsen (2019), it is clear that combinations of CHAMP and Swarm measurements are beginning to resolve some weaker tidal constituents, such as the larger lunar elliptic semi-diurnal constituent N 2 and the lunar diurnal constituent O 1 , and so these have now been included in the CI parameterization.
The development of CM6 has coincided with the call for candidates of the "International Geomagnetic Reference Field: the 13th Generation" (IGRF-13). Therefore, GSFC has submitted candidates for the Definitive Geomagnetic Reference Field: 2015 (DGRF2015) and the International Geomagnetic Reference Field: 2020 (IGRF2020) that are directly extracted from, in the case of the former, and linearly extrapolated from, in the case of the latter, the CM6 model. This paper reports on the development of the CM6 model by first presenting a description of the data selection procedure in "Data selection" section followed by a brief overview of the CI algorithm in "Methodology" section , including model parameterization and the estimation procedure, and ends with a discussion of the results in "Results and discussion" section, including details on the newly derived N 2 and O 1 oceanic tidal components and the IGRF-13 candidates.

Data selection
Both scalar-and vector-based measurements were used as provided by the Ørsted, SAC-C, CHAMP, and Swarm satellite missions and vector measurements from groundbased observatory hourly means.

Temporal selection
The temporal selection criteria for satellite scalar direct measurements and for their along-track (mainly in the north-south direction due to the polar orbits) sums and differences from the Ørsted, SAC-C, CHAMP, and Swarm missions were a 30-s cadence with alongtrack sums and differences computed from samples 15 s apart . The time periods used were 15 March, 1999   through 23 February, 2011 for Ørsted, 23 January,  2001 through 4 December, 2004 for SAC-C, 30 June,  2000 through 4 September, 2010 for CHAMP, and  25 November, 2013 through 31 December, 2018 for  Swarm. Vector direct measurements and their along-track sums and differences from the CHAMP and Swarm missions were taken at the same periods using a 30-s cadence with along-track sums and differences computed from samples 15 s apart. From the two sideby-side flying Swarm Alpha and Charlie satellites, cross-track (mainly east-west due to the orbit configurations) scalar and vector sums and differences were used from the period 17 April, 2014 through 31 December, 2018, at a sampling rate of 30 s. Cross-track sums and differences were computed from data where Alpha and Charlie were nearest in geocentric latitude at a temporal separation of at most 15 s (typically 4-10 s apart).
Ground observatory hourly mean values from 199 stations were used from the period January 1999 through April 2019. These data were checked and cleaned for trends, spikes and other errors (see Macmillan and Olsen 2013).

Environmental selection
The environmental selection criteria were the following: (1) Kp ≤ 2 − (3° for sums and differences); (2) dDst/dt ≤ 3 nT/h (for direct and sums and differences); (3) direct satellite vector and scalar measurements only from dark regions, i.e., the sun at least 10 • below the horizon; (4) satellite direct vector measurements equator-ward of 55 • quasi-dipole (QD) latitude; (5) vector sums and differences from the CHAMP mission were used when both double and single Star Tracker heads were available for attitude determination; and (6) Swarm data diverging from the CHAOS-6 x7 model by more than 100 nT in the scalar field or 500 nT in any of the vector field components were rejected. Data rejection takes place at each iteration through the identification and elimination of gross outliers as determined by an inspection of residual plots (see "Estimation procedure" section).
A stacked bar graph of measurements types and counts per year is shown in Fig. 1 where seven types are broken out: (1) North-South along-track scalar sums and differences; (2) East-West across-track scalar sums and differences; (3) North-South along-track vector sums and differences; (4) East-West across-track vector sums and differences; (5) satellite direct scalar measurements; (6) satellite direct vector measurements; and (7) observatory vector hourly means measurements. Note the absence of satellite measurements in the 2011−2012 year gap between the CHAMP and Swarm missions.

Model parameterization
The parameterization of the CM6 model follows closely that of the CIY4 model (Sabaka et al. 2018) except that now the oceanic tidal magnetic signals of the N 2 and O 1 constituents are included as well as the magnetometer alignment parameters of the CHAMP satellite. A synopsis of the parameterization is listed in Table 1. A brief description of each source field is given below, but more details can be found in Sabaka et al. (2018).

Core and lithospheric fields
The core and lithospheric magnetic fields are expressed as gradients of solid harmonic (SH) expansions of the internal branch solution to Laplace's equation in spherical coordinates corresponding to an SH degree truncation level of N max = 120 , where the first N SV = 18 degrees allow for secular variation (SV) in the form of order-4 B-splines spanning 1999.0 to 2019.5 with knots every 6 months giving a total of 45 parameters per SH coefficient, and for degrees above N SV the coefficients are constant. The expression for the corresponding core/lithospheric potential at time t and position r , corresponding to Earth-Centered Earth-Fixed (ECEF) spherical coordinates of radius, colatitude, and longitude (r, θ , φ) , is given by where the R{·} operator takes the real part of the expression only, (·) * is the complex conjugate of the argument, and Y m n is the surface SH of degree n and order m given by where a is the Earth mean-radius (6371.2 km) and P m n is the Schmidt semi-normalized associated Legendre function of degree n and order m. The complex Gauss coefficients, γ m n (t) , of degree n and order m are given by (2) Y m n (θ, φ) = P m n (cos θ) exp imφ, Fig. 1 Stacked bar graph of measurement counts per year used in deriving CM6, where "NS" and "EW" denote North-South along-track and East-West across-track, respectively, "DF/SF" and "DB/SB" denote scalar and vector differences and sums, respectively, and "SATF", "SATB", and "OHMB" denote satellite scalar, vector, and observatory vector measurements, respectively where b q is the qth cubic B-spline of the expansion and the epoch of the expansion is 2015.0. The general complex Gauss coefficient γ m n (t) is related to the familiar real Gauss coefficients g m n (t) and h m n (t) by γ m n (t) = g m n (t)+ih m n (t) , that is, g m n (t) is the real part and h m n (t) is the imaginary part and i = √ −1.

Oceanic tidal fields
The magnetic field from the oceanic principal lunar semi-diurnal constituent M 2 has been included in CM6 as the gradient of an internal potential (see Sabaka et al. 2015Sabaka et al. , 2016Sabaka et al. , 2018 whose SH degree truncation level is N max = 36 and each coefficient is sinusoidal in time with a 12.42060122 h periodicity with time t rendered with respect to Greenwich phase. The potential at time t and position r in the ECEF system is then where τ m n is the complex coefficient and Y m nω (�t, θ, φ) = P m n (cos θ) exp i(mφ + ω M 2 �t), with ω M 2 = 2π/12.42060122 rads/h. Following the work of Grayver and Olsen (2019), both the larger lunar elliptic semi-diurnal constituent N 2 and the lunar diurnal constituent O 1 were included in the CM6 parameterization in a similar form to that of M 2 , but with ω M 2 replaced by ω N 2 = 2π/12.65834751 rads/h for N 2 and ω O 1 = 2π/25.81933871 rads/h for O 1 .

Ionospheric field
The treatment of the ionospheric and associated induced magnetic fields is identical to that used in Sabaka et al. (2015Sabaka et al. ( , 2018 and is built from QD symmetric basis functions (Richmond 1995;Emmert et al. 2010) that reflect the conductivity structures found in the ionospheric E-region with an induced field due to its interaction with a 3-dimensional (3D) conductivity model where a surface layer containing continents and oceans is underlain by a 1-dimensional (1D) mantle (Kuvshinov and Olsen 2006), and updated in Kuvshinov (2008). More details can be found in Sabaka et al. (2018).
The secondary induced field is expressed in the spectral domain through transfer functions Q(ω) at frequency ω , which in this case are for 24, 12, 8, and 6 h periodicities. These periods are further modulated by an annual and semi-annual periodicity and by scaling from the 3-month running average of the F 10.7 solar radiation index. Longperiod induced variations are approximated by an infinite conductor at depth.

Magnetospheric field
The magnetosphere and associated induced magnetic fields are discretized into 1-h time bins during selected quiet periods within which the fields are treated as static external and internal SH expansions in dipole coordinates, respectively, to degree N max = 1 and order M max = 1 . This results in 114, 185 hourly bins covering 64% of the hours of the model time span from 1 January 1999 to 30 April 2019.

Alignment parameters
Alignment rotations between the vector magnetometer frame (VFM) and the spacecraft common reference frame (CRF) are estimated for Swarm Alpha, Bravo, and Charlie and for CHAMP when two Star Trackers are in use or one of two. These are parameterized in terms of three Euler angles representing rotations around the x-axis of the CRF followed about the new y-axis and then the new z-axis. The angles are treated as static in 10 day intervals.

Estimation procedure
Following Sabaka et al. (2018), the CM6 parameters are determined using the Gauss-Newton (GN) non-linear least-squares estimator (Seber and Wild 2003) with linear equality constraints and smoothing quadratic constraints whose kth step is is the residual vector of the data d with respect to the non-linear model vector is the Jacobian of the model vector evaluated at x k , x k are the adjustments to the current parameter vector x k , and W k ≡ W(x k ) is the data weight matrix. There are N q quadratic constraints, where P j is the jth a priori quadratic constraint matrix that, along with the Lagrange multiplier j , specifies the deviation of the solution from the preferred a priori model vector, which in the case of smoothers, is the zero vector. The matrix G is the linear equality constraint matrix which ensures that the induced SH time series associated with magnetospheric variations and the core field SV temporal basis functions are orthogonal. More details, including the solution formulation, are given in Sabaka et al. (2018). The starting model, x 0 , was obtained by first updating the CIY4 model of Sabaka et al. (2018), based upon the first 4 years of Swarm data, to the CIY5 model by adding the fifth year of Swarm data. The CIY5 model was then used to initialize all of the model parameters with the exception of the core field before 2013.9 and the CHAMP vector magnetometer alignment angles. The CHAOS-6-x8 model of Finlay et al. (2016) was used to initialize the core field from 1999 to 2013.9, but the CIY5 core field epoch of 2015.0 was adopted since it is anchored well within the Swarm data envelope. The CHAMP alignment angles were initialized to zero. A total of four iterations were performed.

Error-covariance
The W k matrix is assembled from error-covariance and error-bias considerations. These numbers are based upon residual analysis from previous models starting from OIFM (Olsen et al. 2000) and successively updated by various generations of the CHAOS and CM models. For vector OHM data, all local times and latitudes were used in the local North-East-Center (NEC) frame and if poleward of ±55 • QD latitude were assigned isotropic sigmas (or uncertainties) of 15 nT, while those equatorward of ±10 • were assigned isotropic sigmas of 7 nT, and elsewhere isotropic sigmas of 4 nT. Nightside scalar data from Ørsted and SAC-C at all latitudes were used and these were given sigmas of 4 nT. Scalar sums and differences from all local times and latitudes were used with the same weighting.
CHAMP scalar data from all latitudes were used, and vector data from the nightside at non-polar latitudes were used when two Star Tracker heads ("2ST") were available. Sigmas of 3 nT are used for scalar data and the vector isotropic factor. The attitude error treatment of Holme (2000), Holme andBloxham (1995, 1996) was applied which resolves the vector measurements into the BP3 ("HB-theory" frame) orthogonal coordinate system where "B" is along the predicted magnetic field direction, "P" is in the n×B direction where n is the unit vector along the CRF z-axis, and "3" completes the system. For 2ST data, 10 arcsec pointing and rotation angle sigmas were used. Scalar North-South sums and differences were used at all local times at all latitudes and vector North-South sums and differences were used at low-mid latitudes at all local times. The sums are assigned a 3 nT sigma isotropic error while differences use 0.3 nT and both use the same attitude error just mentioned. When only one Star Tracker head is available (designated either "AST" or "BST") then only vector data from the nightside at low-mid latitudes and North-South vector sums and differences at low-mid latitudes at all local times are used. The same uncertainty assignments are used for single heads except that the attitude error now comprises 10 arcsec pointing and 40 arcsec rotation angle sigmas.
Data from Swarm Alpha, Bravo, and Charlie were selected in exactly the same way as the dual-head Star Tracker CHAMP data. However, assigned sigmas differ in that 2.2 nT is used for scalar and vector data components and attitude error is not considered. The North-South sums are assigned a 2.2 nT sigma per component while differences use 0.3 nT. The East-West vector sums and differences from Alpha and Charlie are from lowmid latitudes and all local times, while the scalar components are from all latitudes. They too are assigned a 2.2 nT sigma per component for sums while differences use 0.3 nT. The sums and differences are analyzed in the local NEC frame, while the direct vector measurements are analyzed in the BP3 frame even though no attitude error is considered.
All data were weighted by sin θ , where θ is their respective geographic colatitude, to mitigate high-density data in the polar regions.
Robust estimation was employed via Iteratively Reweighted Least-Squares (IRLS) using Huber weighting (Constable 1988) and so for the ith iteration, if the kth residual, e i,k , was within cσ i of the mean of its distribution, in this case defined by a mean of zero, a standard deviation of σ i , and c = 1.5 , then it was treated as Gaussian noise with the previous uncertainty assignment. However, if it was outside of cσ i of the mean, then it was treated as Laplacian noise. Thus, weights were assigned according to To mitigate the effects of error-bias in certain measurement types, the scheme known as "Selective Infinite Variance Weighting" (SIVW) (see Sabaka et al. 2013) was used. It turns out that this can be interpreted as a variant of the Schmidt Type Consider Filter (Bierman 1977) where the variance of considered parameters tends to infinity. The application is based upon measurement type, sun position, and QD latitude range. More details are given in Sabaka et al. (2018), but Table 2 indicates how SIVW is applied in the form of "nominal" and "nuisance" parameters with respect to various data types across the core, lithospheric, and tidal parameter subspaces. Generally, if a particular data subset has a high signal-to-noise ratio for these signals, then they contribute to the "nominal" set, otherwise, to the "nuisance" set.
Table 2 CM6 SIVW application, where the "x" indicates the Quasi-Dipole (QD) latitude and sun position of the data type and which parameters it directly influences, and where "OHM" denotes observatory measurements, "Single", "Diffs" and "Sums" denote single, differenced and summed measurements, respectively, "NEC" and "BP3" denote the local North-East-Center and HB-theory frame for vector data, respectively, and "F" denotes scalar  et al. Earth, Planets and Space (2020) 72:80

Constraints
The number of explicit quadratic smoothing constraints used is N q = 8 along with the linear equality constraint. The ionospheric and magnetospheric/ induced constraints are the same as those used in CIY4 and include the night-side E-region currents, denoted as " P� J eq,MLT:21−05 2 2 � ", the smoothness of the diurnally varying portion of the currents, denoted as " P� ∇ 2 s J eq,p>0,mid−lat 2 2 � ", the Euclidean length of the magnetospheric/induced coefficients in each bin, denoted as " P p mag/ind 2 2 ", and the linear constraint that forces orthogonality of the induced magnetospheric field and the core SV through time, denoted as P�|p ind⊥core | 2 2 � . More details may be found in Sabaka et al. (2018). The lithospheric constraint is on the mean-square field magnitude |B| over Earth's mean surface at 6371.2 km for SH degrees 110 and above and is denoted as " P� B n≥110 2 � ". The biggest difference with CIY4 is that now only the mean-squared third time derivative of the radial component of the magnetic field, B r , at the Core-Mantle Boundary (CMB) at 3480 km radius over the entire time domain of the model, denoted as " P�| ... Br| 2 � ", is constrained as opposed to a mixture of second and third time derivatives. Finally, the first-differences of the Euler angle time series between 10-day bins comprising the magnetometer alignment for CHAMP 2ST, AST, and BST are constrained and are denoted as " P�| E 2ST | 2 � ", " P�| E AST | 2 � ", and " P�| E BST | 2 � ", respectively. Table 3 provides the associated values of for each constraint. Note that is infinite for linear equality constraints since these constraints may be expressed as the limit as →∞ of a related quadratic form.

Residual statistics
The weighted residual statistics for CM6 are shown for Ørsted, SAC-C, CHAMP 2ST, CHAMP AST, Swarm   Table 5 CM6 weighted residual statistics for CHAMP 2ST and AST, where "QD" denotes Quasi-Dipole, "NS" denotes North-South along-track, "DF/SF" denotes scalar differences and sums, respectively, "DB/SB", "DP/SP", and "D3/S3" denote vector differences and sums in the HB-theory frame, respectively, "F" denotes scalar, and "B", "P", and "3" denote vector measurements in the HB-theory frame, respectively Mean and RMS are in units of "nT"  Table 6 CM6 weighted residual statistics for Swarm Alpha and Bravo, where "QD" denotes Quasi-Dipole, "NS" denotes North-South along-track, "DF/SF" denotes scalar differences and sums, respectively, "DN/SN", "DE/SE", and "DC/SC" denote vector differences and sums in the local North-East-Center frame, respectively, "F" denotes scalar, and "B", "P", and "3" denote vector measurements in the HB-theory frame, respectively Mean and RMS are in units of "nT" where N is the number of measurements and e i and w i are the ith residual and Huber weight for a particular component, respectively, at the final iterate.

Table 7 CM6 weighted residual statistics for Swarm Alpha and Charlie, where "QD" denotes Quasi-Dipole, "EW" denotes East-West across-track, "DF/SF" denotes scalar differences and sums, respectively, and "DN/SN", "DE/SE", and "DC/SC" denote vector differences and sums in the local North-East-Center frame, respectively
Mean and RMS are in units of "nT" The Ørsted and SAC-C statistics, as expected, have larger RMSs for high and low versus mid QD latitude ranges, larger RMSs for light versus dark, and much smaller RMSs for differences versus sums. The CHAMP statistics also show similar behavior with respect to QD latitude, sun position, and data type. The single-head Star Tracker data AST and BST show very similar statistics and so the latter has been omitted for brevity. However, while data types involving the "B" component show similar statistics between 2ST and AST/BST, the statistics for the "P" and "3" component data types show consistently large RMSs for the single-head versus the dual-head Start Tracker cases. This makes sense since the attitude error in the bore-sight rotation angle is larger in the single-head cases and this degrades the "P" and "3" components.
As with the CIY4 model, Swarm Alpha and Charlie show very similar statistics as expected since they are in close proximity as the low satellite pair, and so the latter has been omitted for brevity, while Bravo shows slightly higher residuals. Swarm also has similar behavior with respect to QD latitude, sun position, and data type as the other satellites. The "B" and "F" component statistics are also in close agreement, as they are with CHAMP, which is expected since they are nearly in the same direction. The sums and differences statistics also follow the same patterns as seen in CIY4, as expected. It is also interesting that the North-South and East-West statistics are very similar in most cases reflecting a consistent level of measurement calibration between the Alpha and Charlie satellites.
Finally, the OHM statistics also exhibit the same patterns with respect to QD latitude and sun position indicating that these effects are seen at satellite altitudes as well as on ground.

Core and lithospheric fields
The core field of CM6 may be examined both spatially and temporally. The CHAOS-6-x8 model is used for comparison since it provided the starting state for the pre-Swarm core portion of CM6 and so discrepancies are more likely due to the incompatibility of the CHAOS-6-x8 core state with the data and theory incorporated in CM6. Figure 2 shows maps of the core radial field and associated SV at the CMB for SH degrees n = 1−13 at epoch 2016.0 from CM6 and the difference between CM6 and CHAOS. These maps look very similar to those shown for CIY4 in Sabaka et al. (2018) for the same epoch, as expected, particularly the dominant dipolar of the radial field itself. However, the difference maps are quite interesting, particularly for the radial field, which shows a strong zonal harmonic structure in dipole coordinates, perhaps g 0 13 given the number of zero-crossings. This may be related to how induction signals are handled in the two models. For CM6, both ionospheric and magnetospheric induction potentials are parameterized in terms of basic SH functions in dipole coordinates, whereas CHAOS does not consider ionospheric induction. The differences in radial SV appear to be more random. It is noteworthy that the CM6 SV at the CMB, like CHAOS, shows intense structures at high latitude under the Bering Strait and north Eastern Siberia. The differences in both the radial field and SV are at about the 10% level of the CM6 values.
As for the temporal structure of the core field, time series of CM6 and CHAOS-6-x8 Gauss coefficients and their time derivatives for SH degrees n = 1−3 from 1999 to 2019, the satellite portion of the CM6 time domain, are shown in Figs. 3 and 4, respectively. At the scale of the plots, the agreement appears to be very good between the Gauss coefficients. The notable deviations are in the low-degree zonal terms, g 0 1 and g 0 3 , which appear as a constant bias offset. This has been previously noted and attributed to the fact that only nightside data are used in CHAOS and no ionospheric induction is considered, which leads to a zonal signal pattern in QD coordinates in the differences.
The differences are more pronounced when comparing the core SV coefficients, as shown in Fig. 4. While Table 8 CM6 weighted residual statistics for OHMs, where "QD" denotes Quasi-Dipole and "N", "E", and "C" denote vector measurements in the local North-East-Center frame, respectively Mean  the overall patterns are in good agreement, it is clear that the CHAOS SV is smoother than that of CM6, particularly in the zonal terms ġ 0 1 and ġ 0 3 . There are also differences near the end points in the ḣ1 1 and ġ 0 2 terms. Again, this is most likely related to how induction is being handled in each of the models and to how strong the regularization constraints are being applied within both models.
The CM6 lithospheric field can be compared to that of the LCS1 (Olsen et al. 2017) and CIY4 models using the standard three metrics, the Lowes-Mauersberger spectrum, R n (r) , of Lowes (1966) where a and r are the reference and evaluation radii, respectively, and g m n and h m n are the real Gauss coefficients of the SH expansion; the degree correlation between two models (10) R n (r) = (n + 1) a r where g m n,k and h m n,k are the Gauss coefficients of model "k"; and the matrix of normalized coefficient differences (in %), S(n, m) where the subscripted "e" and "r" indicate the Gauss coefficients of the evaluated and reference models, respectively.
The top-left plot of Fig. 5 shows the R n spectra for CM6, LCS1, and CIY4 and the differences between CM6 and the other two models for SH degrees n = 20−100 at radius a = 6,371.2 km . According to this metric, the CM6 lithosphere is in better agreement with CIY4 up to (11)  SH degrees of about 60, but is in better agreement with LCS1 above degrees of about 85. However, the power in the model differences is about 3 orders of magnitude less than the signal strength for all degrees shown indicating generally good agreement between models. The agreement with CIY4 at low degrees is probably due to the commonality of Swarm data in the two models. The agreement with LCS1 above degree 85 is most likely due to the fact that CIY4 constrains its lithospheric expansion at SH degrees n≥85 . This is also corroborated by the ρ n (top-right) and S(n, m) (lower) plots. It may also be the case that CM6 tends to LCS1 because both models analyze CHAMP data, whereas CIY4 does not.

Fig. 5
The R n spectra (upper left) for LCS1, CIY4, and CM6 and the differences between CM6 and LCS1 and between CM6 and CIY4, the degree correlation, ρ n , (upper right) between CM6 and LCS1 and between CM6 and CIY4, and the normalized coefficient differences, S(n, m), of LCS1 (lower left) and CIY4 (lower right) with respect to CM6. All plots are for SH degrees n = 20−100 at radius a = 6,371.2 km Comparing the two plots, which use the same scale, it appears that the CM6 lithospheric field is in good agreement with that of LCS1. The discrepancies are manifested mainly in the auroral zones and along the dip equator, which are delineated by the red curves representing the QD latitudes of ±55 • and 0 • . Although neither CM6 nor LCS1 use dayside data equatorward of QD latitude ±10 • to determine the nominal lithosphere, LCS1 uses only dayside scalar gradients outside this region whereas CM6 uses dayside scalar and vector gradients outside this region and this perhaps may contribute somehow to the discrepancy patterns along the dip equator. In general, the definitive LCS1 model shows that the combination of Swarm and CHAMP measurements provides perhaps the best estimate of the small-scale crustal field at this time and the CMs are now following this course.

Oceanic tidal fields
The oceanic M 2 tidal signal has been detected in satellite data by Grayver and Olsen (2019), Sabaka et al. (2015, 2018), Tyler et al. (2003, while N 2 was crudely Fig. 7 The R n spectra of the time-averaged oceanic tidal magnetic field of M 2 from CM5, CIY4, CM6, and G&O for SH degrees n = 1−28 (top), of N 2 from CM6 and G&O for SH degrees n = 1−12 (lower left), and of O 1 from CM6 and G&O for SH degrees n = 1−12 (lower right), all at radius a = 6,371.2 km detected by Sabaka et al. (2016), but more convincingly by Grayver and Olsen (2019) (denoted as "G&O") who also detected O 1 . It should be noted that previous attempts to extract N 2 and O 1 using the CI approach with CHAMP data failed. However, with the advent of Swarm, these three constituents have now been co-estimated in CM6. To compare the models, a generalization of the classic R n spectrum of Lowes (1966), stated in Eq. 10, was developed by Sabaka et al. (2015Sabaka et al. ( , 2016 and defined for a particular tidal constituent as the mean-square magnitude of the magnetic field at SH degree n over a sphere of radius r and over the tidal constituent period given by where a = 6371.2 km and τ m n is defined in Eq. 4 for a particular constituent. The R n spectra are shown in Fig. 7 at radius a = 6,371.2 km for M 2 SH degrees n = 1−28 (top) from CM5, CIY4, CM6, and G&O, for N 2 SH degrees n = 1−12 (bottom-left) from CM6 and G&O, and for O 1 SH degrees n = 1−12 (bottom-right) from CM6 and G&O. These SH degree ranges were what was estimated by Grayver and Olsen (2019). For M 2 , the models show peak power at either degrees 5 or 6, and of the CI type models, the Swarm only model, CIY4, shows the largest value followed by the CHAMP/ Swarm mixture, CM6, and finally the CHAMP only model, CM5. The G&O model shows a more diminished power than all of the others. The power trends at higher degrees are all in basic agreement, although CM5 power does begin to diverge. For N 2 , CM6 peaks at degree 6 while G&O peak at 5, and for O 1 , CM6 peaks at degree 4 and G&O at 3. The power across most degrees are substantially larger for CM6 than for G&O for all three constituents. These amplitude differences may lie in the treatment of internal induced signals in the mantle or sources in the ionosphere, especially since G&O are based upon CHAOS model residual data. Recall that the CHAOS series does not model the nightside induced ionospheric field, which could leak into its residuals and subsequently be absorbed into the tidal fields. Conversely, if the mantle induction is not accurate in CM6, then this could also contaminate the tidal fields.
The constituents may also be compared in the spatial domain either in the form of amplitude and phase or as the real and imaginary parts of a complex phasor representation. To use the latter form, rewrite the tidal magnetic potential in Eq. 4 as where is a complex phasor whose real and imaginary parts may be examined. Maps of R B r M 2 and I B r M 2 are shown in Figs. 8, 9, and 10 for the M 2 , N 2 , and O 1 constituents, respectively, at a satellite altitude of 430 km and for the corresponding degree ranges shown in Fig. 7 from CM6 and G&O. A visual inspection of the M 2 maps shows very good agreement (17) Fig. 9 Real (left) and imaginary (right) parts of the radial component of the N 2 oceanic tidal magnetic field for SH degrees n = 1−12 at satellite altitude from CM6 (top) and G&O (bottom) in patterns, which is quantified by map correlation coefficients of 0.96 and 0.90 between R B r M 2 and between I B r M 2 , respectively, with CM6 having slightly higher amplitudes, as expected from the R n comparison. The importance of Swarm measurements in the extraction of M 2 is clear when one considers that the CIY4 and CM6 maps are likely trustworthy out to degree 28−36 , whereas that of CM5 are only to about degree 18. For the N 2 maps, the correlations between R B r M 2 and between I B r M 2 are 0.77 and 0.75, respectively, which indicates a modest agreement in patterns. The higher CM6 amplitudes are now very apparent. For the O 1 maps, the correlations between R B r M 2 and between I B r M 2 are 0.55 and 0.47, respectively, which is much worse than the agreement between semi-diurnal constituents, although many of the larger features appear to be centered in similar locations. However, the R B r M 2 map from CM6 shows contamination in the interiors of both Africa and Antarctica that are absent in the G&O map. Again, the CM6 power is much higher. It is interesting that in all constituents, the correlations between the real parts are consistently higher than between the imaginary parts. In general, however, all three constituents are in fair agreement and suggest that oceanic tidal signals can be detected and separated by satellite missions such as CHAMP and Swarm. Figure 11 shows the radial magnetic field at Earth's surface associated with the primary (top) and secondary (bottom) ionospheric current system during vernal equinox centered on local noon for various universal times. The behavior of the CM6 ionospheric field is in good agreement with previous CI models, such as CIY4. Note that the secondary induced field does not vanish during nighttime, but is rather broad-scaled and has odd equatorial symmetry, which is exactly what can map into the low-degree zonal coefficients of the core field in models that do not consider the ionosphere. Furthermore, since the induced ionosphere has power at solar periods of 24.0 and 12.0 h its treatment will no doubt have an effect on the recovery of the nearby periods of the diurnal and semi-diurnal tidal signals such as M 2 , N 2 , and O 1 and may be the reason for the discrepancies between the results from CM6 and G&O. difference between CM6 and CHAOS-7 q 0 1 (red) is also very good and extends over the entire CM6 time period. This is especially promising since CHAOS-7 does not estimate the magnetosphere in discrete time bins, but rather time variation is parameterized by the proxy RC index. Still there remains an unexplained annual variation in the difference. The difference of CM6 with Est (yellow) is much larger and erratic which may be attributable to the well-known baseline instabilities of the Dst index (e.g., Olsen et al. 2014). In addition, the calculation of Est is somewhat different than that used when assessing CIY4, which is due partly to the change of Dst from Quicklook to quasi-definitive and small changes in the conductivity model used to decompose Dst and RC into external and induced parts. Finally, the difference between CM6 and RC e (purple) is also reasonable, but the same annual variation is seen as with CHAOS-7.

IGRF-13 candidates
The International Association of Geomagnetism and Aeronomy (IAGA) working group V-MOD of Divison V requested candidate models for its 13th Generation of the International Geomagnetic Reference Field. These include a definitive model for epoch 2015.0 (DGRF2015) and an operational version for epoch 2020.0 (IGRF2020), both being static internal field SH expansions of degree n = 1−13 that represent the observational "main" magnetic field of the Earth at those epochs. Candidates for both the DGRF2015 and IGRF2020 models were provided by the lead at GSFC and were obtained from CM6. For DGRF2015, CM6 was simply evaluated at epoch 2015.0 and formal errors for the Gauss coefficients provided. The IGRF2020 was derived from a simple linear extrapolation passing through the CM6 Gauss coefficients at 2018.75 and 2019.0 and evaluated at epoch 2020.0, and likewise, the formal error-covariance was extrapolated to 2020.0. This 2018.75−2019.0 interval was chosen because it includes the latest extent of the Swarm data used in CM6 and provides a fair level of continuity in ġ 0 1 . The Swarm data allow for more accurate coefficients of degree n = 1−13 due to good global coverage, as opposed to later times that include only OHM data.

Conclusions
The CM6 model was derived from over 20 years of satellite and observatory ground-based data, essentially covering the "modern era" of satellite magnetometry. This is the first model to combine both CHAMP and Swarm data that are based upon the CI methodology. The resulting core, lithosphere, oceanic tidal M 2 , ionosphere, CM6 -MMA CM6 -CHAOS-7 CM6 + Est -12 nT CM6 + RC e -12 nT Fig. 12 The top panel shows 30 day averages of q 0 1 estimated from CM6 (blue), from CHAOS-7 (red), the Swarm MMA data product (yellow), −Est (purple), and −RC e (green) over the time domain of CM6. The bottom panel shows the differences between the CM6 q 0 1 and the other quantities magnetosphere, and associated induced magnetic field models are in very good agreement with those previously determined in other models. In particular, the core and SV spatial and temporal behavior are very similar to those of CHAOS-6 and CHAOS-7 with the largest discrepancies being well-understood in the context of the treatment of induced fields. Two new magnetic field sources have been included and correspond to the N 2 and O 1 oceanic tidal constituents. These are in generally fair to good agreement with the models obtained in Grayver and Olsen (2019) that were derived from CHAOS residuals and whose discrepancies are also likely related to the treatment of induced fields. In response to the IGRF-13 call, the GSFC candidate models for DGRF2015 and IGRF2020 were extracted from CM6 either directly, as in the former, or through simple linear extrapolation, as in the latter, along with their associated formal error, and submitted for consideration. CM6 will also serve as an anchor for extending CI models back in time to include other historical data types.