Reconstruction of Mercury's internal magnetic field beyond the octupole

The reconstruction of Mercury's internal magnetic field enables us to take a look into the inner heart of Mercury. In view of the BepiColombo mission, Mercury's magnetosphere is simulated using a hybrid plasma code, and the multipoles of the internal magnetic field are estimated from the virtual spacecraft data using three distinct reconstruction methods: the truncated singular value decomposition, the Tikhonov regularization and Capon's minimum variance projection. The study shows that a precise determination of Mercury's internal field beyond the octupole up to the dotriacontapole is possible and that Capon's method provides the same high performance as the Tikhonov regularization, which is superior to the performance of the truncated singular value decomposition.

1 Introduction

The in-depth analysis of planetary magnetic fields is a key element to understand the structure and dynamics of planetary interiors (Glassmeier and Heyner, 2021). In particular, revealing the nature of Mercury's markedly weak internal magnetic field is one of the primary goals of the BepiColombo mission (Benkhoff et al., 2010, 2021). Thereby, the detailed reconstruction of Mercury's internal field from the total measured magnetic field is of major importance for modeling Mercury's internal dynamo process. Especially the determination of higher orders of Mercury's internal multipole field reduces the degrees of freedom within the dynamo models (e.g., Heyner et al., 2021). Due to the plasma physical interaction with the solar wind, the magnetic field around Mercury is composed not only of the planetary internal field (such as the dynamo-generated field and the crustal remanent field), but also of the external component resulting from currents flowing within the magnetosphere (e.g., Glassmeier and Heyner, 2021; Wang et al., 2021). For the reconstruction of Mercury's internal magnetic field, each part of the field has to be parametrized properly.

Besides the parametrization of the magnetic field, the application of a robust inversion method for separating the magnetic field contributions is required. Several inversion methods have successfully been applied to the reconstruction of planetary magnetic fields. For example, the Tikhonov regularization (Tikhonov et al., 1995), also known as L² regularization, has been proposed for the analysis of Mercury's magnetic field. Katsura et al. (2021) separate the internal and external magnetic field contributions by making use of the Tikhonov regularization for estimating Mercury's inner core size. Wardinski et al. (2019, 2021) use a weighted Tikhonov regularization for reconstructing Mercury's internal magnetic field from MESSENGER data and also estimate the size of Mercury's core. On the other hand, Connerney (1981) and Connerney et al. (2018) successfully reconstructed Jupiter's magnetic field using a generalized inverse based on the singular value decomposition. Wang et al. (2021) also use the singular value decomposition for reconstructing Mercury's magnetospheric field. Furthermore, the Capon method (Capon, 1969) is currently being considered as a robust inversion method for Mercury's magnetic field analysis (Toepfer et al., 2020a, b). From the first point of view, it is unknown which inversion method provides more reliable results for the wanted model coefficients, and thus, the existing inversion methods have to be compared and calibrated. For example, Toepfer et al. (2020a) analyzed the performance of Capon's estimator in comparison with the least-squares fit estimator and showed that Capon's method provides more accurate results than the least-squares fit method.

As mentioned in Wardinski et al. (2021), former analyses of Mercury's magnetic field are based on potential field methods for the internal and external magnetic field contributions (Oliveira et al., 2015; Thébault et al., 2018; Wardinski et al., 2019) or on the analysis of magnetic equator crossings of the spacecraft (Anderson et al., 2012). For a detailed review of internal field models of Mercury, the reader is referred to the paper of Heyner et al. (2021). The analysis of the MESSENGER data provided an axisymmetric internal dipole field with a dipole moment of $-\mathrm{190}\mathrm{nT}{R}_{\mathrm{M}}^{\mathrm{3}}$, which is shifted northward by 0.2 R_M, where R_M=2440 km denotes the planetary radius of Mercury (Anderson et al., 2012). This field can equivalently be described as a superposition of multipole fields with the internal Gauss coefficients $g_{\mathrm{1}}^{\mathrm{0}}\approx -\mathrm{190}\mathrm{nT}$ for the dipole field, $g_{\mathrm{2}}^{\mathrm{0}}\approx -\mathrm{78}\mathrm{nT}$ for the quadrupole field, $g_{\mathrm{3}}^{\mathrm{0}}\approx -\mathrm{20}\mathrm{nT}$ for the octupole field, $g_{\mathrm{4}}^{\mathrm{0}}\approx -\mathrm{6}\mathrm{nT}$ for the hexadecapole field and $g_{\mathrm{5}}^{\mathrm{0}}\approx -\mathrm{2}\mathrm{nT}$ for the dotriacontapole field (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019). The resulting coefficients are slightly varying for different analyzed data sets of the whole MESSENGER mission (e.g., data from different altitude ranges). Plattner and Johnson (2021) and Wardinski et al. (2021) detected magnetic anomalies and nonaxisymmetric internal magnetic field contributions in the Northern Hemisphere. However, due to the geometry of the orbits, the MESSENGER mission only provided detailed information about the field in the Northern Hemisphere so that the reconstructed coefficients are correlated and covarying. The symmetrical distribution of the planned Mercury Planetary Orbiter (MPO) trajectories enables a more objective analysis of the internal magnetic field (e.g. Heyner et al., 2021).

In preparation for the analysis of the BepiColombo magnetic field data, the goal of the present study is the comparison of Capon's method with the truncated singular value decomposition and the Tikhonov regularization in application to simulated Mercury magnetic field data. Although the simulated data have to be regarded as a proxy for the MPO data that are not yet available, the great advantage of the application to simulated data lies in the fact that the exact solution of the inversion problem is known from the input of the simulation, which enables the comparison of the estimators with the ideal solution. The plasma interaction of Mercury with the solar wind is simulated using a hybrid plasma code. The resulting magnetic field data are parametrized by making use of a combination of the Gauss representation (Gauss, 1839; Glassmeier and Tsurutani, 2014) with the Mie representation (Backus, 1986; Backus et al., 1996; Olsen, 1997), called the Gauss–Mie representation, that has successfully been applied to reconstruct Mercury's internal magnetic field up to the octupole term (third degree in the multipole expansion) (Toepfer et al., 2021a). In the current study, the parametrization is extended up to Mercury's internal hexadecapole and dotriacontapole contributions, i.e., degrees 4 and 5. Afterwards, the truncated singular value decomposition and the Tikhonov regularization as well as Capon's method are compared, and the inversion methods are applied to the simulated magnetic field data for reconstructing Mercury's internal field up to the fifth degree.

2 The Gauss–Mie representation

The fluxgate magnetometer on board the Mercury Planetary Orbiter (MPO) (Glassmeier et al., 2010; Heyner et al., 2021) is built to measure the magnetic field B around Mercury along elliptic orbits. These orbits are conceptually covered by a spherical shell with inner radius a, outer radius c and mean radius $b=(a+c)/\mathrm{2}$. Since the shell in general includes current-carrying regions with ${\partial }_{\mathit{x}}\times \mathit{B}\ne \mathrm{0}$, the most commonly used Gauss representation (Gauss, 1839; Glassmeier and Tsurutani, 2014) does not yield a proper characterization of Mercury's magnetospheric magnetic field. By virtue of the superposition principle, the total measured field B can be decomposed into Bⁱ, B^e and B^sh. Bⁱ is the irrotational field contribution internally generated by currents flowing beneath the shell (r<a), and B^e is the irrotational field contribution externally generated by currents flowing above the shell (r>c). Briefly, we call Bⁱ the internal field and B^e the external field. B^sh is the rotational field generated by currents flowing within the shell. Considering the Mie representation (Backus, 1986; Backus et al., 1996; Olsen, 1997; Toepfer et al., 2021a), the currents flowing within the shell in the region $a<r<c$ generate toroidal ${\mathit{B}}_{\text{T}}^{\text{sh}}$ and poloidal ${\mathit{B}}_{\text{p}}^{\text{sh}}$ magnetic field structures that superpose the existing internal and external fields, delivering

$$\begin{array}{}\text{(1)}& \mathit{B}={\mathit{B}}^{\text{i}}+{\mathit{B}}^{\text{e}}+{\mathit{B}}_{\text{T}}^{\text{sh}}+{\mathit{B}}_{\text{p}}^{\text{sh}},\end{array}$$

which is called the Gauss–Mie representation of the magnetic field (Toepfer et al., 2021a).

