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 5 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.


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 10 is one of the primary goals of the BepiColombo mission (Benkhoff et al., 2010). Thereby, the detailed reconstruction of Mercury's internal field from the total measured magnetic field is of major importance for modelling 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. Due to the plasma physical interaction with the solar wind, the magnetic field around Mercury is contributed not only by the planetary internal field (such as the dynamo-generated field and the crustal remanent field) but 15 also by the external component resulting from currents flowing within the magnetosphere. 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 20 magnetic fields. For example, the Tikhonov regularization (Tikhonov et al., 1995), also known as L 2 -regularization has been 2 The Gauss-Mie representation The fluxgate magnetometer on board the Mercury Planetary Orbiter (MPO) (Glassmeier et al., 2010;Heyner et al., 2021) is 45 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)/2. Since the shell in general includes current-carrying regions with ∂ x × B = 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 i , B e and B sh . B i is the irrotational field contribution internally generated 50 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). Shortly, we call B i 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, , 1996Olsen, 1997;Toepfer et al., 2021a) the currents flowing within the shell in the region a < r < c generate toroidal B sh T and poloidal B sh P magnetic field structures that superpose with the existing internal and external fields, delivering which is called the Gauss-Mie representation of the magnetic field (Toepfer et al., 2021a).
Since both the internal and external fields B i and B e are purely poloidal and irrotational, these parts can be parametrized via the Gauss representation. Therefore, there exist scalar potentials Φ i and Φ e with and Using the body-fixed, planetary-centered spherical coordinates with radius r ∈ [R M , ∞), azimuth angle λ ∈ [0, 2π] and colatitude θ ∈ [0, π], the scalar potentials can be expanded into spherical harmonics, yielding and where R M is the planetary radius of Mercury and P m l are 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 70 epoch-making work by Gauss (1839), also revisited by Glassmeier and Tsurutani (2014) with the contemporary English translation. The potential Φ i is determined by the internal Gauss coefficients g m l and h m l , whereas the external Gauss coefficients q m l and s m l describe the external field contributions.
In general, the shell includes current-carrying regions and therefore, the contributions B sh T and B sh P cannot be described by 75 the Gauss representation. Making use of the Mie representation, one may introduce scalar functions Ψ sh T and Ψ sh P 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: and 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 Ψ sh T and Ψ sh P into spherical harmonics. In contrast to the scalar potentials Φ i 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 85 the spherical shell (Toepfer et al., 2021a), providing and where the variable ρ = (r − b)/R M denotes the signed radial distance from the mean shell radius b in units of the planetary 90 radius and the expansion coefficients are given by Due to the geometry of the MPO orbits the application of the thin shell approximation (Backus, 1986(Backus, , 1996Olsen, 1997;Toepfer et al., 2021a) is a valid assumption. 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 B sh P generated by toroidal currents flowing within the shell can be neglegted compared to the internal B i and external B e poloidal fields. Combining the Gauss representation with the Mie representation by making use 100 of the thin shell approximation, the magnetic field can be rewritten in the linear algebraic form 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.

105
It should be noted that the internal field is canonically described in a Mercury-Body-Fixed corotating coordinate system (MBF), whereas the external field is canonically described in a Mercury-Solar-Orbital system (MSO) with the x-axis orientated towards the sun, the z-axis orientated parallel to the rotation axis, i.e. antiparallel to the internal dipole moment, and the y-axis completes the right-handed system (Toepfer et al., 2021a;Heyner et al., 2021). In the present study, simulated Mercury 110 magnetic field data are evaluated. Therefore, the internal and external fields are expressed in a Mercury-Anit-Solar-Orbital coordinate system (MASO), where the x-axis is orientated towards the nightside of Mercury (away from the sun), the z-axis is orientated 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.

115
Thus, the resulting shape matrix H is rectangular and its inverse H −1 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.
For the reconstruction of Mercury's internal magnetic field, various kinds of data inversion techniques are available such as the 120 least square 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.

Least square fit
The most prominent inversion method for linear inversion problems is the least square fit method (LSF) (e.g., Haykin, 2014).
The method minimizes the quadratic deviation 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 square fit estimator where the dagger † denotes the Hermitian conjugate. Thus, the LSF method solely weights the data by the shape matrix H.

130
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).

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 H ∈ R m×n with rank n into the form where U ∈ R m×m and V ∈ R n×n are orthonormal transformations, so that U † = U −1 as well as V † = V −1 is valid (e.g., Connerney (1981); Haykin (2014); Connerney et al. (2018)). The matrix where Σ 1 ≥ · · · ≥ Σ n ≥ 0 (sorted in descending order), contains the so-called singular values Σ i of the shape matrix H, which 140 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 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 is called the "pseudo-inverse" or generalized inverse of the matrix H, where the matrix Σ + is given by 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 square fit estimator 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 of the shape matrix increases κ (H) 1 for small singular values and thus, the inversion problem (Eq. 10) is said to be 155 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 160 k < n, are considered within the solution and therefore, the shape matrix H is approximated by a shape matrix 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 is achieved. Due to the modification of the shape matrix, the solution g S transfers onto which will be called the TSVD estimator.

170
On one hand, the decreased condition number improves the solvability of the inversion method. 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)

175
Because of g T SV D = R g, where g denotes the ideal (imagenary) 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 exist 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).

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.

185
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) 190 Therefore, minimizing the energy spectrum W l is equivalent to the minimization of the norm of the wanted model coefficients g 2 and thus it is useful to search for solutions having the minimal norm. This constraint is known as the Tikhonov regularization, L 2 -regularization or minimum norm solution (e.g., Tikhonov et al. (1995); Haykin (2014)). Using this additional constraint, the original least square fit problem (cf. Eq. 11) can be extended as 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 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.