Since both the internal and external fields Bⁱ and B^e are purely poloidal and irrotational, these parts can be parametrized via the Gauss representation. Therefore, there are scalar potentials Φⁱ and Φ^e with

$$\begin{array}{}\text{(2)}& {\mathit{B}}^{\text{i}}=-{\partial }_{\mathit{x}}{\mathrm{\Phi }}^{\text{i}}\end{array}$$

and

$$\begin{array}{}\text{(3)}& {\mathit{B}}^{\text{e}}=-{\partial }_{\mathit{x}}{\mathrm{\Phi }}^{\text{e}}.\end{array}$$

Using the body-fixed, planetary-centered spherical coordinates with radius $r\in \left[R_{\mathrm{M}},\mathrm{\infty }\right)$, azimuth angle $\mathit{\lambda }\in \left[\mathrm{0},\mathrm{2}\mathit{\pi }\right]$ and co-latitude $\mathit{\theta }\in \left[\mathrm{0},\mathit{\pi }\right]$, the scalar potentials can be expanded into spherical harmonics, yielding

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathrm{\Phi }}^{\text{i}}=R_{\mathrm{M}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\sum _{m=\mathrm{0}}^l{\left({\displaystyle \frac{R_{\mathrm{M}}}{r}}\right)}^{l+\mathrm{1}}\left[g_l^m\cos \left(m\mathit{\lambda }\right)+h_l^m\sin \left(m\mathit{\lambda }\right)\right]\\ & P_l^m\left(\cos \left(\mathit{\theta }\right)\right)\end{array}\end{array}$$

and

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathrm{\Phi }}^{\text{e}}=R_{\mathrm{M}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\sum _{m=\mathrm{0}}^l{\left({\displaystyle \frac{r}{R_{\mathrm{M}}}}\right)}^l\left[q_l^m\cos \left(m\mathit{\lambda }\right)+s_l^m\sin \left(m\mathit{\lambda }\right)\right]\\ & P_l^m\left(\cos \left(\mathit{\theta }\right)\right),\end{array}\end{array}$$

where R_M is the planetary radius of Mercury, and $P_l^m$ denotes the Schmidt-normalized associated Legendre polynomials of degree l and order m (Abramowitz and Stegun, 1972). Spherical harmonic expansion of the scalar potentials was introduced by the epoch-making work by Gauss (1839), also revisited by Glassmeier and Tsurutani (2014) with the contemporary English translation. The potential Φⁱ is determined by the internal Gauss coefficients $g_l^m$ and $h_l^m$, whereas the external Gauss coefficients $q_l^m$ and $s_l^m$ describe the external field contributions.

In general, the shell includes current-carrying regions, and therefore, the contributions ${\mathit{B}}_{\text{T}}^{\text{sh}}$ and ${\mathit{B}}_{\text{p}}^{\text{sh}}$ cannot be described by the Gauss representation. Making use of the Mie representation, one may introduce scalar functions ${\mathrm{\Psi }}_{\text{T}}^{\text{sh}}$ and ${\mathrm{\Psi }}_{\text{p}}^{\text{sh}}$ in the spirit of potential functions by defining the magnetic field through the curl of the scalar functions (multiplied by the position vector) as follows:

$$\begin{array}{}\text{(6)}& {\mathit{B}}_{\text{T}}^{\text{sh}}={\partial }_{\mathit{x}}\times \left({\mathrm{\Psi }}_{\text{T}}^{\text{sh}}\mathit{r}\right)\end{array}$$

and

$$\begin{array}{}\text{(7)}& {\mathit{B}}_{\text{p}}^{\text{sh}}={\partial }_{\mathit{x}}\times \left[{\partial }_{\mathit{x}}\times \left({\mathrm{\Psi }}_{\text{p}}^{\text{sh}}\mathit{r}\right)\right],\end{array}$$

where r=r e_r, and e_r is the unit vector in radial direction. Because of the underlying spherical geometry, it is straightforward to expand the scalar functions ${\mathrm{\Psi }}_{\text{T}}^{\text{sh}}$ and ${\mathrm{\Psi }}_{\text{p}}^{\text{sh}}$ into spherical harmonics. In contrast to the scalar potentials Φⁱ and Φ^e, the exact radial dependence of the corresponding expansion coefficients in the Mie representation is unknown. Thus, it is useful to perform a special Taylor series expansion for the coefficients with respect to the radius r in the vicinity of the mean radius b of the spherical shell (Toepfer et al., 2021a), providing

$$\begin{array}{}\text{(8)}& {\mathrm{\Psi }}_{\text{T}}^{\text{sh}}={\displaystyle \frac{R_{\mathrm{M}}}{r}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\sum _{m=\mathrm{0}}^l\left[{\mathit{\alpha }}_l^m+{\mathit{\alpha }}_l^{\prime m}\mathit{\rho }+\mathcal{O}\left({\mathit{\rho }}^{\mathrm{2}}\right)\right]P_l^m\left(\cos \left(\mathit{\theta }\right)\right)\end{array}$$

and

$$\begin{array}{}\text{(9)}& {\mathrm{\Psi }}_{\text{p}}^{\text{sh}}={\displaystyle \frac{R_{\mathrm{M}}^{\mathrm{2}}}{r}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\sum _{m=\mathrm{0}}^l\left[{\mathit{\beta }}_l^m+{\mathit{\beta }}_l^{\prime m}\mathit{\rho }+\mathcal{O}\left({\mathit{\rho }}^{\mathrm{2}}\right)\right]P_l^m\left(\cos \left(\mathit{\theta }\right)\right),\end{array}$$

where the variable $\mathit{\rho }=(r-b)/R_{\mathrm{M}}$ denotes the signed radial distance from the mean shell radius b in units of the planetary radius, and the Big O notation 𝒪 summarizes the orders of the Taylor series expansion higher than (or equal to) ρ². The expansion coefficients are given by

$$\begin{array}{ll}{\displaystyle }{\mathit{\alpha }}_l^m& {\displaystyle }=a_l^m\cos \left(m\mathit{\lambda }\right)+b_l^m\sin \left(m\mathit{\lambda }\right),\\ {\displaystyle }{\mathit{\alpha }}_l^{\prime m}& {\displaystyle }=a_l^{\prime m}\cos \left(m\mathit{\lambda }\right)+b_l^{\prime m}\sin \left(m\mathit{\lambda }\right),\\ {\displaystyle }{\mathit{\beta }}_l^m& {\displaystyle }=c_l^m\cos \left(m\mathit{\lambda }\right)+d_l^m\sin \left(m\mathit{\lambda }\right),\\ {\displaystyle }{\mathit{\beta }}_l^{\prime m}& {\displaystyle }=c_l^{\prime m}\cos \left(m\mathit{\lambda }\right)+d_l^{\prime m}\sin \left(m\mathit{\lambda }\right).\end{array}$$

Due to the geometry of the MPO orbits, the application of the thin shell approximation (Backus, 1986; Backus et al., 1996; Olsen, 1997; Toepfer et al., 2021a) is a valid assumption. Whereas the locally generated poloidal field ${\mathit{B}}_{\text{p}}^{\text{sh}}$ cannot be distinguished from the internal and external contributions Bⁱ and B^e within the reconstruction procedure (Toepfer et al., 2021a), the thin shell approximation enables us to separate the poloidal field into its internal and external parts. When the shell thickness is smaller than the length scale on which the currents change in radial direction, the poloidal field ${\mathit{B}}_{\text{p}}^{\text{sh}}$ generated by toroidal currents flowing within the shell can be neglected compared to the internal Bⁱ and external B^e poloidal fields. Combining the Gauss representation with the Mie representation by making use of the thin shell approximation, the magnetic field can be rewritten in the linear algebraic form

$$\begin{array}{}\text{(10)}& \mathit{B}=-{\partial }_{\mathit{x}}{\mathrm{\Phi }}^{\text{i}}-{\partial }_{\mathit{x}}{\mathrm{\Phi }}^{\text{e}}+{\partial }_{\mathit{x}}\times \left({\mathrm{\Psi }}_{\text{T}}^{\text{sh}}\mathit{r}\right)=\mathbf{H}\mathit{g},\end{array}$$

where the terms of the series expansions are arranged into the shape matrix H, which solely contains known information about the measurement positions. The corresponding expansion coefficients describing the amplitude of the field are summarized into the vector g.

It should be noted that the internal field is canonically described in a Mercury body-fixed (MBF) corotating coordinate system, whereas the external field is canonically described in a Mercury solar-orbital system (MSO), with the x axis oriented towards the sun, the z axis oriented parallel to the rotation axis, i.e., antiparallel to the internal dipole moment, and the y axis completing the right-handed system (Toepfer et al., 2021a; Heyner et al., 2021). In the present study, simulated Mercury magnetic field data are evaluated. Therefore, the internal and external fields are expressed in a Mercury anti-solar-orbital (MASO) coordinate system , where the x axis is oriented towards the nightside of Mercury (away from the sun), the z axis is oriented parallel to the rotation axis and the y axis completes the right-handed system (Toepfer et al., 2021a).

In general, the number of data points from an orbital mission is much larger than the number of wanted model coefficients. Thus, the resulting shape matrix H is rectangular, and its inverse H⁻¹ does not exist (Toepfer et al., 2020b). Therefore, Eq. (10) describes an overdetermined system of linear equations. A unique solution for the desired model coefficients g is not available, so the coefficients have to be estimated from the measurements by making use of a suitable inversion method.

3 Inversion methods

For the reconstruction of Mercury's internal magnetic field, various kinds of data inversion techniques are available, such as the least-squares fit method, the singular value decomposition, the Tikhonov regularization and Capon's minimum variance method (e.g., Haykin, 2014). The construction of these inversion techniques is reviewed along with merits and demerits in this section.

3.1 Least-squares fit

The most prominent inversion method for linear inversion problems is the least-squares fit method (LSF) (e.g., Haykin, 2014). The method minimizes the quadratic deviation

$$\begin{array}{}\text{(11)}& \mathrm{minimize}|\mathbf{H}\mathit{g}-\mathit{B}{|}^{\mathrm{2}}\end{array}$$

between the model H g and the measurements B with respect to the set of model parameters g with unknown values, resulting in the least-squares fit estimator

$$\begin{array}{}\text{(12)}& {\mathit{g}}_{\text{L}}={\left[{\mathbf{H}}^{†}\mathbf{H}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}\mathit{B},\end{array}$$

where the dagger † denotes the Hermitian conjugate. Thus, the LSF method solely weights the data by the shape matrix H. When the measurements are completely described by the underlying model, the LSF estimator provides the most accurate estimation for the wanted model coefficients (Toepfer et al., 2020a).

3.2 Singular value decomposition

The singular value decomposition (SVD) generalizes the LSF method (e.g., Haykin, 2014). The method is based on the decomposition of the shape matrix $\mathbf{H}\in {\mathbb{R}}^{m\times n}$ with rank n into the form

$$\begin{array}{}\text{(13)}& \mathbf{H}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{†},\end{array}$$

where $\mathbf{U}\in {\mathbb{R}}^{m\times m}$ and $\mathbf{V}\in {\mathbb{R}}^{n\times n}$ are orthonormal transformations, so that ${\mathbf{U}}^{†}={\mathbf{U}}^{-\mathrm{1}}$ as well as ${\mathbf{V}}^{†}={\mathbf{V}}^{-\mathrm{1}}$ is valid (e.g., Connerney, 1981; Haykin, 2014; Connerney et al., 2018). The matrix

$$\begin{array}{}\text{(14)}& \mathbf{\Sigma }=\left[\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{1}}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\dots }& {\mathrm{\Sigma }}_n\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}\end{array}\right]\in {\mathbb{R}}^{m\times n},\end{array}$$

where ${\mathrm{\Sigma }}_{\mathrm{1}}\ge \mathrm{\dots }\ge {\mathrm{\Sigma }}_n\ge \mathrm{0}$ (sorted in descending order), contains the so-called singular values Σ_i of the shape matrix H, which are defined as the square roots of the eigenvalues of the matrix H^†H. Inserting the singular value decomposition of the matrix H (Eq. 13) into Eq. (10) delivers the singular value estimator

$$\begin{array}{}\text{(15)}& {\mathit{g}}_{\text{S}}={\mathbf{H}}^{+}\mathit{B}=\mathbf{V}{\mathbf{\Sigma }}^{+}{\mathbf{U}}^{†}\mathit{B}=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{{\mathrm{\Sigma }}_i}}{\mathit{v}}_i{\mathit{u}}_i^{†}\mathit{B}\end{array}$$

for the wanted model coefficients g, where the vectors u_i and v_i are the corresponding column vectors of the matrices U and V, respectively (e.g., Haykin, 2014). The matrix

$$\begin{array}{}\text{(16)}& {\mathbf{H}}^{+}=\mathbf{V}{\mathbf{\Sigma }}^{+}{\mathbf{U}}^{†}=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{{\mathrm{\Sigma }}_i}}{\mathit{v}}_i{\mathit{u}}_i^{†}\end{array}$$

is called the “pseudo-inverse” or generalized inverse of the matrix H, where the matrix Σ⁺ is given by

$$\begin{array}{}\text{(17)}& {\left({\mathbf{\Sigma }}^{+}\right)}_{ij}=\left\{\begin{array}{l}\frac{\mathrm{1}}{{\mathrm{\Sigma }}_i},i=j\\ \mathrm{0},\mathrm{else}.\end{array}\right.\end{array}$$

In the case of a full rank, the matrix H^†H is invertible, and thus, the estimator resulting from the singular value decomposition restores the least-squares fit estimator

$$\begin{array}{}\text{(18)}& {\mathit{g}}_{\text{S}}={\mathbf{H}}^{+}\mathit{B}={\left[{\mathbf{H}}^{†}\mathbf{H}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}\mathit{B}={\mathit{g}}_{\text{L}}.\end{array}$$

Conferring to Eq. (15), the solution g_S (as well as the solution g_L) is determined by the inverse of singular values. In the case of small or even vanishing singular values, i.e., rank(H)<n, the solution diverges. Furthermore, the condition number

$$\begin{array}{}\text{(19)}& \mathit{\kappa }\left(\mathbf{H}\right)={\displaystyle \frac{\underset{i}{max}{\mathrm{\Sigma }}_i}{\underset{i}{min}{\mathrm{\Sigma }}_i}}={\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{1}}}{{\mathrm{\Sigma }}_n}}\end{array}$$

of the shape matrix increases κ(H)≫1 for small singular values, and thus, the inversion problem (Eq. 10) is said to be ill-posed. A large condition number impairs the solvability of the inversion problem and furthermore increases the error propagation (Toepfer et al., 2021b).

One of the most commonly used techniques for reducing the condition number is the low-rank approximation or truncated singular value decomposition (TSVD) (Eckart and Young, 1936). Within this approximation, only k singular values, where k<n, are considered within the solution, and therefore, the shape matrix H is approximated by a shape matrix

$$\begin{array}{}\text{(20)}& {\mathbf{H}}_k={\mathbf{U}}_k{\mathbf{\Sigma }}_k{\mathbf{V}}_k^{†},\end{array}$$

where the matrices U_k and V_k are composed of the first k column vectors of the matrices U and V, respectively. The diagonal elements of the matrix Σ_k are given by the largest k singular values of the shape matrix H. Therefore, the matrix H_k has a lower rank, i.e., rank(H_k)=k, and a lower condition number

$$\begin{array}{}\text{(21)}& \mathit{\kappa }\left({\mathbf{H}}_k\right)={\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{1}}}{{\mathrm{\Sigma }}_k}}<{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{1}}}{{\mathrm{\Sigma }}_n}}=\mathit{\kappa }\left(\mathbf{H}\right)\end{array}$$

is achieved. Due to the modification of the shape matrix, the solution g_S transfers onto

$$\begin{array}{}\text{(22)}& {\mathit{g}}_{\text{TSVD}}={\mathbf{H}}_k^{+}\mathit{B}=\sum _{i=\mathrm{1}}^k{\displaystyle \frac{\mathrm{1}}{{\mathrm{\Sigma }}_i}}{\mathit{v}}_i{\mathit{u}}_i^{†}\mathit{B},\end{array}$$

which will be called the TSVD estimator.

On the one hand, the decreased condition number improves the solvability of the inversion problem. On the other hand, it should be noted that the elimination of singular values causes a lack of information, which can decrease the performance of the data analysis. This lack of information can be quantitatively estimated via the model resolution matrix (Menke, 2012; Heyner et al., 2021)