200
Inserting the singuar value decomposition of the matrix H into Eq. (26), delivers where the matrix will be called the Tikhonov inverse and

215
Thus, the original shape matrix H transfers onto the shape matrix where which has a lower condition number Therefore, the modification of the shape matrix improves the solvability of the inversion problem in analogy to the TSVD.
Comparison of the SVD estimator 230 with the Tikhonov estimator 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 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 non-trivial resolution matrix and thus, in general there exist model covariances.

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 245 robust inversion method for the analysis of Mercury's internal magnetic field (Toepfer et al., 2020a(Toepfer et al., , 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 (cf. Eq. 10), non-parametrized parts v as well as measurement noise n, so that is valid. The measurement noise is assumed to be Gaussian with variance σ 2 n 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 with respect to the distortionless constraint (also referred to as the unit gain constraint) where tr w † M w is the trace of the matrix w † M w and I is the identity matrix. The matrix M = B • B denotes the data covariance matrix. Capon's estimator realizing the minimal output power results in The comparison of Capon's estimator with the LSF estimator (cf. Eq. 12) shows that Capon's method can be interpreted as a special case of the weighted least square fit method (Toepfer et al., 2020b). The robustness of the method can be improved by the diagonal loading technique M → M + σ 2 d I, where σ d is the diagonal loading parameter (Toepfer et al., 2020b). Inserting the diagonally loaded data covariance matrix M = B • B + σ 2 I, where σ 2 = σ 2 n + σ 2 d , into the output power and making 265 use of g C = w † B , delivers Thus, Capon's method minimizes the energy spectrum in analogy to the Tikhonov regularization. In contrast to the Tikhonov 270 regularization, Capon's method additionally weights the data by the inverse data covariance matrix (cf. 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 275 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 280 is valid. Due to the application of the TSVD, Capon's filter matrix is modified to on the cost of violating the unit gain condition (Toepfer et al., 2021a), i.e.
4 Application to simulated Mercury magnetic field data 285 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 290 techniques for the virtual spacecraft data. The estimators for the spectrum are compared to the internal Gauss coefficients implemented in the simulation.

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., 295 Müller et al. (2011); Exner et al. (2018Exner et al. ( , 2020). The internal Gauss coefficients g 0 1 = −190 nT (dipole field), g 0 2 = −78 nT (quadrupole field), g 0 3 = −20 nT (octupole field), g 0 4 = −6 nT (hexadecapole field) and g 0 5 = 8 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 0 5 = 8 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 orientated along the vector (x, y, z) T = (0, 0, 1) T in the Mercury-Anti-  (Milillo et al., 2020). The solar wind proton density number is chosen to n sw = 30 cm −3 and the electron T e and proton temperatures T p are implemented to T e = T p = 2.5 × 10 5 K (Exner et al., 2020;Milillo et al., 2020). The resulting magnitude of the magnetic field in the x-z-plane is displayed in Figure 1.
The geometry of Mercury's magnetosphere is mainly dominated by the internal dipole field. Also the internal quadrupole 305 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.

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, magnetic field B sh T 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 Ψ sh P of the poloidal magnetic field B sh P 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 a priori known as the inputs to the 320 simulation, in the following discussion we will focus on the reconstructed 25 internal coefficients.
The optimal regularization parameter for the application of the Tikhonov regularization results in α ≈ 0.9, so that κ (H T ) ≈ 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 330 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 0 5 of the dotriacontapole field was chosen from a synthetic Mercury magnetic field model of Thébault et al. (2018). The relative errors result in g T SV D − g / g ≈ 3.9%, g C − g / g ≈ 2.6% and g T − g / g ≈ 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 335 Tikhonov estimator, since both methods incorporate the constraint of a minimum norm solution. Especially, 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

Coefficient Input in nT TSVD in nT Capon in nT Tikhonov in nT
340 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 R C = R = I. In this case, the deviation between Capon's estimator and the ideal solution results in g C − g / g ≈ 3.3% for σ = 420 nT and g L − g / g ≈ 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 345 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).

355
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 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). 190.0 -195.5 -190.9 -191.5 The relative errors result in g T SV D −g / g ≈ 3.3%, g C −g / g ≈ 2.2% and g T −g / g ≈ 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 360 performance of the Tikhonov estimator. Due to the smaller value of the internal coefficient g 0 5 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.

365
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)). For example, disturbing the simulated data by a normally distributed error of the width of 1 nT and zero mean and simultanously disturbing the measurement positions by a normally distributed error of the width of 10 km and zero mean, the relative estimation errors result in 370 g T SV D − g / g ≈ 3.3%, g C − g / g ≈ 2.3% and g T − g / g ≈ 2.4%. Since these errors are of the same order as the errors resulting from the undisturbed data, the inversion methods may be declared as robust.
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 375 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 380 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.
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 one hand, each 385 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 390 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, 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 hexadecapole is expectable. The 395 quality of the reconstructed internal dotriacontapole field depends on the magnitude of the value. In the case of g 0 5 = 8 nT, the reconstructed and the implemented coefficients are in very good agreement. For g 0 5 = 2 nT the deviation is of the order of about 1 nT. However, due to the symmetric distribution of the planned MPO orbits around Mercury , it is worthwile to consider the reconstruction of higher degrees of the internal field in future studies.

400
The comparision of the inversion methods shows that Capon's method and the Tikhonov method provide a comparative performance. Since both 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 405 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).

410
In view of the BepiColombo mission (Benkhoff et al., 2010), 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. Toepfer, S., Narita, Y., Heyner, D., Motschmann, U.: Error propagation of Capon's minimum variance estimator,Front. Phys.,9:684410,