$$\begin{array}{}\text{(23)}& \mathbf{R}={\mathbf{H}}_k^{+}\mathbf{H}.\end{array}$$

Because of g_TSVD=R g, where g denotes the ideal (imaginary) solution for the model coefficients, the model resolution matrix identifies the resolution of each model coefficient. In the case of R≠I, the estimated coefficients are determined by a linear combination of the ideal solution, and thus, there are model covariances. If R=I, each coefficient is resolved for 100 % (Menke, 2012). To achieve the highest performance of the estimator, the condition number and the lack of information have to come to a compromise (Connerney, 1981; Toepfer et al., 2021a).

3.3 Tikhonov regularization

Within the application of the LSF method and the singular value decomposition, the wanted model coefficients are required to satisfy the highly overdetermined system of linear equations B=H g. From a physical point of view, an additional constraint for the solution (for example, minimum energy) can be incorporated to reduce the degrees of freedom within the estimation.

Considering the analysis of Mercury's internal magnetic field, the wanted model coefficients describe the amplitude of the magnetic field. Thus, it is obvious that the currents flowing within the magnetosphere as well as inside of Mercury generate magnetic fields with a minimal energy. The energy spectrum W_l of each spherical harmonic degree l of the internal magnetic field at Mercury's surface is represented by the Mauersberger–Lowes spectrum (Mauersberger, 1956; Lowes, 1966)

$$\begin{array}{}\text{(24)}& W_l=(l+\mathrm{1})\sum _{m=\mathrm{0}}^l\left[(g_l^m{)}^{\mathrm{2}}+(h_l^m{)}^{\mathrm{2}}\right].\end{array}$$

Therefore, minimizing the energy spectrum W_l is equivalent to the minimization of the norm of the wanted model coefficients $|\mathit{g}{|}^{\mathrm{2}}$, and thus it is useful to search for solutions having the minimal norm. This constraint is known as the Tikhonov regularization, L² regularization or minimum norm solution (e.g., Tikhonov et al., 1995; Haykin, 2014). Using this additional constraint, the original least-squares fit problem (see Eq. 11) can be extended as

$$\begin{array}{}\text{(25)}& \mathrm{minimize}|\mathbf{H}\mathit{g}-\mathit{B}{|}^{\mathrm{2}}+\mathit{\alpha }|\mathit{g}{|}^{\mathrm{2}},\end{array}$$

where α is the so-called regularization parameter, which describes the corresponding Lagrange multiplier of the additional constraint. The Tikhonov estimator for the wanted model coefficients results in

$$\begin{array}{}\text{(26)}& {\mathit{g}}_{\text{T}}={\left[{\mathbf{H}}^{†}\mathbf{H}+\mathit{\alpha }\mathbf{I}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}\mathit{B},\end{array}$$

where I denotes the identity matrix. Thus, the additional constraint modifies the shape matrix and, therefore, modifies the location dependence of the measurement points with respect to the underlying model.

Inserting the singular value decomposition of the matrix H into Eq. (26) delivers

$$\begin{array}{ll}{\displaystyle }{\mathit{g}}_{\text{T}}& {\displaystyle }={\left[{\mathbf{H}}^{†}\mathbf{H}+\mathit{\alpha }\mathbf{I}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}\mathit{B}\\ {\displaystyle }& {\displaystyle }={\left[{\left(\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{†}\right)}^{†}\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{†}+\mathit{\alpha }\mathbf{I}\right]}^{-\mathrm{1}}{\left(\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{†}\right)}^{†}\mathit{B}\\ {\displaystyle }& {\displaystyle }={\left[\mathbf{V}{\mathbf{\Sigma }}^{†}{\mathbf{U}}^{†}\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{†}+\mathit{\alpha }\mathbf{I}\right]}^{-\mathrm{1}}\mathbf{V}{\mathbf{\Sigma }}^{†}{\mathbf{U}}^{†}\mathit{B}\\ {\displaystyle }& {\displaystyle }={\left[\mathbf{V}{\mathbf{\Sigma }}^{†}\mathbf{\Sigma }{\mathbf{V}}^{†}+\mathit{\alpha }\mathbf{V}{\mathbf{V}}^{†}\right]}^{-\mathrm{1}}\mathbf{V}{\mathbf{\Sigma }}^{†}{\mathbf{U}}^{†}\mathit{B}\\ {\displaystyle }& {\displaystyle }={\left[\mathbf{V}\left({\mathbf{\Sigma }}^{†}\mathbf{\Sigma }+\mathit{\alpha }\mathbf{I}\right){\mathbf{V}}^{†}\right]}^{-\mathrm{1}}\mathbf{V}{\mathbf{\Sigma }}^{†}{\mathbf{U}}^{†}\mathit{B}\\ {\displaystyle }& {\displaystyle }=\mathbf{V}{\left[{\mathbf{\Sigma }}^{†}\mathbf{\Sigma }+\mathit{\alpha }\mathbf{I}\right]}^{-\mathrm{1}}{\mathbf{V}}^{†}\mathbf{V}{\mathbf{\Sigma }}^{†}{\mathbf{U}}^{†}\mathit{B}\\ {\displaystyle }& {\displaystyle }=\mathbf{V}{\left[{\mathbf{\Sigma }}^{†}\mathbf{\Sigma }+\mathit{\alpha }\mathbf{I}\right]}^{-\mathrm{1}}{\mathbf{\Sigma }}^{†}{\mathbf{U}}^{†}\mathit{B}\\ {\displaystyle }& {\displaystyle }=\mathbf{V}{\mathbf{\Sigma }}_{\text{T}}^{+}{\mathbf{U}}^{†}\mathit{B}\\ \text{(27)}& {\displaystyle }& {\displaystyle }={\mathbf{H}}_{\text{T}}^{+}\mathit{B},\end{array}$$

where the matrix

$$\begin{array}{}\text{(28)}& {\mathbf{H}}_{\text{T}}^{+}=\mathbf{V}{\mathbf{\Sigma }}_{\text{T}}^{+}{\mathbf{U}}^{†}\end{array}$$

will be called the Tikhonov inverse and

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& {\mathbf{\Sigma }}_{\text{T}}^{+}={\left[{\mathbf{\Sigma }}^{†}\mathbf{\Sigma }+\mathit{\alpha }\mathbf{I}\right]}^{-\mathrm{1}}{\mathbf{\Sigma }}^{†}\\ & =\left[\begin{array}{ccccccc}\frac{{\mathrm{\Sigma }}_{\mathrm{1}}}{{\mathrm{\Sigma }}_{\mathrm{1}}^{\mathrm{2}}+\mathit{\alpha }}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \frac{{\mathrm{\Sigma }}_{\mathrm{2}}}{{\mathrm{\Sigma }}_{\mathrm{2}}^{\mathrm{2}}+\mathit{\alpha }}& \mathrm{\ddots }& \mathrm{⋮}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{0}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \frac{{\mathrm{\Sigma }}_n}{{\mathrm{\Sigma }}_n^{\mathrm{2}}+\mathit{\alpha }}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\end{array}\right]\in {\mathbb{R}}^{n\times m}.\end{array}\end{array}$$

Thus, the original shape matrix H transfers onto the shape matrix

$$\begin{array}{}\text{(30)}& {\mathbf{H}}_{\text{T}}=\mathbf{U}{\mathbf{\Sigma }}_{\text{T}}{\mathbf{V}}^{†},\end{array}$$

where

$$\begin{array}{}\text{(31)}& {\mathbf{\Sigma }}_{\text{T}}=\left[\begin{array}{cccc}\frac{{\mathrm{\Sigma }}_{\mathrm{1}}^{\mathrm{2}}+\mathit{\alpha }}{{\mathrm{\Sigma }}_{\mathrm{1}}}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& \frac{{\mathrm{\Sigma }}_{\mathrm{2}}^{\mathrm{2}}+\mathit{\alpha }}{{\mathrm{\Sigma }}_{\mathrm{2}}}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\ddots }& \mathrm{0}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \frac{{\mathrm{\Sigma }}_n^{\mathrm{2}}+\mathit{\alpha }}{{\mathrm{\Sigma }}_n}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}\end{array}\right]\in {\mathbb{R}}^{m\times n},\end{array}$$

which has a lower condition number

$$\begin{array}{ll}{\displaystyle }\mathit{\kappa }\left({\mathbf{H}}_{\text{T}}\right)=\mathit{\kappa }\left({\mathbf{H}}_{\text{T}}^{+}\right)& {\displaystyle }={\displaystyle \frac{\underset{i}{max}\frac{{\mathrm{\Sigma }}_i^{\mathrm{2}}+\mathit{\alpha }}{{\mathrm{\Sigma }}_i}}{\underset{i}{min}\frac{{\mathrm{\Sigma }}_i^{\mathrm{2}}+\mathit{\alpha }}{{\mathrm{\Sigma }}_i}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\underset{i}{max}\left({\mathrm{\Sigma }}_i+\frac{\mathit{\alpha }}{{\mathrm{\Sigma }}_i}\right)}{\underset{i}{min}\left({\mathrm{\Sigma }}_i+\frac{\mathit{\alpha }}{{\mathrm{\Sigma }}_i}\right)}}\\ {\displaystyle }& {\displaystyle }\le {\displaystyle \frac{\underset{i}{max}{\mathrm{\Sigma }}_i+\mathit{\alpha }\underset{i}{max}\frac{\mathrm{1}}{{\mathrm{\Sigma }}_i}}{\underset{i}{min}{\mathrm{\Sigma }}_i+\mathit{\alpha }\underset{i}{min}\frac{\mathrm{1}}{{\mathrm{\Sigma }}_i}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{1}}+\frac{\mathit{\alpha }}{{\mathrm{\Sigma }}_n}}{{\mathrm{\Sigma }}_n+\frac{\mathit{\alpha }}{{\mathrm{\Sigma }}_{\mathrm{1}}}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{1}}{\mathrm{\Sigma }}_n+\mathit{\alpha }}{{\mathrm{\Sigma }}_n}}{\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{1}}}{{\mathrm{\Sigma }}_{\mathrm{1}}{\mathrm{\Sigma }}_n+\mathit{\alpha }}}\\ \text{(32)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathrm{\Sigma }}_{\mathrm{1}}}{{\mathrm{\Sigma }}_n}}=\mathit{\kappa }\left(\mathbf{H}\right).\end{array}$$

Therefore, the modification of the shape matrix improves the solvability of the inversion problem in analogy to the TSVD.

Comparison of the SVD estimator

$$\begin{array}{}\text{(33)}& {\mathit{g}}_{\text{S}}=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{\mathrm{1}}{{\mathrm{\Sigma }}_i}}{\mathit{v}}_i{\mathit{u}}_i^{†}\mathit{B}\end{array}$$

with the Tikhonov estimator

$$\begin{array}{}\text{(34)}& {\mathit{g}}_{\text{T}}={\mathbf{H}}_{\text{T}}^{+}\mathit{B}=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{{\mathrm{\Sigma }}_i}{{\mathrm{\Sigma }}_i^{\mathrm{2}}+\mathit{\alpha }}}{\mathit{v}}_i{\mathit{u}}_i^{†}\mathit{B}\end{array}$$

shows that the original singular values are modified by α. In the case of small or even vanishing singular values, the solution g_T remains finite. Furthermore, within the Tikhonov solution, n (modified) singular values are considered, whereas within the TSVD solution (Eq. 22), only k<n singular values are incorporated. Making use of the TSVD, the information that is included within the truncated singular values is eliminated. Within the application of the Tikhonov regularization, small singular values are not truncated but shifted. Thus, the information which is contained in the shifted singular values is still present in a modified form. In the case of α=0, the Tikhonov estimator transfers onto the LSF estimator.

In analogy to the TSVD, the modification of the shape matrix results in a nontrivial resolution matrix

$$\begin{array}{}\text{(35)}& {\mathbf{R}}_{\text{T}}={\mathbf{H}}_{\text{T}}^{+}\mathbf{H}\ne \mathbf{I},\end{array}$$

and thus, in general there are model covariances.

3.4 Capon's method

Capon's minimum variance projection is broadly established in the analysis of seismic and plasma waves (Capon, 1969; Motschmann et al., 1996; Narita, 2012). In view of the BepiColombo mission, the method is currently being considered as a robust inversion method for the analysis of Mercury's internal magnetic field (Toepfer et al., 2020a, b, 2021a).

Due to the complexity of Mercury's magnetosphere, the entire parametrization of the magnetic field contributions, generated by currents flowing within the magnetosphere, is unrealizable. Thus, it is useful to decompose the magnetic field B into parametrized parts H g (see Eq. 10) and non-parametrized parts v as well as measurement noise n, so that

$$\begin{array}{}\text{(36)}& \mathit{B}=\mathbf{H}\mathit{g}+\mathit{v}+\mathit{n}\end{array}$$

is valid. The measurement noise is assumed to be Gaussian with variance ${\mathit{\sigma }}_n^{\mathrm{2}}$ and zero mean (〈n〉=0). Since the shape matrix H is not invertible and the non-parametrized parts are unknown, the exact solution for the wanted model coefficients g is not available in general. Capon's method delivers an estimator g_C for the ideal solution g. The method is based on the construction of a filter matrix w, that minimizes the output power

$$\begin{array}{}\text{(37)}& P=\mathrm{tr}\left[{\mathbf{w}}^{†}\mathbf{M}\mathbf{w}\right]\end{array}$$

with respect to the distortionless constraint (also referred to as the unit gain constraint)

$$\begin{array}{}\text{(38)}& {\mathbf{w}}^{†}\mathbf{H}=\mathbf{I},\end{array}$$

where tr [w^†M w] is the trace of the matrix w^†M w, and I is the identity matrix. The matrix $\mathbf{M}=\langle \mathit{B}\circ \mathit{B}\rangle$ denotes the data covariance matrix. Capon's estimator realizing the minimal output power results in

$$\begin{array}{}\text{(39)}& {\mathit{g}}_{\text{C}}={\mathbf{w}}^{†}\langle \mathit{B}\rangle ={\left[{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}\mathbf{H}\right]}^{-\mathrm{1}}{\mathbf{H}}^{†}{\mathbf{M}}^{-\mathrm{1}}\langle \mathit{B}\rangle .\end{array}$$

The comparison of Capon's estimator with the LSF estimator (see Eq. 12) shows that Capon's method can be interpreted as a special case of the weighted least-squares fit method (Toepfer et al., 2020b). The robustness of the method can be improved by the diagonal loading technique $\mathbf{M}\to \mathbf{M}+{\mathit{\sigma }}_d^{\mathrm{2}}\mathbf{I}$, where σ_d is the diagonal loading parameter (Toepfer et al., 2020b). Inserting the diagonally loaded data covariance matrix $\mathbf{M}=\langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle +{\mathit{\sigma }}^{\mathrm{2}}\mathbf{I}$, where ${\mathit{\sigma }}^{\mathrm{2}}={\mathit{\sigma }}_n^{\mathrm{2}}+{\mathit{\sigma }}_d^{\mathrm{2}}$, into the output power and making use of ${\mathit{g}}_{\text{C}}={\mathbf{w}}^{†}\langle \mathit{B}\rangle$ delivers

$$\begin{array}{ll}{\displaystyle }P& {\displaystyle }=\mathrm{tr}\left[{\mathbf{w}}^{†}\mathbf{M}\mathbf{w}\right]=\mathrm{tr}\left[{\mathbf{w}}^{†}\langle \mathit{B}\rangle \circ \langle \mathit{B}\rangle \mathbf{w}\right]+{\mathit{\sigma }}^{\mathrm{2}}\mathrm{tr}\left[{\mathbf{w}}^{†}\mathbf{w}\right]\\ {\displaystyle }& {\displaystyle }=\mathrm{tr}\left[{\mathit{g}}_{\text{C}}\circ {\mathit{g}}_{\text{C}}\right]+{\mathit{\sigma }}^{\mathrm{2}}\mathrm{tr}\left[{\mathbf{w}}^{†}\mathbf{w}\right]\\ \text{(40)}& {\displaystyle }& {\displaystyle }=|{\mathit{g}}_{\text{C}}{|}^{\mathrm{2}}+{\mathit{\sigma }}^{\mathrm{2}}\mathrm{tr}\left[{\mathbf{w}}^{†}\mathbf{w}\right].\end{array}$$

Thus, Capon's method minimizes the energy spectrum in analogy to the Tikhonov regularization. In contrast to the Tikhonov regularization, Capon's method additionally weights the data by the inverse data covariance matrix (see Eq. 39).

In the case of small singular values of the shape matrix H, the performance of Capon's method decreases in analogy to the LSF method. Since Capon's method already minimizes the energy spectrum P, the method cannot be improved by making use of the Tikhonov regularization, and thus, the TSVD should be applied to reduce the condition number. Within many applications it is unknown which singular values have to be considered within the solution and which singular values can be dropped. On the other hand, the application of the Tikhonov regularization modifies the condition number of the original shape matrix (κ(H_T)≤κ(H)), and thus, the Tikhonov regularization delivers the “optimal” condition number for solving the problem. Therefore, the number of singular values that have to be considered within the solution can be estimated by choosing k, so that

$$\begin{array}{}\text{(41)}& \mathit{\kappa }\left({\mathbf{H}}_k\right)\approx \mathit{\kappa }\left({\mathbf{H}}_{\text{T}}\right)\end{array}$$

is valid. Due to the application of the TSVD, Capon's filter matrix is modified to

$$\begin{array}{}\text{(42)}& {\mathbf{w}}_k^{†}={\left[{\mathbf{H}}_k^{†}{\mathbf{M}}^{-\mathrm{1}}{\mathbf{H}}_k\right]}^{-\mathrm{1}}{\mathbf{H}}_k^{†}{\mathbf{M}}^{-\mathrm{1}},\end{array}$$

on the cost of violating the unit gain condition (Toepfer et al., 2021a), i.e.,

$$\begin{array}{}\text{(43)}& {\mathbf{R}}_{\text{C}}={\mathbf{w}}_k^{†}\mathbf{H}\ne \mathbf{I}.\end{array}$$

4 Application to simulated Mercury magnetic field data

In the following, the above-presented inversion methods are applied to simulated Mercury magnetic field data for reconstructing Mercury's internal multipole field up to the dotriacontapole term. The internal field is modeled as being generated by the dynamo field (crustal fields are not considered here), and the multipole spectrum model is taken from the MESSENGER results (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019). The magnetic field is sampled along the trajectories of virtual spacecraft (representing the MPO spacecraft), and the multipole spectrum is estimated using the different inversion techniques for the virtual spacecraft data. The estimators for the spectrum are compared to the internal Gauss coefficients implemented in the simulation.

4.1 Hybrid simulation of Mercury's magnetosphere

The magnetic field resulting from the plasma interaction of Mercury with the solar wind is simulated with the hybrid code AIKEF (Müller et al., 2011), that has successfully been applied to several problems in Mercury's plasma interaction (e.g., Müller et al., 2011; Exner et al., 2018, 2020). The internal Gauss coefficients $g_{\mathrm{1}}^{\mathrm{0}}=-\mathrm{190}\mathrm{nT}$ (dipole field), $g_{\mathrm{2}}^{\mathrm{0}}=-\mathrm{78}\mathrm{nT}$ (quadrupole field), $g_{\mathrm{3}}^{\mathrm{0}}=-\mathrm{20}\mathrm{nT}$ (octupole field), $g_{\mathrm{4}}^{\mathrm{0}}=-\mathrm{6}\mathrm{nT}$ (hexadecapole field) and $g_{\mathrm{5}}^{\mathrm{0}}=\mathrm{8}\mathrm{nT}$ (dotriacontapole field) (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019) are implemented in the simulation code. For the first validation, the value $g_{\mathrm{5}}^{\mathrm{0}}=\mathrm{8}\mathrm{nT}$ was chosen from a synthetic Mercury magnetic field model of Thébault et al. (2018). The interplanetary magnetic field with a magnitude of B_IMF=20 nT is oriented along the vector $(x,y,z{)}^{\text{T}}=(\mathrm{0},\mathrm{0},\mathrm{1}{)}^{\text{T}}$ in the Mercury anti-solar-orbital (MASO) frame. The magnetic field data are simulated under stationary solar wind conditions, where the solar wind velocity of v_sw=400 km s⁻¹ points along the x axis (Milillo et al., 2020). The solar wind proton density number is chosen to $n_{\mathrm{sw}}=\mathrm{30}{\mathrm{cm}}^{-\mathrm{3}}$, and the electron T_e and proton temperatures T_p are implemented to $T_{\text{e}}=T_{\text{p}}=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{K}$ (Exner et al., 2020; Milillo et al., 2020). The resulting magnitude of the magnetic field in the x−z plane is displayed in Fig. 1.

Figure 1 Simulated magnitude of the magnetic field B in the x−z plane, where internal multipoles are implemented from the dipole to the dotriacontapole field. The white lines represent the magnetic field lines, and the grey circle of radius 1 R_M symbolizes Mercury.

The geometry of Mercury's magnetosphere is mainly dominated by the internal dipole field. Also the internal quadrupole field in terms of the apparently northward shifted dipole field is visible. The influence of the octupole, hexadecapole and dotriacontapole fields is not visually noticeable within the figure.

4.2 Reconstruction of the multipole coefficients

The internal Gauss coefficients implemented in the simulation code represent the ideal solution g for the data inversion. For the reconstruction of the coefficients, the magnetic field data are evaluated along meridional elliptic orbits around Mercury, representing the planned MPO orbits. The orbital plane is rotated about the rotation axis (z axis) from $-\mathrm{50}{}^{\circ }$ (afternoon/post-midnight sector) over 0^∘ (noon/midnight, x−z plane) to 50^∘ (morning/pre-midnight sector). Figure 2 displays the distribution of the synthetically generated measurement positions in the x−z plane. The internal and external potentials (Φⁱ, Φ^e) are expanded into spherical harmonics up to the fifth degree, representing the internal/external dipole, quadrupole, octupole, hexadecapole and dotriacontapole field. The terms of the degrees l≤4 are expanded up to the order m=l. The dotriacontapole (l=5) is only expanded up to the zeroth order (m=0) for simplicity. The scalar function ${\mathrm{\Psi }}_{\text{T}}^{\text{sh}}$ of the toroidal magnetic field ${\mathit{B}}_{\text{T}}^{\text{sh}}$ is expanded into spherical harmonics up to the second degree and order and additionally into a first-order Taylor series for the radial distance. The scalar function ${\mathrm{\Psi }}_{\text{p}}^{\text{sh}}$ of the poloidal magnetic field ${\mathit{B}}_{\text{p}}^{\text{sh}}$ is neglected by making use of the thin shell approximation. Thus, the magnetic field is described by 66 expansion coefficients with 25 internal Gauss coefficients, 25 external Gauss coefficients and 16 toroidal coefficients. These coefficients are to be determined from the data by the inversion techniques. Since the exact solution for the internal Gauss coefficients is known a priori as the inputs to the simulation, in the following discussion we will focus on the reconstructed 25 internal coefficients.

Figure 2 Synthetically generated measurement positions (red dots) in the x−z plane ($\mathit{\lambda }=\mathrm{0}{}^{\circ }$). The orbital plane is rotated about the rotation axis (z axis) from $-\mathrm{50}{}^{\circ }$ over 0 to 50^∘. The grey circle of radius 1 R_M symbolizes Mercury.

The optimal regularization parameter for the application of the Tikhonov regularization results in α≈0.9, so that κ(H_T)≈86. Incorporating all the 66 singular values of the shape matrix results in a condition number of κ(H)≈540 for the original shape matrix. Considering 60 singular values within the reconstruction procedure, the condition number decreases to κ(H₆₀)≈90. Thus, for the calculation of Capon's estimator and for the TSVD estimator, 60 singular values are incorporated. The optimal diagonal loading parameter for the application of Capon's method results in σ=590 nT. The optimal regularization parameter α as well as the optimal diagonal loading parameter σ can be determined by minimizing the deviation between the corresponding estimator and the ideal solution, which is known from the input of the simulation. Within the practical application of the methods in future satellite experiments, the exact solution is unknown, and thus, the regularization parameter as well as the diagonal loading parameter can be estimated by making use of the L-curve technique (Toepfer et al., 2020b). The estimators resulting from the TSVD method, Capon's method and the Tikhonov regularization are displayed in Table 1.

Table 1 Implemented and reconstructed internal Gauss coefficients for the dipole, quadrupole, octupole, hexadecapole and dotriacontapole field. The implemented value $g_{\mathrm{5}}^{\mathrm{0}}$ of the dotriacontapole field was chosen from a synthetic Mercury magnetic field model of Thébault et al. (2018). The errors are derived by synthetically disturbing the simulated data and measurement positions as discussed in the text.

The deviations between the reconstructed and the implemented coefficients result in $|{\mathit{g}}_{\text{TSVD}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{3.9}$ %, $|{\mathit{g}}_{\text{C}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{2.6}$ % and $|{\mathit{g}}_{\text{T}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{2.6}$ %, so that the reconstructed and the implemented coefficients are in good agreement. Since the TSVD method only weights the data by the position information, the TSVD estimator yields the largest deviation. The deviation of Capon's estimator equals the deviation of the Tikhonov estimator, since both methods incorporate the constraint of a minimum norm solution. In particular, Mercury's internal hexadecapole and dotriacontapole fields can be reconstructed from the data to a good precision. Furthermore, the model resolution matrices of Capon's estimator and the Tikhonov regularization are of the same order and close to the identity matrix

$$\begin{array}{}\text{(44)}& {\mathbf{R}}_{\text{T}}\approx {\mathbf{R}}_{\text{C}}\approx \mathbf{I},\end{array}$$

providing a high model resolution. Thus, the eliminated singular values do not contain any additional information that may improve the inversion. Taking into account 66 singular values within the calculation of Capon's estimator and the TSVD estimator, the TSVD estimator transfers onto the LSF estimator, so that ${\mathbf{R}}_{\text{C}}=\mathbf{R}=\mathbf{I}$. In this case, the deviation between Capon's estimator and the ideal solution results in $|{\mathit{g}}_{\text{C}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{3.3}$ % for σ=420 nT and $|{\mathit{g}}_{\text{L}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{6.2}$ % for the LSF estimator. Thus, the errors can be reduced by making use of the TSVD for the shape matrix. Furthermore, it should be noted that the application of the Gauss–Mie representation improves the inversion results significantly in comparison to the sole parametrization of the field via the Gauss representation (Toepfer et al., 2021a). Additionally, the reconstructed Gauss coefficients of the external irrotational field are in agreement with the values reconstructed in former works (Wardinski et al., 2019; Toepfer et al., 2021a).

Table 2 Implemented and reconstructed internal Gauss coefficients for the dipole, quadrupole, octupole, hexadecapole and dotriacontapole field. The implemented multipole spectrum is taken from the MESSENGER results (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019).

In the present study, simulated magnetic field data and synthetically generated measurement positions are evaluated. Within the practical application to in situ spacecraft data, it is expectable that the measurement positions as well as the measurements will be determined defectively, resulting in estimation errors (e.g., Toepfer et al., 2021b). As a proof of concept, the errors are incorporated by synthetically disturbing the simulated data and the measurement positions. For example, disturbing the data by a normally distributed error of the width of 1 nT and zero mean results in a relative error between the disturbed $\stackrel{\mathrm{̃}}{\mathit{B}}$ and the undisturbed data B of $|\stackrel{\mathrm{̃}}{\mathit{B}}-\mathit{B}|/\left|\mathit{B}\right|\approx \mathrm{1}\mathit{\%}$. Disturbing the measurement positions by a normally distributed error of the width of 10 km and zero mean results in a defective shape matrix $\stackrel{\mathrm{̃}}{\mathbf{H}}$ (Toepfer et al., 2021b). The relative error between the defective and the ideal shape matrix H is given by $\left|\right|\stackrel{\mathrm{̃}}{\mathbf{H}}-\mathbf{H}|{|}_{\mathrm{2}}/|\left|\mathbf{H}\right|{|}_{\mathrm{2}}\approx \mathrm{3}\mathit{\%}$, where $\left|\right|.|{|}_{\mathrm{2}}$ denotes the spectral norm. The measurement errors and measurement position errors transfer onto disturbed estimators ${\stackrel{\mathrm{̃}}{\mathit{g}}}_{\text{TSVD}}$, ${\stackrel{\mathrm{̃}}{\mathit{g}}}_{\text{C}}$ and ${\stackrel{\mathrm{̃}}{\mathit{g}}}_{\text{T}}$ (Toepfer et al., 2021b). The unsigned difference between the coefficients resulting from the undisturbed data and measurement positions and the coefficients resulting from the disturbed data and measurement positions is used as an estimation for the error of the coefficients listed in Table 1. For all the coefficients together, the deviations between the undisturbed and the disturbed estimators result in

$$\begin{array}{}\text{(45)}& {\displaystyle \frac{|{\stackrel{\mathrm{̃}}{\mathit{g}}}_{\text{TSVD}}-{\mathit{g}}_{\text{TSVD}}|}{\left|{\mathit{g}}_{\text{TSVD}}\right|}}\approx {\displaystyle \frac{|{\stackrel{\mathrm{̃}}{\mathit{g}}}_{\text{C}}-{\mathit{g}}_{\text{C}}|}{\left|{\mathit{g}}_{\text{C}}\right|}}\approx {\displaystyle \frac{|{\stackrel{\mathrm{̃}}{\mathit{g}}}_{\text{T}}-{\mathit{g}}_{\text{T}}|}{\left|{\mathit{g}}_{\text{T}}\right|}}\approx \mathrm{1}\mathit{\%},\end{array}$$

so that the inversion methods may be declared as robust. However, the influence of specific measurement errors such as offsets, gains resulting from thermal variations and spacecraft magnetic disturbances (Narita et al., 2021) on the reconstruction procedure should be analyzed in future works.

Table 3 Implemented and reconstructed internal Gauss coefficients for the dipole, quadrupole, octupole, hexadecapole and dotriacontapole field. For the implemented quadrupole and octupole coefficients, the lower boundary from the MESSENGER results (Wardinski et al., 2019) is implemented.

The analysis of MESSENGER magnetic field data provided a value of $g_{\mathrm{5}}^{\mathrm{0}}\approx \mathrm{2}\mathrm{nT}$ (Thébault et al., 2018). Thus, the inversion methods are furthermore applied to simulated magnetic field data with the internal Gauss coefficients $g_{\mathrm{1}}^{\mathrm{0}}=-\mathrm{190}\mathrm{nT}$ (dipole field), $g_{\mathrm{2}}^{\mathrm{0}}=-\mathrm{78}\mathrm{nT}$ (quadrupole field), $g_{\mathrm{3}}^{\mathrm{0}}=-\mathrm{20}\mathrm{nT}$ (octupole field), $g_{\mathrm{4}}^{\mathrm{0}}=-\mathrm{6}\mathrm{nT}$ (hexadecapole field) and $g_{\mathrm{5}}^{\mathrm{0}}=\mathrm{2}\mathrm{nT}$ (dotriacontapole field). The data are again evaluated along the planned MPO trajectories. The following analyses will be restricted to undisturbed simulated data.

The optimal regularization parameter results in α≈0.77, so that κ(H_T)≈93. Thus, for the calculation of Capon's estimator and for the TSVD estimator, again 60 singular values are considered. The optimal diagonal loading parameter results in σ=670 nT. The estimators of the TSVD method, Capon's method and the Tikhonov regularization are displayed in Table 2.

Table 4 Implemented and reconstructed internal Gauss coefficients for the dipole, quadrupole, octupole, hexadecapole and dotriacontapole field. The implemented axisymmetric coefficients are taken from the MESSENGER results (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019). For the implemented nonaxisymmetric coefficients, arbitrary values are chosen.

The relative errors result in $|{\mathit{g}}_{\text{TSVD}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{3.3}$ %, $|{\mathit{g}}_{\text{C}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{2.2}$ % and $|{\mathit{g}}_{\text{T}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{2.3}$ %, so that the results are in agreement with the coefficients presented in Table 1. The performance of Capon's estimator again is as competitive as the performance of the Tikhonov estimator. Due to the smaller value of the internal coefficient $g_{\mathrm{5}}^{\mathrm{0}}$ implemented in the simulation, the relative deviation between the reconstructed and the implemented internal dotriacontapole results in about 50 %. Thus, if the true value of Mercury's internal dotriacontapole is of the order of 2 nT, uncertainties within the reconstruction procedure are expectable.

Besides the uncertainty concerning the internal coefficient $g_{\mathrm{5}}^{\mathrm{0}}$, the analysis of the MESSENGER data provided a range for the internal Gauss coefficients of the dipole up to the octupole field (Wardinski et al., 2019). Thus, it is worthwhile to analyze the performance of the estimators by evaluating magnetic field data resulting from the lower boundaries of the internal coefficients along the planned MPO trajectories. The inversion methods are applied to (undisturbed) simulated magnetic field data with the internal Gauss coefficients $g_{\mathrm{1}}^{\mathrm{0}}=-\mathrm{190}\mathrm{nT}$ (dipole field), $g_{\mathrm{2}}^{\mathrm{0}}=-\mathrm{57}\mathrm{nT}$ (quadrupole field), $g_{\mathrm{3}}^{\mathrm{0}}=-\mathrm{16}\mathrm{nT}$ (octupole field), $g_{\mathrm{4}}^{\mathrm{0}}=-\mathrm{4}\mathrm{nT}$ (hexadecapole field) and $g_{\mathrm{5}}^{\mathrm{0}}=\mathrm{8}\mathrm{nT}$ (dotriacontapole field). The implemented value of the hexadecapole field is chosen arbitrarily and the dotriacontapole field is again taken from the Mercury magnetic field model of Thébault et al. (2018). The optimal regularization parameter results in α≈1.57, so that κ(H_T)≈64, and thus, for the calculation of Capon's estimator and for the TSVD estimator, 58 singular values are considered. The optimal diagonal loading parameter results in σ=515 nT. The estimators of the TSVD method, Capon's method and the Tikhonov regularization are displayed in Table 3, where the relative estimation errors result in $|{\mathit{g}}_{\text{TSVD}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{6.6}\mathit{\%}$, $|{\mathit{g}}_{\text{C}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{3.6}\mathit{\%}$ and $|{\mathit{g}}_{\text{T}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{3.4}\mathit{\%}$. In analogy to the examples discussed above, the performance of Capon's estimator again equals the performance of the Tikhonov estimator. Due to the smaller numerical values of the implemented coefficients, the errors are slightly larger than the deviation between the reconstructed and the ideal coefficients presented in Tables 1 and 2, respectively.

Although Mercury's internal magnetic field is dominated by an axisymmetric geometry (Anderson et al., 2012; Wardinski et al., 2019), recent studies indicate the existence of nonaxisymmetric field contributions (Wardinski et al., 2021). Thus, the inversion techniques should be capable of reconstructing nonaxisymmetric internal fields to guarantee an unbiased analysis of the BepiColombo data. Therefore, the inversion methods are furthermore applied to simulated data resulting from arbitrary chosen nonaxisymmetric internal Gauss coefficients with m≠0. The optimal regularization parameter results in α≈1.39, so that κ(H_T)≈68. Thus, for the calculation of Capon's estimator and for the TSVD estimator, again 58 singular values are considered. The optimal diagonal loading parameter results in σ=560 nT. The estimators of the TSVD method, Capon's method and the Tikhonov regularization as well as the implemented values are displayed in Table 4. The relative estimation errors result in $|{\mathit{g}}_{\text{TSVD}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{6.0}$ %, $|{\mathit{g}}_{\text{C}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{3.7}$ % and $|{\mathit{g}}_{\text{T}}-\mathit{g}|/\left|\mathit{g}\right|\approx \mathrm{3.5}$ %. For an objective assessment of the reconstruction quality, the estimation errors of the axisymmetric and the nonaxisymmetric internal field contributions should be discussed separately. The relative deviation of the axisymmetric coefficients with m=0 results in 5.4 % for the TSVD estimator, 2.8 % for Capon's estimator and 2.0 % for the Tikhonov estimator, whereas the relative deviation of the nonaxisymmetric coefficients with m≠0 results in 34 % for the TSVD estimator, 32 % for Capon's estimator and 38 % for the Tikhonov estimator, which is much larger due to the implemented smaller numerical values. Thus, it is worthwhile to analyze the performance of the estimators for different combinations and numerical values of nonaxisymmetric internal Gauss coefficients in future studies. However, in case mother nature surprises us with nonaxisymmetric internal field contributions at Mercury as suggested by Wardinski et al. (2021), the inversion methods presented here enable the determination of the corresponding internal Gauss coefficients from the MPO data.

As newly established, Capon's method provides the same performance as the Tikhonov regularization. For the comparison of the two methods, a more general comment is appropriate. Both the methods incorporate the constraint of a minimum norm solution and, therefore, deliver superior results than the TSVD method. Within the derivation of the Tikhonov estimator g_T, the constraint is included synthetically, whereas the nature of Capon's method is based on the minimization of the output power (Toepfer et al., 2020a, b), which corresponds to the norm of the estimator. Furthermore, Capon's method weights the measurements by the data covariance and the measurement positions. Thus, the weighting of Capon's method is adaptive, since the data determine the weighting by themselves, whereas the Tikhonov method weights all data equally. In view of the practical application of the inversion methods, both the methods provide nearly comparable results. However, Capon's method incorporates more information from the experiment (position and data) than the Tikhonov regularization.

5 Summary and outlook

The detailed characterization of Mercury's internal magnetic field is expected to play an important role in understanding the origin of the field. Due to the interference of the internal parts with the magnetospheric field contributions, on the one hand, each part of the field has to be parametrized properly. On the other hand, a robust inversion method for reconstructing the wanted model coefficients is required.

In preparation for the analysis of the magnetic field measurements provided by the magnetometer on board the MPO, the plasma interaction of Mercury with the solar wind is simulated numerically. The resulting magnetic field data are parametrized by a combination of the Gauss representation with the Mie representation, called the Gauss–Mie representation, and the corresponding expansion coefficients are reconstructed from the data using the truncated singular value decomposition and the Tikhonov regularization as well as Capon's method. The reconstructed internal Gauss coefficients of the dipole, quadrupole, octupole, hexadecapole and dotriacontapole fields are in very good agreement with the coefficients implemented into the simulation code, and thus, a high-precision determination of Mercury's internal magnetic field up to the dotriacontapole is expectable. The quality of the reconstructed internal coefficients depends on the magnitude of the values. For example, in the case of an internal dotriacontapole coefficient of $g_{\mathrm{5}}^{\mathrm{0}}=\mathrm{8}\mathrm{nT}$, the reconstructed and the implemented coefficient are in very good agreement. For $g_{\mathrm{5}}^{\mathrm{0}}=\mathrm{2}\mathrm{nT}$ the deviation is of the order of about 1 nT. However, due to the symmetric distribution of the planned MPO orbits around Mercury (Heyner et al., 2021), it is worthwhile to consider the reconstruction of higher degrees of the internal field in future studies. Furthermore, the methods presented here are capable of reconstructing potentially existing nonaxisymmetric magnetic field contributions at Mercury, which enables an objective analysis of the prospective MPO data.

The comparison of the inversion methods shows that Capon's method and the Tikhonov method provide a comparative performance. Since both the methods incorporate the constraint of a minimum norm solution, which is equivalent to minimum energy, Capon's estimator as well as the Tikhonov estimator deliver superior results than the truncated singular value decomposition. It should be noted that the constraint of a minimum norm solution is included synthetically within the derivation of the Tikhonov estimator. Since Capon's method is based on the minimization of the output power, which corresponds to the norm of the estimator, the constraint of a minimum norm solution is naturally implemented in the method. Furthermore, Capon's method weights the data adaptively, since the weighting is determined by the measurement positions and the measurements themselves, whereas the Tikhonov method weights all data equally. Besides the constraint of a minimum norm solution, further physically based constraints, for example, at the core–mantle boundary (Wardinski et al., 2019), can be incorporated to improve the inversion results (Holme and Bloxham, 1996; Heyner et al., 2021).

Within the analyses presented here, simulated stationary magnetic field data resulting from the plasma interaction of Mercury with the solar wind under constant external conditions, i.e., constant solar wind density, velocity and interplanetary magnetic field orientation, are evaluated. Due to the fast temporal variability of the Hermean environment (e.g., Slavin et al., 2021) and changing solar wind conditions, a suitable filtering procedure should be applied before analyzing in situ data. Thereby, data from calm solar wind conditions are preferable. However, in view of the BepiColombo mission (Benkhoff et al., 2010, 2021), the present work provides an overview of the most commonly used inversion methods and shows that Capon's method as well as the Tikhonov method enable a high-precision determination of Mercury's internal magnetic field. The subsequent comparison of the reconstructed internal field resulting from the analysis of the BepiColombo data with the MESSENGER results will open the door to discuss the evidence of yet undetected secular variations at Mercury (e.g. Philpott et al., 2014; Oliveira et al., 2019; Heyner et al., 2021).