Geodesy, Geophysics and Fundamental Physics Investigations of the BepiColombo Mission

Abstract

In preparation for the ESA/JAXA BepiColombo mission to Mercury, thematic working groups had been established for coordinating the activities within the BepiColombo Science Working Team in specific fields. Here we describe the scientific goals of the Geodesy and Geophysics Working Group (GGWG) that aims at addressing fundamental questions regarding Mercury’s internal structure and evolution. This multidisciplinary investigation will also test the gravity laws by using the planet Mercury as a proof mass. The instruments on the Mercury Planetary Orbiter (MPO), which are devoted to accomplishing the GGWG science objectives, include the BepiColombo Laser Altimeter (BELA), the Mercury orbiter radio science experiment (MORE), and the MPO magnetometer (MPO-MAG). The onboard Italian spring accelerometer (ISA) will greatly aid the orbit reconstruction needed by the gravity investigation and laser altimetry. We report the current knowledge on the geophysics, geodesy, and evolution of Mercury after the successful NASA mission MESSENGER and set the prospects for the BepiColombo science investigations based on the latest findings on Mercury’s interior. The MPO spacecraft of the BepiColombo mission will provide extremely accurate measurements of Mercury’s topography, gravity, and magnetic field, extending and improving MESSENGER data coverage, in particular in the southern hemisphere. Furthermore, the dual-spacecraft configuration of the BepiColombo mission with the Mio spacecraft at higher altitudes than the MPO spacecraft will be fundamental for decoupling the internal and external contributions of Mercury’s magnetic field. Thanks to the synergy between the geophysical instrument suite and to the complementary instruments dedicated to the investigations on Mercury’s surface, composition, and environment, the BepiColombo mission is poised to advance our understanding of the interior and evolution of the innermost planet of the solar system.

1 Introduction

The internal structure of the planet Mercury preserves crucial information regarding the formation and evolution of the Solar System. A thorough characterization of its interior is one of the major goals of the ESA/JAXA mission BepiColombo that will start orbiting the innermost planet of the Solar System in 2025 (Benkhoff et al. 2010). Mercury was previously visited by two NASA missions only. Mariner 10 flew by Mercury three times in 1974–1975, providing a detailed imaging and temperature mapping of the planet’s surface, and unprecedented and surprising evidences of the presence of a weak internal magnetic field, and of a tenuous atmosphere (Murray 1975). The MErcury Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER) spacecraft was launched in August 2004 (Solomon et al. 2007), and after three Mercury flybys in 2008–2009 and ∼7 years of cruise, it was inserted in a highly eccentric polar orbit around the planet in March 2011. The orbital configuration was maintained with a pericenter altitude of ∼200–500 km in the northern hemisphere (∼70^∘N) and an apocenter altitude of ∼10,000–15,000 km for the first three years. In April 2014, a final extended mission (XM2) was approved to lower the pericenter altitude to 5-25 km for several weeks until Mercury impact on 28 April 2015. The science phase of the MESSENGER mission enabled outstanding findings on the surface composition, interior structure, and magnetic environment (Solomon et al. 2018). However, high-resolution results were limited to the planet’s northern hemisphere because of the high-eccentric spacecraft orbit with its pericenter in the northern hemisphere.

The BepiColombo mission consists of two spacecraft that will survey Mercury from two different orbits. The Mercury Magnetospheric Orbiter (MMO, or Mio) was built by JAXA to study the interaction between Mercury’s magnetosphere and the solar wind from a 590×11,640 km polar orbit (Kasaba et al. 2020). The Mercury Planetary Orbiter (MPO) was designed by ESA to operate in a less eccentric 480×1,500 km orbit that will yield a uniform global coverage of Mercury’s surface (Benkhoff et al. 2010). The MPO orbit configuration was conceived after the first MESSENGER results of Mercury’s gravity field that showed large values of the low–degree zonal harmonics (Smith et al. 2012; Genova et al. 2013). The precession of the MPO pericenter due to the gravitational perturbations will enable, for the first time, the exploration of Mercury’s southern hemisphere from altitudes lower than 500 km (Imperi et al. 2018).

The instruments onboard the MPO spacecraft will provide accurate measurements of the properties of Mercury’s interior, surface, exosphere, and magnetic field (Benkhoff et al. 2010). Geodesy and geophysical investigations will aim at addressing the open questions regarding Mercury’s internal structure by accurately observing the topography, magnetic, and gravity field of the planet. The BepiColombo laser altimeter (BELA) will measure Mercury’s topographic relief, rotational state (e.g., spin rate and amplitude of the physical librations) and tidal deformations (Thomas et al. 2007). The MPO magnetometer (MPO–MAG) will provide data of Mercury’s magnetic environment, focusing on the structure of the internal field (Glassmeier et al. 2010). The Mercury Orbiter Radio science Experiment (MORE) will allow retrieving Mercury’s gravity field, rotation, and ephemeris through the extremely precise orbit determination of the MPO spacecraft (Iess et al. 2009). An accurate knowledge of the planet’s orbit around the Sun will also yield tests of theories of gravitation including Einstein’s theory of General Relativity (Einstein 2019). MORE investigations will be supported by the onboard Italian Spring Accelerometer (ISA) that will acquire measurements of the non-gravitational forces (e.g., solar radiation pressure) (Iafolla et al. 2010). The synergy between BELA, MPO–MAG, MORE, and ISA will be fundamental for a comprehensive understanding of Mercury’s internal structure.

This paper focuses on the joint efforts of the BepiColombo science teams involved in the Geodesy and Geophysics Working Group (GGWG). In Sect. 2, we present the science objectives of the GGWG activities, including the fundamental physics experiment conducted by the MORE team (Iess et al. 2009). In Sect. 3, we present the instruments that will acquire geodetic and magnetic measurements for the geophysical investigations outlined in Sect. 4. Finally, we discuss the resulting geophysical constraints on Mercury’s internal structure and evolution in Sect. 5.

2 Scientific Objectives of the Multidisciplinary Investigations

The MESSENGER and BepiColombo missions were conceived to address key scientific questions regarding Mercury’s origin and evolution, and its surrounding environment (Solomon et al. 2007; Benkhoff et al. 2010). The measurements acquired by the MESSENGER spacecraft enabled the accomplishment of the NASA mission goals, but raised also fundamental questions (Solomon et al. 2018), which will be investigated by the BepiColombo mission. Table 1 shows the BepiColombo GGWG science themes, questions, and objectives. The MPO will provide extremely highly accurate measurements with the onboard science instruments, and a uniform planet coverage through the MPO lower eccentric orbit than MESSENGER. A global view of Mercury and the improved MPO data qualities will help answer the fundamental questions that are still open.

Table 1 Science themes, questions, and objectives of the BepiColombo Geodesy and Geophysics Working Group investigations

The origin and evolution of Mercury stand out from the major themes of the BepiColombo multidisciplinary investigations. A better understanding of the planetary formation will be achieved by determining properties of Mercury’s interior and surface. Mariner 10 (Murray 1975), Earth-based radar (Margot et al. 2007), and MESSENGER observations (Smith et al. 2012) demonstrated the presence of a large core that is consistent with Mercury’s high ratio of metal to silicate. Different scenarios have been presented to describe the processes that led to the measured metal/silicate ratio (Weidenschilling 1978; Cameron 1985; Fegley Jr. and Cameron 1987). These cases, however, result in different predictions of the properties of Mercury’s outer silicate layer. A significant refinement of surface chemistry and mineralogy is one of the main BepiColombo science goals. The gravity and topography investigations will yield Mercury’s global crustal thickness and bulk density. The Mercury Imaging X-ray Spectrometer (MIXS) and Mercury Gamma-ray and Neutron Spectrometer (MGNS) will determine the elemental composition of Mercury’s crust (Fraser et al. 2010; Mitrofanov et al. 2010). Multispectral imaging and spectroscopic data collected by the Spectrometer and Imagers for MPO BepiColombo Integrated Observatory System (SIMBIO-SYS) will reveal crustal differentiation, space weathering, and rock minerals abundances (Flamini et al. 2010). The synergetic analysis of these complementary datasets will enhance our knowledge of the physical parameters of the outer layers, including particle size, strength, and porosity.

High-resolution gravity and topographic global maps will provide crucial information on Mercury’s geological history. The surface of the planet hosts important records of the past endogenic and exogenic activities. The comparison between gravity and topography is fundamental for estimating the level of internal compensation, which informs the evolution of the planetary crust. Mercury’s surface, which will be mapped by the MPO SIMBIO-SYS images (Flamini et al. 2010), consists of heavily cratered regions, intercrater plains, hilly and lineated terrain (e.g., antipodal to the Caloris basin), and smooth plains (Spudis and Guest 1988). A global view of these areas will allow us to determine the time frame when surface regions formed. By observing wrinkle ridges and lobate scarps (e.g., Enterprise Rupes above the Rembrandt basin), the timing and the amount of secular internal cooling will be constrained by using BepiColombo data. Thermal evolution models suggest a planetary radius contraction of ∼4-10 km (Solomon 1977; Van Hoolst and Jacobs 2003; Dombard and Hauck 2008), which is in contrast with the results based on the imaging data from Mariner 10 (Watters et al. 1998) and MESSENGER flybys (Di Achille et al. 2012). A more recent analysis of images collected by MESSENGER during its orbital mission yielded a value of global contraction closer to the expectations (∼7 km) (Byrne et al. 2014). BepiColombo high-resolution images combined with topographic altimetry profiles will enable a refined mapping of the tectonic features in the southern hemisphere to determine the radius contraction. Extensive features will also be studied through gravity gradiometry, which is based on the computation of the second spatial derivatives of the gravitational potential (Andrews-Hanna et al. 2013). The distribution and orientation of the gravity gradient anomalies indicate the stress state induced by the radius contraction.

Magnetic field measurements may provide further constraints on Mercury’s evolution. MESSENGER data showed that the crust is strongly magnetized (Johnson et al. 2015), which suggests that Mercury once may have had a stronger magnetic field than Earth. These measurements were obtained during the XM2 low-altitude campaign of the mission when the spacecraft was closer than $\sim130$ km to the surface. Measuring crustal magnetization will be more challenging for the MPO spacecraft because of its high initial pericenter. This will decline over the course of the mission—how fast depends on the gravity field which has to be further explored as described here. Also, it remains to be seen how long MPO will survive in the Hermean harsh thermal environment. However, MPO–MAG will attempt to constrain first crustal magnetic anomalies in the southern hemisphere at least at large spatial scales.

A detailed characterization of Mercury’s internal magnetic field through the MPO–MAG measurements will allow us to constrain the internal dynamo process that operates in the planet’s core. A better understanding of Mercury’s deep interior will provide critical insights into the thermal evolution of the planet. To determine accurately the internal structure of Mercury, gravity and altimetry investigations will measure the dimensionless polar moments of inertia of the whole planet $\left (\frac{C}{MR^{2}}\right )$ and the fractional polar moment of inertia of the silicate outer layers $\left (\frac{C_{cr+m}}{MR^{2}}\right )$ (Peale et al. 2002). These geophysical quantities will be retrieved by estimating two Mercury’s rotational parameters: the pole obliquity, $\epsilon$, and the amplitude of the physical longitudinal librations, $\phi _{0}$. Furthermore, gravity and altimetry data will allow us to measure Mercury’s gravitational (i.e., Love number $k_{2}$) and radial surface (i.e., Love number $h_{2}$) tidal responses, respectively. An accurate estimation of both Love numbers $k_{2}$ and $h_{2}$ will provide a strong constraint on the size of the outer liquid core (Padovan et al. 2014) and of the solid inner core (Steinbrügge et al. 2018). The retrieval of the gravitational phase lag, which informs on the level of internal dissipation, will also enhance our knowledge of the rheology (i.e., viscosity and rigidity) of the mantle, which plays a major role in Mercury’s thermochemical evolution (Tosi et al. 2013). Therefore, the GGWG joint analysis of altimetry, gravity, and magnetic results will be crucial to precisely recover the thermal state and configuration of Mercury’s internal structure.

Mercury’s orbit around the Sun and its 3:2 spin–orbit resonance help understand how the planet formed and evolved (Correia and Laskar 2004; Wieczorek et al. 2012). Refinements of Mercury’s ephemeris will indicate possible departures from this minimum-energy condition. The MORE radio science instrument will enable a precise determination of Mercury’s orbit by using extremely accurate range measurements. A thorough investigation of Mercury’s ephemeris also represents a great opportunity to conduct fundamental physics and heliophysics experiments. Because of the planet’s proximity to the Sun, Einstein’s theory of General Relativity (GR) (Einstein 2019) must be accounted for to accurately model its orbital evolution (i.e., perihelion precession). The sub-meter precision of the MORE range data will be well-suited to estimate the parameterized post-Newtonian (PPN) parameters $\beta$, $\gamma$, $\alpha _{1}$, and $\alpha _{2}$ to test possible GR violations. Einstein’s theory of GR is based on the assumption that the ratio between gravitational ($m_{g}$) and inertial ($m_{i}$) masses is equal to 1. MORE range data will enable the detection of potential discrepancies in this equality condition due to the self-gravitational energy of the planets in the solar system (i.e., Strong Equivalence Principle, SEP) (Nordtvedt 1968). An enhanced knowledge of Mercury’s orbital motion will also inform on the interior structure of the Sun by measuring the solar $GM$, $J_{2}$, and $\frac{\dot{GM}}{GM}$, which depends on the time variation of the gravitational constant, $G$, and the solar mass loss rate due to solar radiance and wind.

3 BepiColombo Science Instrumentation

The MPO spacecraft hosts four main instruments dedicated to the geodetic and geophysical investigations of the BepiColombo mission. BELA, MORE, ISA, and MAG are briefly described in Sects. 3.1, 3.2, 3.3, and 3.4, respectively. Further details of these instruments are reported in the papers by Thomas et al., Iess et al., Santoli et al., and Heyner et al. of this issue.

3.1 BEpiColombo Laser Altimeter (BELA)

The MPO laser altimeter, BELA, will acquire range measurements to precisely determine the relative distance between the spacecraft and Mercury’s surface from altitudes <1055 km. The estimate of this upper bound altitude is based on instrument specifications providing global surface coverage. It will be refined in orbit around Mercury since it depends on several unknown parameters (e.g., short-scale surface roughness of Mercury). BELA will provide a global surface coverage, including first altimetric observations of the southern hemisphere. The range data is computed by accurately measuring the time–of–flight of a short (∼5 ns) laser pulse emitted from the instrument to the planetary surface and back-scattered to the receiver. The start time is recorded by transferring a small fraction of the emitted pulse directly to the receiver through fiber optics. From a precise two-way time-of-flight measurement between the transmitted and received wave-package, which is ∼3.2 ms and ∼7 ms for altitudes of 480 km and 1055 km, respectively, topographic data can be retrieved when combined with spacecraft position and attitude data. Topographic models of the planetary surface at global, regional, and local scales will be derived from several hundred million laser shots collected during the BepiColombo nominal mission.

The onboard software of BELA is capable of analyzing the return pulse by using polynomial fits to approximate the pulse shapes. If requested, the fully digitized pulse can be returned to Earth. Furthermore, the shape of the return pulse provides information on the surface albedo at the laser wavelength and on the roughness of the surface on the scale of the laser footprint (∼16 to 53 m diameter, depending on spacecraft altitudes).

BELA is equipped with two redundant Nd:YAG-lasers, capable of generating 50 mJ laser pulses at 1064 nm wavelength. The lasers can be operated from 1 to 10 Hz. The receiver is a Cassegrain-type telescope with an aperture of 20 cm and a field of view of 247.5 μrad (half cone). For detection of the reflected laser pulses an APD (Avalanche Photo Diode) is used. Details of the instrument design can be found in Thomas et al. (this issue) and in Thomas et al. (2019).

The signal is transmitted via the Analog Electronics Unit to the Range Finder Module (RFM), where the laser pulses are processed and transmitted to the Digital Processing Module (DPM), the control board of the instrument and the interface to the spacecraft platform. The transmitted and returned pulses are sampled with a bin-size of 12.5 ns which would correspond to a range resolution of 1.875 m. However, due to filter-matching algorithms within the range finder electronics, a sub-sampling accuracy smaller than 1.5 ns corresponding to a range resolution of better than 20 cm can be achieved under optimum conditions. The range error is also affected by surface slope and roughness. Extremely steep terrains (i.e., ∼40^∘ slopes) cause range errors that are still below 80 cm even for measurements collected at the estimated detection threshold (i.e., 1050–km spacecraft altitude) (Steinbrugge et al. 2018).

The precision of the BELA measurements relies on the calibration of the following error sources: (a) small misalignments of the transmitter with respect to the spacecraft reference frame, which induce instrument pointing errors; (b) electronics and clock drifts; and (c) orbital errors that affect the georeferentiation of the altimeter data. The pointing errors will be significantly reduced by adopting measurement techniques that enable a refined calibration of the transmitter pointing with respect to the SIMBIO–SYS imaging system (Stark et al. 2017). Since both transmitted and received signals, undergo the same electronic chain, the electronic delays cancel out. Furthermore, the range finder clock drifts on long time scales (much longer than laser pulse time–of–flight) will be calibrated by using a precise onboard pulse per second signal (PPS) over the mission duration. The MPO orbital errors will be significantly mitigated by the MORE team, which will provide precise trajectory reconstructions by processing the X/X/Ka Doppler data. By assuming uncorrelated errors, which also account for possible mismodeling of Mercury’s rotational state, the overall range measurement error is less than 10 m with assumed roughness values of 12.1 m at 200-m baseline and 6.4 m at 50-m baseline, and a mean albedo of 0.19 (Steinbrugge et al. 2018). By analyzing the energy of the transmitted and returned pulses BELA will also be sensitive to Mercury’s surface albedo at the laser wavelength of 1064 nm.

After the successful launch of the BepiColombo mission, three BELA instrument check-outs have been performed (Near-Earth Commissioning Phase and the first two cruise checkouts). Since BELA is facing the Mercury Transfer Module (MTM), firing the laser is not possible due to the enhanced risk of harmful back-reflections into the instrument. Therefore, functional check-outs focused on the receiver chain and general housekeeping data. Dark noise levels of the detector were recorded and the first data for calibrating the long–term drift of the RFM clock with respect to the on-board clock were obtained.

3.2 Mercury Orbiter Radio Science Experiment (MORE)

The MPO radio science experiment, MORE, will enable the precise orbit determination of the spacecraft to accurately estimate the physical quantities responsible for the dynamical evolution of both the MPO and Mercury’s trajectories. The radio science data are acquired by Earth’s ground stations during dedicated radio tracking passages. The spacecraft receives a signal from one of the ESA’s Deep Space Antennas (DSA) and sends it back to the same DSA (i.e., two–way link) or another ground station (i.e., three–way link) to establish telecommunication operations. The radio tracking measurements that are acquired by the DSA station are range and range-rate observables. The time delay and the Doppler shift of the received signal measure the spacecraft relative distance and velocity in the line–of–sight between the spacecraft and the ground station, respectively. The radio tracking data are then processed in orbit determination software to reconstruct the spacecraft trajectory and to accurately adjust the parameters of scientific interest that affect spacecraft and/or central body orbital motion.

The standard configuration of a radio science instrument is usually based on the Telemetry, Tracking and Command (TT&C) subsystem, which includes a transponder (i.e., Deep Space Transponder, DST) for a single X-band uplink (∼7.2GHz) and a two coherent downlinks in X- (∼8.4 GHz) and Ka-band (∼32 GHz), respectively (Asmar et al. 2005). The MORE instrument includes a 2-m-diameter steerable high-gain antenna and a dedicated transponder, the Ka-Transponder (KaT) (De Tiberis et al. 2011; Ciarcia et al. 2013), that enables Ka-band (∼34GHz–∼32 GHz) up- and down-link functionalities. Since the plasma is a dispersive medium (i.e., waves of different frequencies travel at different velocities), the MORE multi–frequency X/X/Ka configuration allows calibrating for charged particle effects, which significantly reduces the level of noise of standard radio tracking systems (Fig. 1). The requirements of range-rate and range accuracies were set to 3 μm s⁻¹ at 1000-s integration time (i.e., 12 μm s⁻¹ at 60 s) and 20 cm, respectively. Tests of the MORE radio system carried out in May 2019 showed better than expected accuracies. Range data attained a sub-cm accuracy with 4-s integration time when the spacecraft was at 0.3 AU (Cappuccio et al. 2020).

Fig. 1

Block diagram of BepiColombo radio tracking scheme

The radio science instrument of the MESSENGER mission was designed for a single X-band radio link (Srinivasan et al. 2007). This configuration was significantly affected by fluctuations of the solar plasma in proximity of superior solar conjunctions (Iess et al. 2014). MESSENGER radio data acquired at Sun–probe–Earth (SPE) angles larger than 30^∘ showed accuracies of ∼1-2 m and ∼50 μm s⁻¹ at 60-s integration time for range and range-rate, respectively (Genova et al. 2018, 2019). This level of noise was mainly due to thermal effects induced by the telecommunication system. At lower SPE angles, the plasma noise caused larger errors in both range and range-rate data.

The high quality of the MORE radio tracking data will guarantee significantly enhanced accuracies of orbit and gravity determination compared to the MESSENGER spacecraft. To compensate undesired effects due to mismodeling of perturbing forces, the MPO hosts the ISA instrument (Iafolla et al. 2010; Santoli et al. 2018), which is described in Sect. 3.3. The joint analysis of radio and accelerometer data will yield an extremely accurate estimation of the geophysical parameters and an excellent knowledge of the MPO orbit, which will be used in BELA and SIMBIO–SYS data processing.

3.3 Italian Spring Accelerometer (ISA)

The ISA instrument is a high-sensitivity three-axis accelerometer devoted to providing highly accurate measurements of the MPO non-gravitational perturbations (Iafolla et al. 2010, 2011). The dynamical evolution of the spacecraft orbit will be strongly affected by Mercury’s gravity field, and by non-conservative forces, including solar, planetary albedo and thermal infrared radiation pressures (Lucchesi and Iafolla 2006). An accurate knowledge of these non–gravitational accelerations is fundamental to achieve precise orbit and gravity determination through the processing of MORE radio tracking data.

The assembly of the ISA instrument consists of three mono–axial accelerometers, arranged to form an orthogonal reference frame within the MPO spacecraft (Fiorenza et al. 2016). Each sensor is manufactured from a single piece of aluminium Al7075, carved by a milling machine, to obtain a proof mass (i.e., the sensing element) suspended on an external frame through a tiny foil-shaped spring (∼140 μm in thickness) as shown in Fig. 2. This spring-mass system is a mechanical oscillator with a natural frequency of ∼3.6 Hz. ISA is designed to measure acceleration signals in the frequency range of 3×10⁻⁵ – 10⁻¹ Hz and with a maximum amplitude of 3×10⁻⁶ m s⁻², enabling accuracies up to 10⁻⁸ m s⁻² (Iafolla et al. 2010; Fiorenza et al. 2016).

Fig. 2

(left) Configuration of ISA sensors, and (right) naked accelerometer and foil-shaped spring (Images courtesy of Thales Alenia Space)

Two pairs of symmetric plates face the central proof mass and realize four capacitors. A pair of capacitors, named pick-up plates, is used to measure the displacement of the proof–mass from its equilibrium position when the frame undergoes an acceleration (i.e., capacitive transducer). The capacitive bridge is biased with a modulated signal of $f_{p} = 10$ kHz and it is decoupled by an isolation transformer. Any acceleration at frequency $f_{s}$ ≪ $f_{p}$ induces a movement of the proof mass, and hence a modulation of the bias voltage: at the output of the bridge the signal is seen at the two side bands $f_{\pm }=f_{p} \pm f_{s}$. The signal is digitized by an ADC (Analog to Digital Converter) and then demodulated. A second pair of capacitors, named actuators, is used to apply electrostatic forces to the sensing mass (Fig. 2). Actuators have three functions: to recenter the mass at its working position, to provide a calibration signal used in-flight to calibrate the transduction factor of pick-up chain, and finally to damp the resonance of the mechanical oscillator to reach an amplification factor at resonance frequency.

3.4 Magnetometer (MPO-MAG)

The instruments of the BepiColombo mission include a dual-sensor magnetometer onboard each spacecraft. Mio hosts the MMO/MGF to study Mercury’s magnetosphere, and interplanetary solar wind. A similar instrument, MPO-MAG, is onboard the MPO spacecraft to provide measurements more relevant for geophysical applications. The MPO-MAG consists of two fluxgate magnetometers, which measure the three magnetic field components. The two sensors per spacecraft are required to distinguish between natural signals and magnetic disturbances internally originated by the probe subsystems (e.g., reaction wheels and other instruments) as discussed by e.g., Ness et al. (1971). The working principle of the magnetometers is described in the work by Glassmeier et al. (2010). The primary science goal of the MPO-MAG team is to improve our knowledge of the global internal field of Mercury (Heyner et al. 2021). Due to the low amplitude of the planetary magnetic field, the expected signals are very weak. Therefore, an enhanced understanding of the spacecraft generated disturbances and the magnetosphere, which will be investigated by the MMO/MGF instrument onboard the Mio spacecraft (Baumjohann et al. 2010), is required.

The two MPO magnetometers are mounted on a boom, which was deployed in the Near-Earth Commissioning Phase. Since then, the MPO-MAG instrument has been fully operational and has collected data continuously in order to characterize the spacecraft magnetically.

Table 2 gives an overview of the experiment characteristics and the spacecraft effects on the instrument as verified in space so far. During cruise, the offsets can be corrected on a routine basis, using the approach established by Hedgecock (1975).

Table 2 MPO-MAG instrument characteristics. The instrument noise stated here is for a temperature of $T=180\;^{\circ }$C, which has not been reached so far. The values for static and dynamic fields related to internal spacecraft sources are derived from the limited dataset collected so far. These measurements enable a first assessment only, since few instruments are already operative. The magnetic disturbance levels, furthermore, exclude solar electric propulsion where disturbances up to 160 nT are expected. The sensor alignment has also been checked during the Earth flyby of the BepiColombo composite by using the well-known terrestrial magnetic field

Given the predicted trajectory of the MPO spacecraft and the instrument performance in Table 2, we can estimate the accuracies of the internal planetary magnetic field inversion. The resulting uncertainties mainly depend on the MPO orbit configuration, offset determination (i.e., instrument readings in zero ambient field), and orientation knowledge. A conservative inversion estimate yields a solution of the internal magnetic field in spherical harmonics to degree and order 6 as it is described in detail by Heyner et al. in this issue. The predicted error in the Gauss coefficients becomes as large as the field coefficients beyond degree 6. More robust inversion techniques and data cleaning algorithms could enhance the model resolution.

4 Geophysical Models and Measurements

4.1 Shape and Topography

Topographic data are essential for understanding local and regional processes that have shaped the planetary surface. The BELA instrument will significantly improve our knowledge of Mercury’s shape and topography at different scales by retrieving a global network of laser tracks. Benefitting from global coverage, BELA will also provide refined estimates of Mercury’s rotational parameters. Mercury’s spin rate, obliquity, as well as amplitude and phase of the physical librations will be adjusted by using techniques that were developed for the MESSENGER mission (Stark et al. 2015). Time-varying perturbations will also be estimated, including a key objective as the tidal Love number $h_{2}$ that informs on the deformation of the planet due to tides raised by the strong gravity field of the Sun (Steinbrugge et al. 2018; Thor et al. 2020).

4.1.1 Mercury’s Topographic Map

The analysis and processing of the BELA measurements will enable an accurate mapping of Mercury’s topographic relief. The global coverage and horizontal resolution of Mercury’s topography will be mainly constrained by the MPO trajectory and the BELA performances, which are modeled through the probability of false detection (PFD). The PFD relies on the instrument characteristics, the spacecraft altitude, and the physical properties of the surface within the laser footprint, including roughness and albedo. By assuming BELA flight model tests and Mercury’s surface characteristics, our current performance models indicate that the PFD is close to zero (i.e., 100% successful measurements) when the MPO will be at altitudes lower than 1400 km, 1000 km, and 700 km over terrains with slopes of $0^{\circ }$, $20^{\circ }$, and $40^{\circ }$, respectively (Steinbrugge et al. 2018). These results suggest that BELA will provide a uniform global coverage at mission completion. To estimate the horizontal resolution of the topographic map, we considered a PFD of $< 20$% and BELA operations for the entire nominal and one-year extended mission. Figure 3 shows that the horizontal resolution varies from ∼3 km at the equator down to less than ∼250 m at latitudes above 80^∘ and below −80^∘. This spatial resolution corresponds to a maximum degree 1100 in spherical harmonics. However, the topographic elevation model after the nominal mission will show a lower resolution because of the longitudinal gaps in BELA surface coverage.

Fig. 3

Horizontal resolution of the topographic map based on BELA profiles obtained after two years of operation in Mercury orbit (Steinbrugge et al. 2018)

The large amount of BELA profiles will allow us to accurately map high- and mid-latitude geological features, including lobate scarps and wrinkle ridges (Watters et al. 2015). A vertical resolution of <1 m in optimum conditions will be fundamental for an accurate determination of the height of geological features and units (e.g., central peak, hollows, and rim of impact craters). A detailed mapping of these features also depends on the gaps between individual laser spots. The diameter of BELA footprint will be 24 and 53 m at the MPO altitudes of 480 and 1055 km, respectively. By assuming 10 Hz shot frequency and the nominal MPO orbit configuration, the gaps between contiguous laser spots will vary between 170 and 250 m, enabling uniform along-track coverage for the orbital ground-tracks. Due to the orbit, BELA will obtain the densest coverage in the north and south polar regions. However, geological features in the equatorial and mid-latitude regions (also south), which have been poorly explored before will be mapped accurately.

By sampling and analyzing the digitized return pulse, the BELA instrument will have the capability to determining the pulse broadening that is indicative of slope and roughness at the footprint scale of ∼50 m. If the effect of the slope is subtracted from a sequence of laser spots, the pulse-spreading is a measure of the surface roughness on the footprint spatial scale. Correlation (or anti-correlation) with geological units will provide crucial information on the processes that have shaped the surface.

Complementary measurements of Mercury’s topography over the entire surface will be obtained through stereo photogrammetric analysis of the SIMBIO-SYS data (Flamini et al. 2010). The Stereo Imaging Channel (STC) of SIMBIO-SYS will provide a ground sampling resolution of 40-150 m accordingly to the MPO pericenter latitude (Slemer et al. 2018). The combination of BELA and SIMBIO-SYS measurements will yield extremely accurate high-resolution maps of Mercury’s digital elevation models.

Mercury’s shape will be determined by BELA with a lateral resolution of <2.5 km (Fig. 3). An accurate knowledge of the shape’s orientation and the offset between the center-of-mass and the center-of-figure of the planet will inform on Mercury’s non-hydrostatic state. Analyses of the MESSENGER data provided estimations of this offset of 140 m (Perry et al. 2015) and $185\pm 45$ m (Stark et al. 2017). An equatorial rotation of the degree-2 shape relative to the principal axes of $\sim 17^{\circ }$ was also measured, suggesting asymmetries in the deep compensation (Perry et al. 2015). Mercury’s gravity field retrieved by the processing of MESSENGER radio science data, furthermore, indicates significant deviations from the hydrostatic state for the current tidal forces and rotational state (Smith et al. 2012). MORE and BELA investigations will provide global coverage and high resolution of both gravity and altimetry data, leading to a better characterization of the center-of-mass and the center-of-figure offset, and the relative orientation of the shape with respect to the principal axes.

4.1.2 Radial Tidal Deformation

Tidal forces exerted by the Sun cause time-varying deformations of Mercury’s surface. Because of the 3:2 spin–orbit resonance, the main tidal cycle is the 88-days orbital period. The tidal effects are parameterized by the second-degree Love numbers $h_{2}$, $l_{2}$, and $k_{2}$. These parameters describe the surface and gravitational response of the planet to external tidal forces, and depend on the properties of the planet’s interior, including its internal density structure and the rheology (e.g., rigidity and viscosity). The Love number $k_{2}$ describes the change in the gravitational potential due to the re-distribution of mass in the planet’s interior, and will be estimated by MORE (Sect. 4.2.2). The Love numbers $h_{2}$ and $l_{2}$ measure the radial and lateral surface displacement, respectively. Whereas the estimation of the Love-Shida number $l_{2}$ is only possible with a landed element, constraints on the Love number $h_{2}$ can be obtained by analyzing BELA data as a function of time. Figure 4 shows the maximum radial deformation (peak-to-peak) as a function of longitude and latitude for each point on Mercury’s surface. The tidal deformation is proportional to the Love number $h_{2}$, which is assumed here accordingly to Mercury’s interior structure and rheological modeling. These surface deformations are 60 cm and 200 cm at the poles and equatorial regions, respectively.

Fig. 4

Maximum tidal amplitudes (peak-to-peak) on Mercury’s surface over one tidal cycle of 89 days, i.e., one revolution around the Sun. Here a typical Love number $h_{2} = 0.85$ is assumed (Hussmann and Stark 2020)

To accurately determine the Love number $h_{2}$, BELA will probably have to collect data over two years. Because of orbital and operational constraints, an extended one-year mission will enable spatial and temporal distributions of the BELA data that are well-suited to detect the surface tides (Steinbrügge et al. 2018; Thor et al. 2020). Altimetric measurements collected over intersecting ground-tracks at different tidal phases (i.e., crossovers) will be processed to precisely determine the tidal radial deformation. The MPO orbit configuration will provide ∼60M crossovers after two years of operations. By processing these measurements, the Love number $h_{2}$ will be determined with an accuracy of ∼0.14, which is $\sim 18\%$ of the previously assumed a priori $h_{2} = 0.8$ (Steinbrügge et al. 2018). A complementary approach to determine the Love number $h_{2}$ is based on a global inversion of the altimetric data by using cubic B-splines to model the local static topography for a measurement point (Koch et al. 2010; Thor et al. 2020). Thor et al. (2020) show that the processing of the BELA measurements through this novel technique will yield an $h_{2}$ accuracy of 0.012, which may be subject to uncertainties in periodic misalignment behavior of the instrument.

4.1.3 Outer Layer Orientation

A precise modeling of Mercury’s crust and mantle rotation allows constraining the physical state, density, and size of the outer core. Earth-based observations of Mercury’s surface by Margot et al. (2007) revealed that the rotation of the outer layers is decoupled from the outer core. An amplitude of the longitudinal libration, $\phi$, of ∼400 m was observed. Libration amplitude measurements based on MESSENGER imaging and altimetry data confirmed this result with refined estimates (Stark et al. 2015). Comparable accuracies of the libration amplitude (i.e., ∼1-2 arcsec) are expected from the analysis of the BELA data (Koch et al. 2008; Rosat et al. 2008), and images from SIMBIO-SYS high-resolution imaging channel (HRIC) (Pfyffer et al. 2011; Aboudan et al. 2014). The combination of MESSENGER altimetric and imaging data, and the BELA measurements will allow us to estimate the amplitude of the long-period librations that are induced by the other planets in the solar system (Yseboodt et al. 2013). A first measurement of the long-period librations may significantly constrain the properties of the solid inner core. The obliquity of Mercury’s pole also provides information on the deep interior, and it will be estimated through the analysis of the BELA measurements.

As a consequence of Mercury’s small obliquity, the solar incidence angle at Mercury’s polar regions is close to zero. In topographic lows (e.g., crater floors) near the poles, the terrain can remain in permanent shadow. Due to the lack of an atmosphere and highly insulating regolith, H₂O ice and other volatiles brought to the inner solar system by comets can be stable over millions of years in these permanently shadowed regions. Both polar regions will be prime targets for altimetry and reflectivity measurements (Neumann et al. 2013; Chabot et al. 2014a) of the BELA instrument at the 1064 nm wavelength.

4.2 Gravity

An accurate knowledge of a planetary gravity field enables an in-depth characterization of the planet’s interior. The internal mass distribution induces gravitational anomalies that cause dynamical perturbations on the trajectory of spacecraft in orbit around the planet. By processing the MPO radio science data, the MORE team will provide an extremely precise estimation of Mercury’s gravity field from local to global scales. Mercury’s gravitational anomalies preserve information on the structure and properties of its core, mantle, and crust. The deep interior will be also investigated by adjusting gravitational tides, and rotational parameters.

4.2.1 Gravitational Field

The modeling of Mercury’s gravity field is based on the following spherical harmonic expansion (Kaula 2000)

$$ U=\frac{GM}{r}\left \{ 1+\sum _{l=2}^{l_{\mathrm{max}}}\left ({\frac{R}{r}} \right )^{l}\left [\sum _{m=0}^{l}\left (\overline{C}_{lm}\cos {m\phi }+ \overline{S}_{lm}\sin {m\phi }\right )\overline{P}_{lm}(\cos \theta ) \right ]\right \} , $$

(1)

where $GM$ and $R$ are the gravitational constant and radius of the planet, respectively; $l$ and $m$ are the degree and order of the normalized spherical harmonic coefficients $\overline{C}_{lm}$ and $\overline{S}_{lm}$; $\overline{P}_{lm}$ are the normalized associated Legendre functions; and $\phi$, $\theta$, and $r$ are longitude, colatitude, and relative distance from the center of the planet, respectively. The normalization factor adopted in this representation is $\sqrt{(2-\delta _{m0})(2l+1)\frac{(l-m)!}{(l+m)!}}$ with $\delta _{m0}$, the Kronecker delta, equal to 0 and 1 for $m\ne 0$ and $m = 0$, respectively. The coefficients $\overline{C}_{lm}$ and $\overline{S}_{lm}$ are adjusted in the gravity solution to determine the gravity anomalies associated with Mercury’s internal mass distribution. The maximum degree, $l_{\mathrm{max}}$, of the estimated spherical harmonic coefficients constrains the spatial resolution of the gravity field. Higher degrees provide crucial information on the finer spatial scales of the gravity field. A better understanding of the gravity anomalies at different spatial scales enables a more comprehensive geophysical investigation of Mercury’s internal structure from the outer silicate layers to the core.

The radio science team of the MESSENGER mission provided very accurate models of Mercury’s gravity field, named HgM008, and orientation (Genova et al. 2019). These geophysical results were retrieved by processing the entire MESSENGER radio science dataset, which includes the measurements collected during the low-altitude campaign of the extended mission. The HgM008 gravity model reports the global field with the associated uncertainties in spherical harmonics to degree and order 100 (Fig. 5). MESSENGER gravity mapping provided an uneven coverage of the planetary surface because of the highly eccentric orbit of the spacecraft. Therefore, the resolution in spherical harmonics to degree 100, which corresponds to a spatial resolution of ∼80 km, is only achieved in the northern hemisphere. To determine the local resolution of the HgM008 gravity field, a degree strength map was determined by using the gravity model and its covariance matrix (Konopliv et al. 1999). The retrieved gravity field enables the computation of the following expected radial accelerations

$$ a_{l}=\frac{GM}{R^{2}}\left (l+1\right )\sum _{m=0}^{l}\left ( \overline{C}_{lm}\cos {m\phi }+\overline{S}_{lm}\sin {m\phi }\right ) \overline{P}_{lm}(\cos \theta ). $$

(2)

To yield a profile of these accelerations that relies on the gravity degree only, a Kaula power rule is introduced as follows

$$ C_{l}=\sqrt{\frac{1}{2l+1}\sum _{m=0}^{l}(\overline{C}_{lm}^{2}+ \overline{S}_{lm}^{2})}=\frac{A_{k}\times 10^{-5}}{l^{2}}. $$

(3)

The coefficient $A_{k}$ is directly determined by the measured gravity field. The power spectrum of the HgM008 gravity model is consistent with a coefficient $A_{k} = 4$ of the Kaula rule. The RMS magnitude spectrum of the predicted radial accelerations is given by:

$$ \left (a_{l}\right )_{RMS} = \frac{GM}{R^{2}}\sqrt{\frac{2}{n}}\left (A_{k} \times 10^{-5}\right ). $$

(4)

This value is then compared to the acceleration uncertainty, which is stored in the gravity covariance matrix. Profiles of the acceleration uncertainties, $\sigma (a_{2,l})$, are computed by accounting for the covariance matrix from degree 2 to $l$ ($\boldsymbol{P}_{2,l}$). The acceleration uncertainty is:

$$ \sigma (a_{2,l}) = \frac{\partial a_{2,l}}{\partial \boldsymbol{G}_{2,l}}^{T} \boldsymbol{P}_{2,l} \frac{\partial a_{2,l}}{\partial \boldsymbol{G}_{2,l}}. $$

(5)

The vector $\boldsymbol{G}_{2,l}$ includes all the normalized gravity coefficients from degree 2 to $l$. To compute the acceleration uncertainty at degree $l$, we subtract the contribution of all the degrees up to $l-1$ $\left (i.e., \sigma (a_{l}) = \sigma (a_{2,l}) - \sigma (a_{2,l-1}) \right )$. The intersection between the predicted acceleration, $\left (a_{l}\right )_{RMS}$, and the retrieved uncertainty, $\sigma (a_{l})$, provided the maximum spherical harmonic degree for each latitude and longitude (Konopliv et al. 1999). The degree strength map of HgM008 shows accuracies close to degree 90 in regions of the northern hemisphere, where MESSENGER had a pericenter altitude of ∼5-20 km. The equatorial region and the southern hemisphere show poorer resolutions close to degree $l=15$ (i.e., spatial resolution of ∼500 km).

Fig. 5

Root mean squared power spectra of MESSENGER HgM008 gravity solution. Thin colored lines show the formal uncertainties of HgM008 (blue) and MORE solution after one–year of BepiColombo nominal mission (red)

The MORE radio science investigation will enable an accurate mapping of Mercury’s gravity field in these regions. Extremely precise radio tracking data and the spacecraft lower altitudes in the southern hemisphere compared to MESSENGER will provide unprecedented measurements of gravitational accelerations associated with internal mass anomalies. By comparing the measured gravity signal (i.e., free-air gravity) with the expected gravity from topography (i.e., Bouguer correction), we will improve our understanding of the processes that led to the formation and evolution of surface features (e.g., Rembrandt crater) revealing important information on Mercury’s geological history. Figure 6 shows the degree strength map of Mercury’s gravity field solution retrieved through the processing of the MORE radio science data simulated over the entire 1–year nominal mission, and the extra–year for the extended mission. The resulting gravity map of the MORE investigation will enable spatial resolutions in the southern mid-latitudes of ∼170–190 km ($l=40$–45). This refined gravity field of Mercury will allow revealing the properties of geological units in the southern hemisphere.

Fig. 6

Degree strength map for the MORE gravity solution after one-year nominal mission (top) and a second year of the extended mission (bottom) over shaded global Digital Elevation Model (Becker et al. 2016) in a Mollweide projection

4.2.2 Gravitational Tides

A detailed characterization of the properties of Mercury’s mantle will be achieved by estimating the Love number $k_{2}$ and its phase lag. The retrieval of these parameters through the analysis of the MESSENGER radio science data was significantly limited by orbital perturbations due to non-conservative forces (e.g., solar radiation pressure). Therefore, MESSENGER gravity solutions only provided an accurate recovery of the amplitude of planet’s gravitational tidal response $k_{2}=0.569\pm$0.025 (Genova et al. 2019). This value unambigously confirms the presence of a liquid inner core but also allows some conclusions on the mantle rheology. The bulk composition of Mercury’s mantle is assumed to be close to enstatite chondrite or bencubbinite chondrite meteorites (Malavergne et al. 2010). Based on this assumption, different rheological models have been conceived (Padovan et al. 2014) suggesting unrelaxed rigidities between 59 and 71 GPa. However, two significant unknowns are the iron content and the grain size. For the former, current rheological models assume no iron content which is justified by the low surface abundance of FeO (Nittler et al. 2011). While the effect of small amounts of iron-rich minerals on the rheological properties of the mantle would be small for most minerals, olivine shows a strong dependence on the iron content (Zhao et al. 2009) and could possibly affect the assumed mantle rigidity. The other unknown, the grain size, directly affects the viscosity of the mantle. Lower viscosities lead to higher Love numbers. General ranges are assumed to range from mm to cm scale. Within this range, the other significant factor influencing the viscosity is the mantle temperature. The most recent measurement of $k_{2}$ supports the presence of a warm and weak mantle (Padovan et al. 2014).

Our results of the numerical simulations of the MORE experiment show a dramatic enhancement in the estimation of the Love number $k_{2}$ (Table 3). The full compensation of the non-conservative forces through the processing of the ISA measurements lead to improved determination of Mercury’s gravitational tides, including its phase lag (Table 3). The joint solutions of $k_{2}$ amplitude and phase lag will lead to strongly constrain the rheological properties of Mercury’s mantle.

Table 3 Formal uncertainties (3 standard deviations, 3–$\sigma$) of a set of geophysical parameters estimated in the HgM008 model (Genova et al. 2019), and by the MORE investigation after BepiColombo nominal and extended mission

4.2.3 Deep Interior Orientation

A better coverage of the equatorial regions will provide highly accurate measurements of Mercury’s rotation and orientation. Table 3 shows the resulting formal uncertainties of the rotational parameters estimated through the numerical simulations of the MORE experiment. These results, which were obtained by simulating BepiColombo nominal and extended mission, are compared to the HgM008 gravity model. The analysis of the MORE radio science data will enable refined estimates of the pole coordinates (i.e., right ascension, $\alpha$, and declination, $\delta$), which constrain the planet’s obliquity ($\epsilon$). MESSENGER low–altitude campaign provided radio science data that were well-suited to enhance our understanding of Mercury’s orientation. The obliquity of the HgM008 solution is $\epsilon = 1.968\pm$0.027 arcmin (Genova et al. 2019) that is fully consistent with the Cassini state, which represents the main assumption to determine the dimensionless polar moment of inertia ($C/MR^{2}$) as a function of $\epsilon$ (Peale et al. 2002).

The MORE gravity investigation will provide an accurate estimate of Mercury’s obliquity with 3–$\sigma _{\epsilon }$ formal uncertainties of 0.007 and 0.004 arcmin after the nominal and extended mission, respectively. This level of accuracy of the planet’s obliquity will improve our knowledge of $C/MR^{2}$ (Sect. 5.2.1). A full characterization of the properties of Mercury’s internal structure requires a precise determination of the fractional polar moment of inertia of the solid crust and mantle ($C_{cr+m}/C$). This geophysical quantity depends on the amplitude of the physical longitudinal librations, $\phi _{0}$, which will be measured by MORE with outstanding 3–$\sigma$ formal uncertainties of 2-5 m, which are 0.5–1% of the total amplitude (i.e., ∼400 m; Margot et al. 2007). Highly accurate solutions of both $C/MR^{2}$ and $C_{cr+m}/C$ will result in sophisticated modeling of Mercury’s deep interior including size and status of its core (Sect. 5.3).

4.3 Planetary Magnetic Field

Planetary magnetic fields consist of internal as well as external parts. Internal fields are the dynamo field, the field of the magnetized crust, induced fields from external field variations, and, as a unique characteristic of Mercury, fields from subsurface currents that close vertical currents from the magnetosphere. In this section, we discuss various approaches to map and interpret the different contributions.

4.3.1 Measurements and Models of Mercury’s Magnetic Field

The magnetic field of Mercury was investigated by the NASA missions Mariner 10 (Ness et al. 1975) and MESSENGER (Anderson et al. 2012). Two of the three flybys of the Mariner 10 spacecraft were close enough to detect a global planetary magnetic field in interaction with the solar wind. A much more complete survey of its magnetic environment was accomplished by the MESSENGER spacecraft during the 4 years of the nominal and extended mission.

To separate the dynamo generated internal field from the other contributions, some assumptions are required. In the absence of local currents, the magnetic field $\underline{B}$ may be described as the gradient of the scalar potential $\varPsi$:

$$ \underline{B} = - \nabla \varPsi . $$

(6)

The representation of the magnetic field as gradient of a scalar potential $\varPsi$ leads to:

$$ \nabla \times \underline{B} = \mu _{0} \mathbf {j} = \nabla \times (- \nabla \varPsi ) \equiv 0 . $$

(7)

Thus, local currents within the magnetosphere cannot be described with the scalar potential initial guess. Because the magnetic field is solenoidal, this potential obeys the Laplace equation:

$$ \Delta \varPsi = \nabla \cdot (\nabla \varPsi ) = \nabla \cdot \underline{B} = 0 . $$

(8)

The magnetic potential may be conveniently represented in the spherical harmonic expansion (i.e., eigenfunctions of the Laplace operator) of the internal ($\varPsi _{\text{int}}$) and external ($\varPsi _{ \text{ext}}$) magnetic potentials:

$$\begin{aligned} \varPsi _{\text{int}} =& R \sum _{\ell =1}^{L}\sum _{m=0}^{ \ell } \left (\frac{R}{r}\right )^{\ell +1} (g_{\ell }^{m} \cos (m \varphi ) + h_{\ell }^{m} \sin (m\varphi ) ) P_{\ell }^{m} (\cos (\theta )) \end{aligned}$$

(9)

$$\begin{aligned} \varPsi _{\text{ext}} =& R \sum _{\ell =1}^{L}\sum _{m=0}^{ \ell } \left (\frac{r}{R}\right )^{\ell } (G_{\ell }^{m} \cos (m\varphi ) + H_{\ell }^{m} \sin (m\varphi ) ) P_{\ell }^{m} (\cos (\theta )) . \end{aligned}$$

(10)

These potentials are defined in a body-fixed, planetocentric, spherical coordinate system with $r$, $\theta$, and $\varphi$ denoting distance to the planet’s center, colatitude and azimuth, respectively. The parameters $R$, $\ell$, $m$, and $L$ stand for the mean planetary radius, and degree, order, and cutoff degree of the spherical harmonics, respectively. The associated Legendre polynomials $P_{\ell }^{m}$ differ from the functions $\overline{P}_{lm}$, which are defined in the spherical harmonic expansion of the gravity field (Sect. 4.2.1), because of a different normalization. Magnetic investigations adopt the Schmidt-normalized coefficient $\sqrt{(2-\delta _{m0})\frac{(l-m)!}{(l+m)!}}$ with $\delta _{m0}$, the Kronecker delta, defined in Sect. 4.2.1. The lower case $g_{\ell }^{m}$ and $h_{\ell }^{m}$ denote the Gauss coefficients of the internal field, while the upper case $G_{\ell }^{m}$ and $H_{\ell }^{m}$ denote the external field contributions. The internal contributions of Mercury’s magnetic environment stem from the dynamo, induced fields, closure currents of field-aligned currents and crustal magnetic fields. Its external fields are generated by the planet’s interaction with the solar wind and are of significant magnitude at the planetary surface. Depending on the location on the surface, the external fields can reach strengths from roughly $1/10$ to unity relative to the internal fields.

The Mauersberger-Lowes spectrum is defined by:

$$ W_{\ell }= (\ell +1) \left ( \frac{R}{r}\right )^{(2\ell +4)} \sum _{m} \left [(g_{\ell }^{m})^{2} + (h_{\ell }^{m})^{2} \right ] , $$

(11)

and represents the energy content for each spherical harmonic degree $\ell$. Summing over all $\ell$ for a given $m$ yields the respective expression for each spherical harmonic order. A unique expansion of the field in spherical harmonic functions requires dense and evenly distributed data. Any departure from a non-uniform distribution of the magnetic data results in high correlations among the Gauss coefficients, which can be defined by using several independent representations. The orbit of the MESSENGER spacecraft allowed collecting planetary magnetic field measurements in the northern hemisphere only. Different field models, therefore, fulfill equally well the observations, which are not well-suited to fully disentangle the equatorially symmetric and anti-symmetric contributions.

The Gauss coefficient correlation due to orbital restrictions may be quantified as follows. The spherical harmonic representation may be written as a linear vector transform:

$$ \underline{\underline{A}} \; \underline{m} = \underline{b} . $$

(12)

Here, $\underline{\underline{A}}$, $\underline{m}$ and $\underline{b}$ stand for the design matrix, which depends on the measurement locations ony, the Gauss coefficient vector, and the magnetic field measurements vector, respectively. The correlation matrix $\underline{\underline{R}}$ is derived from the model parameter covariance matrix $Cov(\underline{m}) =: \underline{\underline{\varSigma }}$. Following Menke (Menke 2018), this matrix is defined as:

$$ \underline{\underline{\varSigma }} = \sigma _{d}^{2} [ \underline{\underline{A}}^{T} \underline{\underline{A}}]^{-1} . $$

(13)

The data a priori variance, $\sigma _{d}$, may be derived, for example, from the measurement uncertainty. Then, Fahrmeir et al. (2003) define the correlation matrix:

$$ \underline{\underline{R}} = \left (\mathit{diag}\left ( \underline{\underline{\varSigma }}\right )\right )^{-1/2} \underline{\underline{\varSigma }} \left (\mathit{diag}\left ( \underline{\underline{\varSigma }}\right )\right )^{-1/2} . $$

(14)

According to this definition, the model parameter correlation matrix depends on the measurement locations only.

Despite a large number of measurements of the magnetic field, strong correlations of the Gauss coefficients are expected because of the orbital geometry of the MESSENGER spacecraft. Johnson et al. (2018) and Thébault et al. (2017) analyzed the model parameter correlation. By assuming the standard spherical harmonic representation of internal and external fields, high correlations (or anti-correlations) are obtained between the following coefficient pairs (up to degree 4):

• $g_{1}^{0} \leftrightarrow g_{2}^{0}$

• $g_{1}^{0} \leftrightarrow g_{3}^{0}$

• $g_{1}^{0} \leftrightarrow g_{4}^{0}$

• $g_{1}^{0} \leftrightarrow G_{1}^{0}$

• $g_{2}^{0} \leftrightarrow G_{1}^{0}$

• $g_{1}^{1} \leftrightarrow G_{2}^{1}$.

Other correlations not stated here may be strong, but at the same time irrelevant as the respective coefficients are negligible in Mercury’s magnetic field models obtained after the MESSENGER mission. Figure 7 (top) shows the correlation matrix of the internal Gauss coefficients to degree and order 6 retrieved by the analysis of the MESSENGER data.

Fig. 7

Correlation matrices of the internal Gauss coefficients retrieved after the MESSENGER mission (top) and through the numerical simulations of the BepiColombo mission (bottom). The order of the coefficients is: $g_{1}^{0}, g_{1}^{1}, h_{1}^{1}, g_{2}^{0}, g_{2}^{1}, \dots$. Terms with $m\neq 0$ are not shown for readability reasons

Several approaches have been proposed to determine the interior magnetic field of Mercury. We refer to a static field model derived from the entire MESSENGER magnetic field dataset as an average main field model. The different modeling approaches are based on independent schemes of data reduction and filtering. Therefore, a direct comparison of the average misfits is not straightforward. Since the MESSENGER magnetic field data were partially acquired on regions that are non current-free and with strong time-variable external fields from the magnetosphere, different techniques have been used to obtain the inversion.

From the analysis of the data taken during Mercury flybys of the Mariner 10 spacecraft, Ness et al. (1974) already concluded that a weak axial dipole dominates the internal magnetic field, but its center is shifted northwards. Ness (1978) discuss several other analysis attempts by assuming different internal field geometries and structures of the external field. The internal axial dipole coefficient was estimated to range from 170 to 350 nT. Strong axial quadrupole contributions were also identified, demonstrating the non-uniqueness of the inversion caused by a limited spatial measurement coverage before the arrival of the MESSENGER mission. After the first two Mercury flybys of the MESSENGER spacecraft, Alexeev et al. (2010) employed a relatively simple magnetospheric model to prevent high correlations among the coefficients of internal and external fields. The combination of the Mariner 10 measurements with the MESSENGER two-flyby dataset yielded an axial dipole moment of −196 nT $R^{3}$ and a northward dipole offset of $0.165 R$. The analysis of early orbital MESSENGER data by Anderson et al. (2011) confirmed the dipole moment with a slightly corrected offset of $0.198 R$.

A precise definition of the offset dipole is important to describe the structure of the internal field. The offset between magnetic and geographic equator was assessed by Anderson et al. (2011) $\rho _{z} = \sqrt{x^{2}+y^{2}}$ > 1.29 $R$. Extrapolation of this offset dipole model towards Mercury’s core is not advisable as harmonic contributions decay with the distance from the planet, and the decay rate increases with the spherical harmonic degree (see equation Eq. (10)). Thus, in general, the magnetic field complexity is expected to increase towards the planet’s interior. When the field is axisymmetric, for very large distances to the rotation axis $\left (\rho _{z} \gg 0\right )$, the magnetic equator shows an offset in the $z$-direction by a constant value $d$:

$$ d\approx \frac{1}{2} R \frac{g_{2}^{0}}{g_{1}^{0}} . $$

(15)

At these large distances, this field then appears to be the field of a pure axial dipole with an offset equal to $d$. The internal field in the dynamo region, however, can be completely different.

Alexeev et al. (2010) show that Eq. (15) can be generalized to include all axial harmonics:

$$ g_{l}^{0} = \gamma _{1}^{0} l \left (\frac{d}{R}\right )^{l-1} . $$

(16)

Here, $\gamma _{1}^{0}$ denotes the Gauss coefficient of the pure offset dipole. Equation (16) describes a field where the offset $d$ is constant at any distance to the planet. In theory, this equation allows us to calculate any axial harmonic for a given axial dipole contribution, once $d$ is known. In practice, the determination of the magnetic equator is not precise enough to constrain the tiny contributions of higher harmonics.

The analysis by Anderson et al. (2011) provided each single magnetic equator crossing and the average $d$. By analyzing all available magnetic data within Mercury’s magnetosphere, deviations from the specific offset dipole series in Eq. (16) were also computed. These results provide evidence of unique properties of Mercury’s dynamo field that have been confirmed by other independent studies. The magnetic field of Mercury is relatively weak and very axisymmetric, with a dipole tilt < $1^{\circ }$. The resulting axisymmetric Gauss coefficients are (in ascending degrees): −190, −74.6, −22.0, 5.7 nT. This is close to the specific offset dipole coefficient series but shows noticeable differences. Non-axisymmetric Gauss coefficients are all below 3 nT. Dipole and quadrupole contributions dominate, and the quadrupole to dipole ratio is relatively large.

Alternative internal field models were developed with respect to the dipole offset (or the quadrupole-to-dipole ratio $g_{2}^{0} / g_{1}^{0}$ as in Eq. (15)). These model approaches, however, provide a certain centered multipole combination only. The Gauss coefficient correlations due to the MESSENGER orbit were accounted for in different ways.

Johnson et al. (2018) adopted a semi-empirical magnetospheric model introduced by Korth et al. (2014) to describe the external field. The resulting internal field, which was modeled as an offset axial dipole, presented an offset of $0.195 R$ and a slightly weaker dipole moment of −188 nT $R^{3}$. The average misfit between this model and observations is only 9.6 nT during magnetically quiet orbits (i.e., minor magnetospheric activity as defined by Anderson et al. (2013)).

Local field models of the northern hemisphere were also introduced to fit MESSENGER magnetic measurements. Oliveira et al. (2015) adapted a crustal field modeling technique (i.e., equivalent source dipoles) to model the core field of Mercury during the first six months of the MESSENGER mission. Thébault et al. (2017) considered another localized approach (i.e., revised spherical cap analysis) to model the field in the northern hemisphere by processing the entire MESSENGER mission dataset. Both analyses used MESSENGER measurements acquired at spacecraft altitudes $<1000$ km. The Gauss coefficients of the model by Oliveira et al. (2015) yield a $g_{2}^{0}/g_{1}^{0} = 0.212$ that is equivalent to $d= 0.106 R$. The local model by Thébault et al. (2017) was converted into a global field in spherical harmonics to degree and order $\ell = 5$, and the resulting quadrupole-to-dipole ratio is 0.27, which also results in a less important equivalent dipole offset of $d=0.135 R$. A similar result was retrieved by Wardinski et al. (2019) through an independent global spherical harmonic approach, which is based on a priori regularization constraint.

The different estimates of Mercury’s magnetic dipole offset suggest a significant model dependency. Furthermore, the pure offset dipole modeling is unable to fully describe the magnetic equator locations. To demonstrate this inconsistency between measurements and theoretical modeling, we linearly fitted the measured $B_{\rho }$ component (1 min averages) in the equatorial vicinity ($200 < Z_{\text{MSO}} < 800$ km) by following the method presented by Anderson et al. (2011). The parameter $Z_{\text{MSO}}$ and $\rho _{z}$ are defined in the Mercury-Solar-Orbital (MSO) planetocentric coordinate system with the $x$-axis pointing from the planet towards the Sun, the $y$-axis pointing against Mercury’s orbital motion tangentially to the orbital plane, and the $z$-axis completes the right-handed system pointing north. Figure 8 shows the equator crossings in the $Z_{\text{MSO}}\rho _{z}$-plane of the measured $B_{\rho }$ compared to the internal field model predictions. The magnetic equator $z$-locations increase with increasing $\rho _{z}$. High altitude crossings beyond $\rho _{z}=1.45R$ are much more disturbed than the low altitude ones, which may be due to the neutral sheet current reacting to upstream solar wind variations (Rong et al. 2018). A restriction of the input data to lower altitudes results in probably too low dipole offsets.

Fig. 8

Location of the Z-coordinate of the magnetic equator as function of distance from the rotational axis ($\rho$). Grey plus signs mark the locations of the magnetic equator in the magnetic field data. The red line shows the averaged equator locations in $\rho$-bins (“histogram”). The green line displays the equator location for a pure offset dipole ($Z_{ \text{offset}} = 0.196R$). The purple, orange and yellow lines depict the equator locations for the internal field models of Anderson et al. (2011) and Thébault et al. (2017)

Local currents may then cause undesired effects in the estimation of the Gauss coefficients. To account for the presence of local currents, the basic assumption that the region of analysis is current free has to be violated by using, for example, a field-aligned current (FAC) system (Anderson et al. 2018). At Earth, the FAC are grouped in three different regions. At Mercury, only the region-1 FAC analogue has been observed so far. The region-1 FAC flow from the dawn-side magnetopause along the magnetic field lines towards the planet and re-emerge on the dusk-side to flow again towards the magnetopause. These local currents are quasi-steady and may produce a magnetic signature pole-ward of $60^{\circ }$N of up to 60 nT. Sub-surface Pederson-like currents within the planet’s interior (primarily the lower mantle and the core-mantle boundary) may provide current closure. The magnetic fields from these closure currents with an amplitude $>=50$ nT are also observable from the spacecraft altitudes. A detailed configuration of the FAC system is currently under debate (Exner et al. 2020). Furthermore, Pederson-like closure currents within the planet’s interior produce magnetic fields that must appear as internal multipoles in the field analysis.

The BepiColombo measurements of Mercury’s magnetic field will help to adequately model the effects of the local current systems. The MPO–MAG instrument will provide crucial magnetic data with a more uniform planet’s coverage than MESSENGER. The measurements collected from the MPO altitudes will enable an accurate characterization of the properties of Mercury’s magnetic field, including a better understanding of the dipole offset. Precise measurements of the magnetic field in the southern hemisphere will be fundamental for breaking the high correlations among the Gauss coefficients obtained through the analysis of the MESSENGER data. Figure 7 (bottom) shows the correlation matrix of the magnetic field inversion retrieved with numerical simulations of the MPO–MAG investigation. The MPO orbit geometry will enable a dramatic decrease of the correlations among the internal parameters, except for $g_{1}^{0} \leftrightarrow g_{3}^{0}$ correlation.

4.3.2 Induced Magnetic Fields

The external fields from the magnetosphere are relatively strong compared to the interior magnetic field. As the magnetosphere is also quite dynamic, these external fields change on short timescales. Furthermore, the magnetosphere changes accordingly to the solar wind conditions along its 88-day orbit around the Sun. Most importantly, the solar wind dynamic pressure, $p_{sw}$, significantly changes with the heliocentric distance, $r_{hel}$. Because of the lack of a significant ionosphere of Mercury, the planet is exposed to these magnetic field changes, and thus, relatively strong induced fields from the planet’s interior are expected.

These induced magnetic fields depend on the internal conductivity structure, which in turn depends on the composition and thermophysical properties of the planet. The analysis of continuous measurements of these induced fields will be, therefore, a very important task of the magnetometer instrument team to provide crucial information on Mercury’s internal structure complementary to geodetic and gravitational investigations. To accomplish this science objective, two-point measurements are needed to independently recover the external (inducing) and internal (induced) fields. The BepiColombo mission will allow disentangling these two contributions by measuring the external field through the Mio spacecraft, and by monitoring the inductive response with the MPO spacecraft.

Induction within the planet’s interior plays an important role by influencing the possible compression of the magnetosphere into the planet (Hood and Schubert 1979). Different approaches have been proposed to account for the induced fields at Mercury. A first estimate of the induced field was retrieved by Grosser et al. (2004) with a simplified two-shell planetary conductivity model. By assuming a core size of 1840 km and a dipole field of $g_{1}^{0} =300$ nT, the induced fields yield 7–12% of the mean surface intensity of the internal planetary field..

The magnetic measurements of the MESSENGER mission enables the modeling of the induced fields without assuming the properties of Mercury’s core. Johnson et al. (2016) modeled the annual induction signal from the core due to annual changes in the magnetosphere. These time-varying effects due to the variability of the solar wind density along Mercury’s orbit induce an internal dipole field, which depends on the size of the core. By assuming ideal magnetohydrodynamic conditions, changes in the solar wind density lead to a departure from the pressure balance law between shocked solar wind pressure outside, $p_{sw}$, and magnetic pressure inside the magnetopause, $p_{B}$. For an axial dipole, this balance at the front nose of the magnetopause is given by:

$$ \kappa p_{sw} (r_{hel})= \frac{1}{2} \kappa \rho v_{sw}^{2} = p_{B} = \frac{B^{2}(r=R_{MP})}{2\mu _{0}} = \frac{ \xi (g_{1}^{0})^{2} }{2\mu _{0} } \left (\frac{R}{R_{{MP}}} \right )^{6} $$

(17)

where $r_{hel}$ and $v_{sw}$ stand for the heliocentric distance and the solar wind speed, respectively. This leads to a sub-solar magnetopause standoff distance:

$$ R_{MP} \propto p_{sw}^{-1/6} (r_{hel}) . $$

(18)

A constant solar wind thermalization factor $\kappa \approx 0.88$ is assumed, and $\xi$ stands for a constant magnetospheric form factor (this must be less than 9), which depends on the magnetopause model in consideration. If the solar wind speed $v_{sw}$ can be considered constant over Mercury’s orbit, $p_{sw}$ varies only with the solar wind density. From continuity of the flow, we expect a dependence on the heliocentric distance $r_{hel}$ as $p_{sw} \propto r_{hel}^{-2}$. Thus, in absence of any induction the magnetopause standoff distance should vary as a function of heliocentric distance (and thus time) as follows:

$$ R_{MP} \propto r_{hel}^{0.33} . $$

(19)

However, with induction in the interior this relationship is expected to change. Johnson et al. (2016) found a different magnetospheric compression exponent of 0.29, which corresponds to a temporal change of $5\%\; (9.5\text{nT} \cdot R^{3})$ of the dipole moment from Mercury’s perihelion to aphelion. This result is slightly lower than the estimate presented by Grosser et al. (2004) probably because of differences in the modeled magnetopause distances and core sizes.

An alternative approach was discussed by Johnson et al. (2016) to the ratio between internal and external fields from orbital subsets at peri- and aphelion. To avoid an influence of field-aligned currents, only the radial magnetic field was used for the inversion. This different method showed a change in the internal axial dipole coefficient of 7.5 nT between perihelion and aphelion. By modeling the interior as simple nested shells of uniform conductivity (iron core: $\sigma =10^{6}$ S/m; silicate mantle+crust: $\sigma =0$ S/m), a simple transfer function (Rikitake 2012; Grosser et al. 2004) of induced $I_{\ell }^{m}$ ($g_{\ell }^{m}$ and $h_{\ell }^{m}$ in Eq. (10)) to inducing field $E_{\ell }^{m}$ coefficients ($G_{\ell }^{m}$ and $H_{\ell }^{m}$ in Eq. (10)) is given by:

$$ \mathcal{R}_{\ell }=\frac{I_{\ell }^{m}}{E_{\ell }^{m}} = \frac{\ell }{\ell +1} q^{2\ell +1} , $$

(20)

which depends on the core size in $q=R_{CMB} / R$. This expression does not explicitly depend on the order $m$ but also describes the induction effect of non-axisymmetric Gauss coefficients. If the whole planet was superconducting, the transfer function for dipolar fields would be $\mathcal{R}_{1} = 0.5$. The transfer function was found to be in the range $0.23 - 0.30$, and the iron core size from both approaches from Johnson et al. (2016) should be in the range $1900 - 2060$ km. This core size estimate is based on the magnetic data only, and it is fully consistent with the results obtained by independent geodetic investigations (Hauck et al. 2013; Genova et al. 2019).

A determination of the magnetopause crossings in the MESSENGER data was also retrieved by Zhong et al. (2015) with a refined modeling of the magnetopause. These results indicate a larger magnetospheric compression exponent of $R_{MP} \propto r_{hel}^{0.42}$. The data interestingly follows quite closely the classical $\propto r_{hel}^{0.33}$ relation at perihelion, and the deviation increases with distance to the Sun. This was interpreted as a higher reconnection efficiency (due to lower plasma $\beta$ and lower Alfvénic Mach number) at perihelion (Zhong et al. 2015). At the magnetopause of Mercury, this process connects magnetic flux tubes from the magnetosphere with flux tubes in the magnetosheath. The solar wind stream there pulls these tubes to the magnetotail, eventually. This leads to a flux erosion of the dayside magnetosphere. To infer the planet’s core size from this result, an accurate knowledge of the interplay between reconnection and induction at Mercury is necessary. Heyner et al. (2016) showed that the reduction in standoff distance during reconnection is enhanced by a neutral current sheet amplification, leading to a more complicated scenario. Observations of Mercury’s magnetic field from a single spacecraft, as MESSENGER, do not provide enough information to relate processes at the dayside with processes happening at the magnetotail. The dual-spacecraft configuration of the BepiColombo mission will be less affected by these limitations, enabling a better characterization of reconnection and induction mechanisms.

An independent estimation of Mercury’s core through the analysis of MESSENGER magnetic data was retrieved by Wardinski et al. (2019). By using an estimated average field model to isolate the time-varying magnetic fields around Mercury, a dipole transfer function of $R_{1} = 0.30$, which suggests a core size of $R_{CMB}=2060$ km, was obtained in agreement with the upper bound found by Johnson et al. (2016). The variation of the internal axial dipole field, on the other hand, is roughly 5 times stronger than the estimates provided by Johnson et al. (2016). At perihelion, the internal as well as the external field are amplified by the magnetic field of induced currents to withstand the compression effect by the increased solar wind pressure. Equation (18) shows that the induced dipole field suggested by Wardinski et al. (2019) would increase the sub-solar standoff distance by about $0.2 R$ at perihelion. This additional standoff distance change, however, was not observed and an explanation for such a high change in internal and external dipole fields has not been found yet. The accurate measurements of BepiColombo magnetometers will help answer these open questions to constrain the interior conductivity profile by using inductive methods.

Induction in the planet’s interior is also important on short time scales. Due to the shallower penetration depth of theses signals, this would eventually enable us to constrain the mantle conductivity structure. Jia et al. (2019) showed that the magnetic field induced in the core works against the compression of the magnetosphere by the temporally high solar wind pressure during a coronal mass ejection passage. During an extreme event, with a reconstructed solar wind ram pressure of 93 nPa, compared to the average values of $10 - 15$ nPa along Mercury’s orbit, the effective dipole moment of the planet was increased up to 246 nT $\cdot R^{3}$. The BepiColombo MPO and Mio spacecraft will enable accurate estimations of the time evolution of higher multipole degrees.

4.3.3 Crustal Magnetic Fields

Crustal magnetic fields may bear the record of dynamo evolution and/or true polar wander events. Planetary crust is magnetized if its constituents or iron-bearing minerals cooled below the Curie temperature in presence of a background magnetic field. The crustal field can be measured by a low-flying spacecraft. On the northern hemisphere, the crustal fields were only unambigously observed at altitudes below about 130 km. Therefore, it is expected that the crustal magnetic fields on the southern hemisphere are also only detectable at these altitudes. Johnson et al. (2018) performed a global analysis for crustal magnetic fields based on MESSENGER data with altitudes below 100 km from the low-altitude campaign of the extended mission. By removing the average magnetospheric field and (band-pass) filtering the data, only fields with spatial scales <520 km and temporal scales <1 Hz were included in this study. The strongest crustal magnetizations (with max. 0.4 A/m for a 10 km thick layer) are found at the Circum Caloris Plains and in the northern part of the Caloris Basin. The analysis by Johnson et al. (2018) of the observed magnetic field strengths over the Suisei Planitia region reveals that, depending on the thickness and the magnetic susceptibility of the crust, the ancient magnetizing field ranged between the present fields strengths of only a few hundred nT and Earth-like fields around $10,000$ nT.

Oliveira et al. (2019) studied 5 magnetic anomalies associated with impact craters and inverted for the locations of the paleopoles to gain insights into the magnetic history of Mercury. Magnetically northern paleopole locations were all found in the geographically southern hemisphere, but the uncertainties in the results are quite large. Further data from the BepiColombo spacecraft during the extended mission, which will enable lower spacecraft altitudes in the southern hemishere, may allow drawing a more complete picture of the crustal magnetization.

Another way to characterize the crustal magnetic fields is to analyze the magnetic field spectrum of the planet. At Earth, the relatively flat spectrum typical for crustal magnetization starts to dominate the spectrum at harmonic degree 14 (see, e.g., Glassmeier and Heyner (2021) for a recent overview). Johnson et al. (2018) estimate that the minimum scale of crustal magnetization at Mercury is <40 km and these magnetizations are only unambiguously visible at altitudes below 100 km, which corresponds to a spherical harmonic degrees in the range of $150< l < 380$. The contribution at the larger scales of typical models for the dynamo field is likely negligible, thus, there should not be a break in the magnetic field spectrum. In addition, at these small scales stated above, it will be challenging to correctly separate local external magnetic fields from crustal anomaly fields. However, the analysis by Hood (2016) and Hood et al. (2018) indicate much larger magnetic anomaly scales, which may be already detected at higher spacecraft altitudes during the BepiColombo nominal mission. A more detailed discussion on the resolvability of crustal magnetic field anomalies at Mercury may be found in e.g. Heyner et al. (2021) and Rothery et al. (2020).

5 Mercury’s Internal Structure

5.1 Crust

Mariner 10 observations suggested a volcanic origin for some of the smooth plains of Mercury, but the nature and extent, both spatially and temporally, of the volcanism remained unknown (Strom et al. 1975). Multispectral imaging, spectroscopy, and topography observations collected by the MESSENGER spacecraft showed that essentially the whole of Mercury’s surface is a secondary crust formed by partial melting in the mantle, ascent of the buoyant melt and deposition on or below the surface by extrusive and intrusive volcanism between 4.2 and 3.5 Gyr ago (see Byrne et al. (2018) for a recent review). The existence of some low reflectance material at the surface may be associated with a primary graphite floatation crust resulting from magma ocean crystallization (Vander Kaaden and McCubbin 2016). Because of the secondary nature of the crust, observations of the surface composition can be linked to the composition of the original source regions in the deeper mantle. The surface mineralogy could not be inferred directly from MESSENGER data, which do not show any diagnostic spectral absorption features (Byrne et al. 2014), but it was approximately determined from crystallization experiments and calculations on candidate compositions based on measured elemental compositions. The results show that the surface consists essentially of plagioclase, sulfides, and Mg endmember pyroxene and olivine (Stockstill-Cahill et al. 2012; Namur et al. 2016a; Vander Kaaden and McCubbin 2016). The mantle is dominantly lherzolitic, with rocks mainly made of olivine (forsterite) and orthopyroxene (enstatite) (Namur et al. 2016a; Vander Kaaden and McCubbin 2016). Because of the relatively low pressure in the thin mantle of Mercury compared to the other terrestrial planets, phase transitions to minerals as in the lower mantle of the Earth do not occur.

BepiColombo will significantly improve our knowledge of the surface and our understanding of the crustal formation processes by extending the data to the southern hemisphere, by providing higher spatial mapping resolution in equatorial regions, and by thermal infrared spectrometer observations (i.e., MERTIS, SIMBIO-SYS) capable of characterizing directly the surface mineralogy (see the paper by Rothery et al. of this issue). Complementary improved topography and gravity data from BELA (Sect. 3.1) and MORE (Sect. 3.2) investigations will further constrain the crustal density and crust–mantle interface. All the data combined will contribute invaluable information on the total volume of the crust, which represents the total amount of igneous processes over Mercury’s history. It provides, therefore, insight into the internal temperature evolution, the relative amount of heat producing elements in the crust compared to the mantle, and the thermal history of Mercury (Sect. 5.5).

Crustal density can be obtained independently from mineralogy and gravity. Surface compositional data are only representative of the outer layers of the crust, and, most importantly, do not account for porosity. Gravity has the advantage of providing an average estimate of the crustal density from the surface to the crust-mantle boundary. The building history of the crust is ill–constrained, but the surface heterogeneity is indicative of a complex formation with different phases of crust production. The best sensitivity to the local density of the crust is obtained with local gravity analyses techniques (e.g., Beuthe et al. 2006; Goossens et al. 2005; James et al. 2018). During the nominal BepiColombo mission phase, the much lower MPO spacecraft altitudes over the southern hemisphere will largely improve the gravity field (Sect. 4.2.1). The information on the density from gravity and spectroscopy can be used as input to interpret gravity anomalies in terms of the structure of the crust and lithosphere, density anomalies in the mantle, and topography of the core–mantle boundary (see Phillips et al. (2018) for a recent review). By subtracting the contribution of the topography (i.e., height above the reference ellipsoid) from the observed gravity, the resulting gravity signal (i.e., the Bouguer gravity) reveals the asymmetric mass distribution below the surface in the crust, mantle, and core. The highest sensitivity is to the mass anomalies closest to the surface at high gravity resolutions.

The combined use of global BepiColombo gravity and topography data will also inform on how surface topography can be maintained and on the degree of compensation in Mercury. Compensation indicates that the free–air gravity variations are weaker than those expected from variations in topography. With compensation, high topographic rises have less mass below the mountains than regions of average topographic height, and topographic lows have excess mass underneath. Several scenarios can lead to this equilibrium state, in particular, variations in crustal thickness (Airy isostasy) or regional variations in composition and density (Pratt isostasy). Airy isostasy, the dominant compensation mechanism on Earth, expresses that any topographic rise is associated with a crustal root in such a way that the pressure at the root bottom is equal to that at the same depth below the geoid but at another horizontal position without topographic relief. The mean crustal thickness has been determined based on the assumption of Airy isostasy by Padovan et al. (2015) and by Sori (2018) to be 35±18 km and 26±11 km, respectively. Significant peak-to-peak variations are also expected (James et al. 2015; Beuthe et al. 2020). These estimates were retrieved by excluding surface regions expected not to be in a state of Airy isostasy, such as large impact craters and smooth planes that might have been emplaced on a thick lithosphere. Furthermore, the assumptions of isostasy could not be fully proved on the remaining regions. Pratt isostasy is unlikely to be a major compensation mechanism (Sori 2018). Crustal thickness variations derived from the Bouguer gravity (Smith et al. 2012; Genova et al. 2013; Mazarico et al. 2014; Genova et al. 2019; Beuthe et al. 2020) are of the order of several tens of km and are generally not isostatic. Additional topographic support is most likely due to mass anomalies below the crust, with possible contributions from lithospheric flexure (James et al. 2015). The nature of these deep mass anomalies is unclear. It can include compositional variations, perturbations of an interface between layers of different density, and lateral variations in temperature. Mantle convection is unlikely to contribute, since the length scale of convection in the thin mantle is too small (James et al. 2015), and also the gravity effects are most likely too small to be detected from spacecraft (Padovan et al. 2015).

Our knowledge of Mercury’s crust is still limited to the northern hemisphere, where the MESSENGER spacecraft provided high resolution gravity and topography data. The global surface coverage through the MPO spacecraft will enable a better understanding of surface features and crust mineralogy, which will allow us to determine the mean thickness of the crust and elastic lithosphere, the crustal thickness variations and mantle structure. These geophysical results will help characterize the history and process of crust formation and its relation with mantle convection.

5.2 Core and Mantle

5.2.1 Obliquity

Compared to the other terrestrial planets, the mean density of Mercury is much larger than expected for its radius. With a similar bulk composition as Mars or the Earth, Mercury would be expected to have a mean density close to 4000 kg m⁻³, ∼25% smaller than the actual mean density of 5427 kg m⁻³. The large density implies that Mercury has a much larger fraction of heavy chemical elements than the other terrestrial planets. Since iron is the only heavy element sufficiently abundant in the Solar System, Mercury must be strongly enriched in iron. As for all Solar System bodies with a radius of at least a few hundred kilometers, the accretion process, and in particular large impacts at the end of Mercury’s formation led to an increase in temperature of at least a significant part of Mercury’s interior to above melting temperatures and to the formation of a magma ocean (Brown and Elkins-Tanton 2009; Elkins-Tanton 2012). The differentiation of the planet was then facilitated by the formation of the core through the sinking of siderophile elements with iron, and of the mantle formed by lithophile elements. The high mean bulk density indicates that the core must be relatively much larger than for the other terrestrial planets. Early studies reported estimates of 1860±80 km (e.g., Spohn et al. 2001) for the iron core radius, which is ∼3/4 of the radius of the planet. The existence of a global–scale magnetic field is an indirect evidence for the existence of Mercury’s iron–rich core, which has to be at least partially liquid and convecting to produce the dynamo action. It also suggests that the solidification of the inner core may have started, since the additional related buoyancy sources would help to maintain a dynamo to the present date (Schubert et al. 1988; Hauck et al. 2004, and see Sect. 5.4).

A further constraint on the radial density profile in Mercury, and thus also on the differentiation of Mercury into a core and silicate shell is given by the mean moment of inertia. For bodies in hydrostatic equilibrium, the moment of inertia can be determined from the degree–two gravity field and shape, but this method cannot be applied to Mercury because of its slow rotation. The flattening into an ellipsoid by rotation for a hydrostatic Mercury is two orders of magnitude lower than the observed flattening, confirming the expectations (Van Hoolst et al. 2008; Matsuyama and Nimmo 2009). Constraints on the core, therefore, require additional data to complement the low–degree shape and gravity field. The most important constraints are given by the rotation, tides, and magnetic field, all of which will be measured with unprecedented accuracy by the BepiColombo mission (Sects. 4.1.1, 4.2.1, and 4.3.1).

For the Earth, and also for Mars (Konopliv et al. 2016), the moment of inertia is very accurately determined from precession combined with the degree–two gravity field. The slow precession rate for those planets is inversely proportional to the polar moment of inertia $C$ (e.g., Smart 1953), and the degree–two gravity field provides two relations for the three principal moments of inertia (i.e., the equatorial principal moments of inertia $A$ and $B$ in addition to $C$), allowing to estimate the moment of inertia ${\mathrm{MOI}}=(A+B+C)/(3MR^{2})$. This method, however, cannot be used as such for Mercury because tidal dissipation drove Mercury to an equilibrium Cassini state on a time scale that is short compared to the age of the Solar System (e.g., Peale 1974; Genova et al. 2019). Tidal dissipation depends on the orbital distance $d$ to the Sun as $d^{-6}$ and is, therefore, strongest for Mercury that is the innermost planet of the Solar System. In this state, both the orbit normal and the rotation axis precess around the normal to the Laplace plane in ∼320 kyr (Yseboodt and Margot 2006). The Laplace plane is the plane about which Mercury’s orbit precesses with nearly constant inclination. Besides the slow precession, the rotation axis of Mercury also precesses about its position in the Cassini state if Mercury does not exactly occupy but is close to the Cassini state (Peale 1974; Ward 1975). If we neglect the much slower orbital precession, the precession due to the solar torque (also called free precession) can be compared to the precession of the Earth about the normal to the ecliptic, but is much faster with a period of ∼500 years (Peale 2005; Van Hoolst 2015). The final obliquity of the Cassini state is not zero, since planetary induced changes in the orbit of Mercury would cause obliquity to become non–zero again. The equilibrium obliquity depends on how well Mercury’s orientation can track the changing orbit, and, therefore, depends on the polar moment of inertia. It can be expressed as (Peale 1981; Baland et al. 2017):

$$ \epsilon = \frac{-C\dot{\varOmega }\sin i}{C\dot{\varOmega }\cos i + 2nMR^{2}G_{201}(e)C_{22} - nMR^{2}G_{210}(e)C_{20}}, $$

(21)

where $n$ is the mean motion, $i$ is Mercury’s orbital inclination with respect to the Laplace plane, $e$ the orbital eccentricity, $\varOmega$ the longitude of the ascending node of the orbital plane with respect to the Laplace plane ($\dot{\varOmega }$ is negative), and $G_{201}(e)$ and $G_{210}$ are Kaula’s eccentricity functions defined as

$$ G_{210}=(1-e^{2})^{-3/2}, $$

(22)

$$ G_{201}=\frac{7}{2}e - \frac{123}{16} e^{3} + \frac{489}{128}e^{5}+O(e^{7}). $$

(23)

Equation (21) shows that a smaller obliquity is associated with a smaller $C$ as expected. Mercury can more easily track the orbital precession when it has a smaller moment of inertia, since its free precession period is then shorter. The dimensionless polar moment of inertia $C/MR^{2}$ is largest and equal to 0.4 for a planet with a uniform density. Mercury’s obliquity has been determined based on Earth–based radar observations (Margot et al. 2012), MESSENGER imaging and laser altimeter data (Stark et al. 2015), and radio science investigation (Mazarico et al. 2014; Verma and Margot 2016; Genova et al. 2019; Konopliv et al. 2020) (Fig. 9). Surface–related measurements and the gravity solution HgM008 (Genova et al. 2019) suggest that Mercury occupies the Cassini state. These solutions are all close to 2 arcmin, but only marginally consistent. The resulting $C/MR^{2}$ derived from the estimated obliquity is between ∼0.32 and ∼0.35, although individual estimates are more precise with errors down to ∼1% (Genova et al. 2019). For comparison, Mercury’s dimensionless polar moment of inertia is much less accurately known than that of Mars, which has an uncertainty of ∼0.1% (Konopliv et al. 2016). Gravity measurements of the BepiColombo mission will allow achieving this level of accuracy for $C/MR^{2}$ by determining the obliquity with an uncertainty of <1 arcsec (Sect. 4.2.3).

Fig. 9

Pole’s coordinates right ascension ($\alpha$) and declination ($\delta$) in the International Celestial Reference Frame at epoch J2000. Square–markers show the estimates of the pole’s coordinates retrieved by direct surface measurements (Margot et al. 2012; Stark et al. 2015), and circle–markers the results from different gravity solutions (Mazarico et al. 2014; Verma and Margot 2016; Genova et al. 2019; Konopliv et al. 2020). The black dashed line is the Cassini state constraint. Mazarico et al. (2014), Verma and Margot (2016), and Konopliv et al. (2020) uncertainties are reported as dashed lines since those studies only reported $\alpha$ and $\delta$ formal uncertainties that were scaled by a certain factor. The error ellipses of the solutions by Margot et al. (2012), Stark et al. (2015), and Genova et al. (2019) account for the correlation amongst these two parameters providing a 95% confidence uncertainty region

Because of the precession of Mercury’s pericenter, Mercury’s orientation will slightly deviate periodically from coplanarity (Peale 1974). The obliquity and instantaneous precession rate, therefore, vary with time. The deviation with respect to coplanarity is of the order of 1 arcsec and varies periodically with a period of about 84 kyr (Peale et al. 2014, 2016; Baland et al. 2017). At the current epoch, the pericenter precession leads to a small change in obliquity of ∼0.2 arcsec and a deviation with respect to the Cassini state of ∼0.9 arcsec (Baland et al. 2017). The deviation is measured with respect to the Cassini state line, which is the intersecting line of the orbit plane with the plane formed by the coplanar spin vector, the orbit normal, and the normal to the Laplace plane. In addition, tidal deformations of Mercury slightly modify Mercury’s angular momentum and decrease the average gravitational torque of the Sun on Mercury, which drives the spin precession. The main effect of the tides is to introduce a constant shift in the mean obliquity of ∼0.4 arcsec and a time-varying deviation of <3 arcsec (Baland et al. 2017). A departure of Mercury’s spin axis from the stationary Cassini state is then expected. However, these predicted deviations are well below the current uncertainty level of the measured obliquity. BepiColombo measurements will enable a highly accurate estimate of the obliquity, and its geophysical interpretation will need to account for deviations due to the effect of the tides. A precise measure of the obliquity will also lead to constrain the ratio $k_{2}/Q$, which informs on the mantle rheology since $1/Q$ is the dissipation factor (Kaula 1964). Our current knowledge of the obliquity already puts an upper limit of 0.02 on Mercury’s $k_{2}/Q$ (Baland et al. 2017).

Equation (21) assumes that the core follows the mantle during the very slow orbital precession with a period of ∼320 kyr (Peale 1974). Peale et al. (Peale et al. 2014) showed that pressure coupling of the core fluid on the degree–two shape of the core–mantle boundary drives the spin axis of the mantle to a direction that deviates from the classical Cassini state by <0.1 arcsec, which lies even below the BepiColombo observational precision. A Poincaré flow is assumed to model deviations of the core flow from its corotation with the mantle (Poincaré 1910). The spin axis of the core is also nearly on the Cassini state line, but can be at an angular distance of several arcmin from the mantle spin axis. This deviation depends on the total amount of core–mantle boundary coupling. Besides the pressure torque of the core fluid on the ellipsoidal shape of the core-mantle boundary, also the dissipative viscous, electromagnetic, tidal and topographic torques contribute to the total torque. The dissipative core–mantle boundary torques can also lead to a deviation of the core spin from the Cassini state by several arcsec (Peale et al. 2014). Although the spin axis of the core can be slightly misaligned with that of the mantle, the estimates of the obliquity from surface observations (e.g., imaging and altimetry) would give the same value as from gravity measurements if the symmetry axes of the core–mantle boundary are aligned with those of the surface. Therefore, if the core of Mercury were to be entirely liquid, independent estimates of the obliquity from surface-related and gravity measurements would be fully consistent, and Eq. (21) could be used in both cases to determine the moment of inertia of the whole planet.

A solid inner core would lead to different obliquities of the core and mantle spin axes, which would deviate from the classical $\epsilon$ of Eq. (21) (Peale et al. 2016). Peale et al. (2016) suggest that the inner core have a larger obliquity than the silicate shell, inducing gravitational coupling between the silicate shell and the core. The obliquity of the silicate shell then tends to be larger than the value predicted by Eq. (21). Larger inner core sizes yield larger silicate shell obliquities. An inner core of $\sim60\%$ of the size of the total core leads to an increase of ∼10 arcsec, which is ∼5 times larger than the uncertainty of current obliquity estimates. Therefore, the presence of a solid inner core could be responsible for the different $\epsilon$ estimates obtained with MESSENGER imaging and altimetry data (Stark et al. 2015), which are sensitive to the silicate shell only, and gravity data (Genova et al. 2019), which measure the spin axis of the whole planet. In the interior models considered by Peale et al. (2016), the use of the classical obliquity equation (21) for the mantle spin can, therefore, lead to an overestimation of the normalized moment of inertia $C/(MR^{2})$ by up to ∼0.02-0.03. This conclusion, however, depends on the shape of the inner core boundary and of the core–mantle boundary, and applies only to the case that Mercury’s interior is in hydrostatic equilibrium and behaves elastically. Deriving the moment of inertia of Mercury from the observed obliquity requires knowledge of the size and shape of the inner core. This makes it more difficult to constrain the interior of Mercury from the observed obliquity, but also offers means to put constraints on the inner core. Independent estimates of the obliquity of the rocky shell from BELA data (Sect. 4.1.3) and SIMBIO-SYS images (see Rothery et al. in this issue), and the mean planet obliquity from the gravity field (MORE, Sect. 4.2.3) will enable a better characterization of Mercury’s interior by disentangling the contributions of the core and the outer layers.

5.2.2 Libration

Mercury experiences a gravitational torque exerted by the Sun that periodically reverses as a result of the 3:2 spin–orbit resonance. Mercury’s spin rate increases and decreases with a period equal to the orbital period of 88 days. Those longitudinal librations have been determined to be between ∼37 arcsec and 40 arcsec (or 435 m and 475 m as measured on the equator), based on independent observations from the Earth–based radar (Margot et al. 2012), the MESSENGER camera and laser altimeter, (Stark et al. 2015) and radio science (Mazarico et al. 2014; Genova et al. 2019). If Mercury were to be entirely solid, the whole planet would participate in the libration and the amplitude of libration would be proportional to the moment of inertia of the whole planet and equal to ∼190 m. The observationally determined value is more than twice as large (Margot et al. 2007), demonstrating that only part of the planet participates in the motion. Therefore, libration provides convincing evidence that Mercury’s core must be liquid, at least partly (Margot et al. 2007, 2012). With an entirely liquid core, only the mantle performs the librational motion (Peale 1974; Van Hoolst et al. 2012). Libration then gives a measure of the moment of inertia of the silicate shell (i.e., $C_{cr+m}/C$). This allows distinguishing between contributions to the moment of inertia from the shell and from the core, and puts strong constraints on the size of the entire core. Different studies support that Mercury’s core radius must be close to 2000 km with an uncertainty of several tens of km (Hauck et al. 2013; Rivoldini and Van Hoolst 2013; Knibbe and van Westrenen 2015; Margot et al. 2018; Genova et al. 2019). Remarkably, Mercury’s core radius is better known than that of Mars, a planet that has been investigated by many more spacecraft but lacks direct information on the mass distribution in the core.

The existence of an inner core complicates the interpretation of the libration amplitude in terms of the moment of inertia of the shell. Gravitational coupling between the inner core and the shell decreases the libration amplitude, as expected since the inner core also librates (Van Hoolst et al. 2012). The libration amplitude decreases for larger inner cores by up to about 20 m (i.e., 2 arcsec) and is almost independent of the flattening of the inner core (Van Hoolst et al. 2012; Dumberry et al. 2013). Libration offers the possibility to constrain the size and density of the inner core. In addition to studying the main libration at 88 days, information on the interior can be gained from libration related to orbital changes caused by other planets (Peale et al. 2007; Dufey et al. 2008; Yseboodt et al. 2013). Those orbital perturbations and resulting librations have periods commensurate with the orbital periods of the perturbing planets and can be resonantly amplified when the periods are close to those of free librational modes of Mercury. A significant perturbation is induced by Jupiter with a period of 11.86 yr (Yseboodt et al. 2010). A solid Mercury has one free librational mode with a period close to 18 years. With an entirely fluid core, the mode describes a long–period oscillation of the axis of minimum moment of inertia of the solid outer shell about the Mercury–Sun line at perihelion and has a period close to 11 years (Peale 2005). With a solid inner core, the period is between those two values. A solid inner core also introduces a second free librational mode, in which the inner core and mantle oscillate out of phase, with a period expected to be between 3 and 5 years (Van Hoolst et al. 2012; Dumberry et al. 2013). The free periods can be close to the period of the planetary perturbations and lead to a resonant enhancement of the corresponding forced libration to observable values (Yseboodt et al. 2013). Observation of such a libration will give information on the inner core. The free librations would also provide valuable information on the inner core, but would require a recent excitation (e.g., Koning and Dumberry 2013).

MESSENGER surface–related (Stark et al. 2015) and gravity measurements (Genova et al. 2019) determined a higher mean rotation rate than that of the exact 3:2 resonant rotation, an offset that could be caused by a long–period libration. Identification of the period and the amplitude is needed to infer properties of the core and requires observations over several years. This might be possible with BepiColombo, which will significantly extend the time series of libration observations. The combination of MESSENGER and BepiColombo datasets will then be fundamental to retrieve long–term variations of Mercury’s rotational state.

5.2.3 Tides

Tides are excellent probes for a liquid planetary core. The Love number $k_{2}$ is currently the main constraint on the Martian core (Yoder et al. 2003; Rivoldini et al. 2011; Genova et al. 2016; Konopliv et al. 2016) and the only direct indication that the core of Venus is at least partially liquid (Konopliv and Yoder 1996). Likewise, the precise determination of the Love numbers $h_{2}$ (Sect. 4.1.2) and $k_{2}$ (Sect. 4.2.2) by BepiColombo will put strong constraints on Mercury’s interior. Thanks to its proximity to the Sun, Mercury’s tides are the largest of the Solar System planets, with a tidal surface displacement of the order of 2 meter (Van Hoolst and Jacobs 2003; Rivoldini et al. 2009; Padovan et al. 2014). The largest tidal waves are sectorial and have periods equal to the orbital period of Mercury (i.e., equal to half a solar day on Mercury) and half that value. The sectorial wave at the orbital period is about 3 times larger than the zonal wave at the same period, which is the third–largest tide (Van Hoolst and Jacobs 2003). The effect of tides on the external gravitational field of Mercury has already been determined with an accuracy of about 5% by analyzing radio science data of the MESSENGER mission. Early estimates based on data before the extended mission low-altitude campaign suggested $k_{2}=0.451\pm 0.014$ (Mazarico et al. 2014) and $k_{2}=0.46\pm 0.02$ (Verma and Margot 2016). The analyses of the entire MESSENGER dataset including data at orbital altitudes lower than 100 km, led to larger independent estimates $k_{2}=0.57\pm 0.03$ (Genova et al. 2019) and $k_{2}=0.53\pm 0.03$ (Konopliv et al. 2020).

The latest estimates of the Love number $k_{2}$ are about an order of magnitude larger than for Mercury models with an entirely solid core and, therefore, confirm that Mercury’s core is at least partially liquid. Tides depend mainly on the size of the core. However, other interior properties, including core and outer layers densities, (an)elastic properties of the mantle, and inner core size and (an)elastic parameters, provide a significant theoretical spread in $k_{2}$ and $h_{2}$ values of ∼0.15 (Van Hoolst and Jacobs 2003; Rivoldini et al. 2009; Padovan et al. 2014). Tidal Love numbers increase with increasing core size, decreasing mantle rigidity, and increasing mantle temperature (e.g., Margot et al. 2018). Since $k_{2}$ and $h_{2}$ depend differently on these geophysical properties, joint estimates of these Love numbers can help disentangling the different contributions. The MESSENGER mission only enabled the recovery of the gravitational tides, limiting this geophysical interpretation. The analysis of BELA and MORE data will allow us to precisely adjust both surface deformation and gravitational tides. By combining the Love numbers $h_{2}$ and $k_{2}$ formal uncertainties that will be determined by the BepiColombo mission, Steinbrügge et al. (2018) showed a method to significantly constrain the size of the inner core if its radius is larger than ∼700 km. Furthermore, the determination of the tidal phase–lag with an accuracy better than ∼0.5^∘, which will be fulfilled by the MORE investigation (Table 3), provides information on the temperature at the core–mantle boundary and the grain size in the mantle.

5.3 Interior Modeling

Planetary geodesy and geophysics provide fundamental constraints on Mercury’s core and silicate shell, but require theoretical interior modeling to tie the observational data with the properties of the internal structure. Obliquity and libration give mainly information on the polar moment of inertia of the planet and the mantle, depending to a lesser extent on the moments of inertia and flattening of the inner and outer core. Love numbers also inform on the mass distribution and the layers physical state (i.e., solid or liquid), and on the deformation. The geodetic data are thus related to the internal density and (visco-)elastic properties of the planetary solid and liquid materials. In the geophysical interpretation of the data, the theoretical interior modeling must be able to self-consistently describe radial profiles of density and other material properties, and to distinguish between solid and liquid layers. Adequate equations of state describing the behavior of the constituent materials, based on experimental data and theoretical results, are, therefore, needed.

Initial results of Mercury’s interior modeling (e.g., Stevenson et al. 1983; Van Hoolst and Jacobs 2003; Hauck et al. 2007; Rivoldini and Van Hoolst 2013) assumed that the core consists of sulfur as the only light element in addition to iron. Sulfur strongly reduces the melting temperature of an iron alloy and can easily explain why Mercury’s core is still (partially) liquid. Because of its high solubility in metal and high abundance in the Solar System, it is also considered to be the main light element in the core of Mars, which has similar pressure and temperature conditions as Mercury. Silicon can also have a high solubility in molten iron over an extended pressure range. Fe-S and Fe-Si alloys are considered the two end-member compositions of Mercury’s core. At oxidizing planetary formation conditions and up to moderately reducing conditions, S behaves siderophile and strongly prevents Si from entering the Fe liquid, in which case S is expected to be the most abundant light element (Cartier et al. 2020; Rivoldini and Van Hoolst 2013). Mercury’s surface composition, however, indicates a highly reducing formation history (Nittler et al. 2011; McCubbin et al. 2012). Under such conditions, Si enters significantly the metal phase and strongly decreases the dissolution of S in the iron alloy. The resulting higher Si and lower S abundances yields models that agree with the inferred moments of inertia (e.g., Mann et al. 2009; Malavergne et al. 2014; Chabot et al. 2014b; Namur et al. 2016b; Margot et al. 2018). Recent melting experiments suggest that the sulfur concentration is at most 1.5 wt$\%$ in the core (Namur and Charlier 2017; Cartier et al. 2020). Higher concentrations of Si are required to maintain Mercury’s core at least partially liquid because of its smaller effect on the melting temperature (e.g., Hauck et al. 2013).

Fe-S-Si alloys show an immiscible behavior at the pressure and temperature conditions in the top part of Mercury’s core (Sanloup and Fei 2004; Morard and Katsura 2010). If significant amounts of S and Si are dissolved in Mercury’s core, immiscibility can occur and result in a core structure with an S-rich alloy on top of a Si-rich alloy. The combined geochemistry and geodesy data, however, suggest that the light elements concentrations in Mercury’s core are too small for immiscibility to occur, making a fully mixed core most likely. This also suggests that a solid FeS layer at the core–mantle boundary, previously predicted by Smith et al. (2012), is unlikely to have formed by exsolution from the Fe-Si rich alloy (Rivoldini and Van Hoolst 2013). In an alternative scenario, an FeS layer could have formed during planetary differentiation from sulfur that does not dissolve into the metallic core or the silicate mantle (Malavergne et al. 2014; Chabot et al. 2014b). For bulk planetary sulfur concentrations below the chondritic S abundance, such a layer can at most be 90 km thick (Namur and Charlier 2017). By using Ti as a tracer of a sulfide melt formation, Cartier et al. (2020) showed that the FeS layer could have a maximum thickness of only 13 km and most likely is completely absent.

Information on Si thermodynamic properties are more limited than for S (e.g., Terasaki et al. 2019). First results show that the likely radius range for the core, estimated from the crust and mantle moment of inertia, only weakly depends on which light elements are present in the core (Hauck et al. 2013; Knibbe and van Westrenen 2015; Margot et al. 2018; Genova et al. 2019). The concentration of a given light element differs significantly between chemical elements because of their different effects on the density of an iron alloy. Obliquity and libration measurements allow then constraining different light element concentrations in Mercury’s core for different chemical elements. As Fe-Si alloys are denser than Fe-S alloys for a given concentration, more Si is needed than S to explain the geodesy data. Whereas cores with sulfur as the only light element are thought to have a bulk light element concentration of about 3 to 8 wt$\%$, silicon needs to be more abundant at 6 to 15 wt$\%$ if it is the only light element (Rivoldini and Van Hoolst 2013; Knibbe and van Westrenen 2015; Margot et al. 2018).

The existence, size, and composition of a possible inner core also strongly depend on the core composition. Whereas S does almost not partition in the solid, Si partitions nearly equally between solid and liquid at the core pressure and temperature of Mercury (Li et al. 2001; Kuwayama and Hirose 2004; Zhang and Fei 2008). If Mercury’s core consists of Fe–Si, its possible inner core would be significantly less dense compared to an Fe–S core. The two end-member compositions have different liquidus temperatures and, therefore, yield different sizes for the inner core. Larger solid inner cores are predicted for a pure Fe–Si alloy than for Fe-S (Knibbe and van Westrenen 2018; Rivoldini et al. 2018). Although the absence of an inner core for both compositions cannot be fully excluded (Rivoldini and Van Hoolst 2013), this scenario seems less likely with the latest obliquity estimate of the HgM008 gravity model (Genova et al. 2019). In general, an inner core has only a small effect on the moments of inertia, which are integrals over the product of density and a distance from the center to the power four. Therefore, the presence of a solid inner core may affect the polar moment inertia if a density contrast between inner and outer core is induced by significant different light element compositions.

Additional information on the inner core can also be obtained from Mercury’s observed radial contraction. Cooling of Mercury and the volume decrease that would be caused by a complete solidification of the core would result in a radial contraction of Mercury of above 20 km (Solomon 1976; Van Hoolst and Jacobs 2003; Grott et al. 2011). Early estimates based on Mariner 10 data of less than 0.5 km of contraction (Watters et al. 2004) have been continuously refined upward as more and more images of the planet’s surface obtained under better illumination conditions became available (e.g., Watters et al. 2009; Di Achille et al. 2012). The latest estimates based on $\sim 6000$ structures retrieved from MESSENGER images suggest an overall contraction of up to 7 km (Byrne et al. 2014), with a rate that may have declined over time (Crane and Klimczak 2017). In addition, up to two more km may have been accommodated by elastic deformation prior to the formation of the observed faults (Klimczak 2015). Thermal evolution models based on earlier estimates of 3–5 km of contraction or even less, suggested that convection in Mercury’s mantle could have stopped around 500–1000 Myr ago (Tosi et al. 2013). After the contraction estimates have been revised upward, this does no longer appear to be the case. These controversial planetary contraction results set an upper bound on the size of the inner core, depending sensitively on the composition and the evolution of temperature in the core and mantle. While the lowest estimate is incompatible with cooling histories of Mercury, the largest estimate would only allow a maximum inner core radius of ∼1000 km (Dombard and Hauck 2008; Grott et al. 2011; Hauck et al. 2004, 2013; Tosi et al. 2013; Dumberry and Rivoldini 2015; Knibbe and van Westrenen 2015).

5.4 Core Dynamo

Models of Mercury’s dynamo should at least strive to explain the three main field characteristics that have been shown by independent analyses of the MESSENGER data (Sect. 4.3.1). Important constraints on the dynamo models rely on models of Mercury’s internal structure, material properties, and thermal evolution discussed in Sect. 5.5. A classical prediction of the magnetic field strength in a dynamo is based on the force balance in the Navier-Stokes equation, which governs the flow dynamics. The magnetic field strength should saturate on a level where Lorentz forces are strong enough to sufficiently modify the flow. This would require that the Lorentz force is comparable to the Coriolis force, which dominates the Navier-Stokes equation. By using the Elsasser number $\varLambda$ as a simple estimate for the ratio of Lorentz to Coriolis force,

$$ \varLambda = \frac{\text{Lorentz Force}}{\text{Coriolis Force}} \approx \frac{B_{rms}^{2} \sigma }{\rho \varOmega } , $$

(24)

we would thus expect $\varLambda \approx 1$.

Recent more in-depth studies indicate that the Lorentz force effects on the flow may be more complex and are not well captured by the simple Elsasser number estimate. Theoretical considerations and dynamo simulations suggest that the rms magnetic field amplitude $B$ and the rms flow amplitude $U$ are primarily determined by the available convective power (Christensen and Aubert 2006; Christensen 2010; Davidson 2013; Wicht and Sanchez 2019). The respective scaling laws indicate that the magnetic energy in a planetary dynamo scales like

$$ B^{2}\big/(\rho \mu _{0}) \sim \left (R P\right )^{2/3} , $$

(25)

while the kinetic energy scales like

$$ U^{2} \sim \left (R^{1/3} P^{4/3}\big/\varOmega \right )^{2/3} . $$

(26)

Here $\rho$ is the mean core density, $\mu _{0}$ the magnetic permeability, $R$ the planetary radius, $P$ the available convective power per mass, and $\varOmega$ the rotation rate. The combination of both laws allows to eliminate $P$ and yields

$$ \varLambda \sim \text{Rm}^{3/2} E_{\lambda }^{1/2}, $$

(27)

where $\text{Rm}$ is the magnetic Reynolds number

$$ \text{Rm} = U\sigma \mu _{0} R_{M}, $$

(28)

and $E_{\lambda }$ the magnetic Ekman number

$$ E_{\lambda }=\left (\sigma \mu _{0} \varOmega R^{2}\right )^{-1}. $$

(29)

By assuming an electrical conductivity of $\sigma =10^{6}$ S/m and a rotation rate of $\varOmega =1.24 \times 10^{-6}$ s⁻¹, the magnetic Ekman number for Mercury is ${\sim} 10^{-7}$.

For an active dynamo, the magnetic Reynolds number must exceed unity since it measures the ratio of magnetic field production to dissipation. Numerical simulations suggest that $\text{Rm}>50$ is a more realistic requirement for planetary dynamos (Christensen and Aubert 2006). The scaling relations derived by Christensen (2010) indicate that $\varLambda$ is larger than $10^{-1}$ for dynamos dominated by the axial dipole component. Weaker multipolar dynamos, on the other hand, have a lower bound of $\varLambda \approx 4\times 10^{-2}$. By assuming a mean outer core density of $\rho =7000$ kg m⁻³, Mercury’s radius of $R=2440$ km, and the axial dipole amplitude as a proxy for the field strength, $\varLambda$ reaches only $5\times 10^{-6}$.

There are several possible reasons why the axial dipole component of the surface field may not be a good proxy for the typical field that should enter the force balance. Simple downward continuation to the top of the core would increase the dipole field by a factor of two. Higher order harmonics and the toroidal field, which cannot be measured at the surface, may actually be stronger than the axial dipole in the dynamo region. Moreover, special configurations like a particularly small (Heimpel et al. 2005) or a particularly large inner core (Stanley et al. 2005) may support the possibility that the field deeper in the dynamo region is notably larger than the field at the core-mantle boundary. Dynamo simulations suggest that the difference may amount to about an order of magnitude, possibly somewhat larger in extreme cases. This would bring the Elsasser number to $\varLambda \approx 10^{-3}$ at best, which is still at least one order of magnitude too small.

Simulations for standard dynamo setups geared to explain the geomagnetic field show that the axial dipole component dominates until the convective driving (Rayleigh number) exceeds a certain value. Beyond this threshold, the axial dipole loses its prominent role and the field becomes multipolar, has high dipole tilts and prominent higher degree field contributions and is also more time-dependent. The transition is best quantified in terms of the local Rossby number

$$ Ro_{\ell }= \frac{U}{\varOmega d}, $$

(30)

where $d$ is a typical flow length scale. Christensen and Aubert (2006) found that the threshold lies at $Ro_{\ell }\approx 0.1$, independent of the other system parameters. Estimates for Mercury suggest that $Ro_{\ell }$ is of order 10 for Mercury (Olson and Christensen 2006), which puts the planet’s dynamo safely in the multipolar regime. This is not in agreement with the observed dominant large scale field and the low dipole tilt.

5.4.1 The Feedback Dynamo

One way of explaining Mercury’s extraordinarily weak magnetic field is the special situation where the solar wind generated field and the internal dynamo process couple in a negative feedback loop. The planet is subject to an intense solar wind and thus relatively strong external fields. Chapman-Ferraro currents flowing within the magnetopause (the outer boundary towards the solar wind) produce a magnetic field that is anti-parallel to the internal dipole field at the core-mantle boundary. The long-time components of this external field may diffuse deeper into the core and interact with the dynamo (Glassmeier et al. 2007). Because of the anti-parallel configuration, the external field will reduce the internal main field and thereby lead to a less efficient dynamo process, the feedback dynamo.

When the convection in the core is driven very weakly and the dynamo starts off with a very weak seed field, for example after an impact or a field reversal, the feedback mechanism can keep the dynamo in the weak configuration (Heyner et al. 2010, 2011a,b). A relative external field strength of $\sim10\%$ is required to yield this effect. For stronger convective driving, the external field appears to influence the dynamo even more easily. In the simulations by Gómez-Pérez and Wicht (2010), the external field continuously decreases the core field until a reversal is triggered. Depending on dynamo parameters, only 2 % external field contribution may suffice to initiate this process. After the dipole reversal, the magnetosphere would reconfigure within $\sim2$ min (much faster than typical dynamo time scales) and re-establish the anti-parallel magnetic field configuration. This situation would in turn induce another reversal. Its influence on the overall $\varLambda$ in the core on longer time scales is still highly uncertain.

Another type of feedback dynamo has been proposed by Vilim et al. (2010). It seems conceivable that Mercury’s core temperature profile crosses the core material solidus twice somewhere in the middle of the core. A stably stratified iron snow region would develop between the two crossings, with a convective region above and a second one below. Both convective regions could harbor dynamos, which tend to favor the anti-parallel configuration, ultimately leading to a weak overall field.

Potential problems of the feedback dynamos are the unrealistic quadrupole-to-dipole (magnetospheric feedback) or octopole-to-dipole (Vilim et al. 2010) ratios. Figure 10 compares the Mauersberger-Lowes spectrum (Eq. 11) for different dynamo models and three different field models and illustrates the excessive quadrupole and/or octopole contributions in the feedback. In addition, the dipole tilt and non-axisymmetric field contributions are unrealistically large for the model by Vilim et al. (2010). Since only a few setups and parameter combinations have been considered, it seems conceivable that further numerical experiments would yield more realistic field configurations. However, there are certain limitations. While, for example, a stronger convective driving will boost the quadrupole to dipole ratio to more Mercury-like values, it will also lead to a higher time dependence, too large dipole tilts and generally stronger non-axisymmetric field contributions (Wicht and Heyner 2014).

Fig. 10

Mauersberger-Lowes spectrum normalized to the dipole part for the three different field models of Anderson et al. (2011), Thébault et al. (2017) (model: M2) and Wardinski et al. (2019) (model: MBF_a-n). These are compared to different dynamo models given in the legend. Thick lines represent field models and thin lines present results from dynamo simulations

5.4.2 A Stably Stratified Layer in Mercury’s Core

An alternative explanation of Mercury’s low field strength is a stably stratified layer in the outer part of the core, just underneath the core-mantle boundary, as sketched in Fig. 11. Several scenarios, which assume different composition and temperature of the core, would support the development of this layer (Sect. 5.3). Since Mercury’s core is still partially liquid, it must contain sufficient amounts of light constituents to reduce the melting temperature of the iron alloy. Oxygen, sulfur or silicon are likely candidates for planetary cores, but since Mercury seems to have formed under highly reducing conditions, oxygen can be ruled out, and also sulfur can only be present in small amounts (Cartier et al. 2020; Namur and Charlier 2017; Knibbe and van Westrenen 2018). Smaller amounts of carbon and hydrogen also seem to be possible. Consequently, large amounts of silicon of the order of 10 wt% are required to keep the core partially liquid. Such a high Si concentration is also required to yield a realistic core density in the absence of sulfur.

Fig. 11

Possible interior structure of Mercury with a stably stratified liquid core layer

The light core constituent can play an important role in driving the dynamo. Some light constituents are not readily incorporated into the solid phase, are at least partially released upon core solidification and give rise to compositional convection (as opposed to thermal convection driven by temperature differences). Since silicon partitions equally into the solid and the liquid phase, it would not support compositional convection. However, Si may increase the latent heat release (Desai 1986) and can help to power the dynamo in this way. Sulfur, on the other hand, promotes compositional convection since it partitions less into the solid phase.

The heat flux allowed to escape the core is fixed by the heat flux through the sluggish mantle. Convection in Mercury’s mantle is currently confined to a thin lower layer or may have stopped altogether. The heat flux through the core-mantle boundary is, therefore, small and likely lies below the adiabatic flux. This means that the upper part of the core is thermally stably stratified (Christensen and Wicht 2008; Knibbe and van Westrenen 2018), though the thickness of such a layer is not well constrained. Compositional convection based on sulfur, however, can potentially destabilize the stably stratified layer (Manglik et al. 2010).

Another complication could arise when the Fe-Si-S alloy is not miscible but rather separates into two immiscible FeSi and FeS constituents. Whether this really happens under Mercury’s conditions remains unclear (Morard and Katsura 2010; Margot et al. 2018). Since FeS is lighter than FeSi, it could accumulate underneath the core-mantle boundary (Hauck et al. 2013; Knibbe and van Westrenen 2015). If furthermore this FeS becomes solid, it would form a compositionally distinct stably stratified layer (Pommier et al. 2019). It cannot be excluded that liquid FeS would also form a stratified layer, but this would depend on the formation history and dynamics, both of which are unclear.

A third option for forming such a layer is iron snow. When the core adiabat crosses the melting curve at the core-mantle boundary, iron crystals form and snow towards the hotter interior, where they remelt. This leaves a lighter residual behind and thus forms a chemical stratified region that grows inward over time (Wicht and Heyner 2014; Breuer et al. 2015). The remelting iron crystals drive compositional convection in the underlying region. The latent heat released when the ‘snow flakes’ form, on the other hand, is released in the stable layer close to the core-mantle boundary and, therefore, cannot contribute to driving the dynamo. Recent ab-initio simulations suggest a core adiabat that seems to support this scenario (Edgington et al. 2019). On the other hand, a large inner core that occupies 50% of the core radius, as suggested by a recent analysis of geodetic data (Genova et al. 2019), strongly suggest that core solidification started at the center while core-mantle boundary solidification kicked in later (Wicht and Heyner 2014). This seems not compatible with the adiabat suggested by Edgington et al. (2019).

Christensen (2006) and Christensen and Wicht (2008) were the first to explore the effect of a stably stratified layer in dynamo simulations. They assumed a purely thermal driving and stratification with an inner core that occupies 50% in radius, which is in accordance with recent interior models (Knibbe and van Westrenen 2018; Genova et al. 2019). When combining a multipolar dynamo ($Ro_{\ell }>3$) with a stratified layer in the outer 28% in radius, the axial dipole reached a realistic amplitude. Dipole tilt and quadrupole to dipole ratio were Mercury-like at times, but also varied strongly during the simulation. Non-axial-dipole contributions remain often too prominent (Wicht and Heyner 2014). A comparison of the surface field with the field inside the dynamo region revealed that the stratified conducting layer acts like a filter that damps the more strongly time-varying field components with a magnetic skin effect. The time scale of spectral field contributions is roughly inversely proportional to the spherical harmonic degree: $\tau _{\ell }\sim \tau _{0}/\ell$ (Christensen and Tilgner 2004; Lhuillier et al. 2011). Smaller scale contributions (higher $\ell$) are therefore stronger damped. The exception is the axial dipole contribution, however, which typically has a much longer time scale than this simple rule would predict. Earth’s axial dipole varies on time scales of millennia, while higher harmonics have time scales of centuries to decades. Azimuthal flows in the stably stratified layer, driven by lateral temperature difference (thermal winds) or by Lorentz forces, can increase the damping of non-axisymmetric components. These considerations indicate that the stratified layer is an important ingredient of Mercury’s dynamo.

5.4.3 Modeling of a Persistent Offset

Three different mechanisms have been proposed for making the quadrupole to dipole ratio and possibly the ‘dipole offset’ more permanent. Two rely on a specific heat flux pattern that Mercury’s mantle imposes on the core. Cao et al. (2014) assumed a heat flux pattern for their numerical model, which allows more heat to escape at low than at high latitudes. This does not break the north-south symmetry per se, but, nevertheless, promotes a convective configuration with much stronger flows in the northern than in the southern hemisphere. Not surprisingly, dynamo action is also stronger in the north than in the south, which results in a persistent Mercury-like quadrupole to dipole ratio. However, the authors did not include a stably stratified layer in their simulation and, therefore, could not explain other characteristics of Mercury’s magnetic field. For a comparison with other dynamo models see Fig. 10.

Tian et al. (2015) show that a combination of a sufficiently thick stratified layer (35% of the core radius) and a core-mantle boundary heat flux pattern where more flux escapes through the north than the south yields very favorable results. Not only is the field strength very Mercury-like, offset values and tilt also remain very realistic at nearly all times.

Somewhat problematic for the models by Cao et al. (2014) and Tian et al. (2015) is the fact that neither core-mantle boundary heat flux pattern seems realistic. Mantle convection models rather predict a small scale pattern. If mantle convection has stopped operating long ago, the time averaged insolation pattern could have diffused down to the core-mantle boundary in $\sim500$ Myr (Tosi et al. 2015). This would suggest a larger flux at high latitudes but a weaker flux at low latitudes, exactly the inverse to what the model by Cao et al. (2014) requires.

A new model by Takahashi et al. (2019) relies on a completely different mechanism to break the north-south symmetry. A strongly subadiabatic but uniform core-mantle boundary heat flux promotes a thick stratified layer that occupies the outer 50% of the core radius. The inner core radius amounts to 20%. Convection in the in-between dynamo region is driven by a mixture of thermal and compositional effects. Both are modeled with separate evolution equations and the thermal diffusivity is assumed to be an order of magnitude larger than the compositional counterpart. The true difference should be several orders of magnitude larger, but their choice is dictated by numerical limitations. Somewhat surprisingly, the resulting field configuration is very Mercury-like and quite similar to the one reported by Tian et al. (2015) (see Fig. 10), but without imposing a north-south symmetry breaking through the core-mantle boundary heat flux pattern.

Helicity $h=\underline{u}\cdot ( \nabla \times \underline{u})$ is an essential ingredient in the dynamo process. In planetary dynamos, the axial helicity component $h_{z}=u_{z}\cdot (\nabla \times \underline{u})_{z}$ is particularly important, where the index $z$ denotes the direction of the rotation axis. Takahashi et al. (2019) shows in dedicated dynamo simulations that the relative $h_{z}$ is significantly stronger in the northern hemisphere than in the southern hemisphere, which could explain the offset of dynamo action towards the north. A simulation at identical parameters but without magnetic field, on the other hand, did not show any noteworthy hemispherical helicity asymmetry. This strongly suggests that the Lorentz force is responsible for the symmetry breaking.

Sreenivasan and Jones (2011) have shown that Lorentz forces can promote local helicity generation. This would in turn intensify the local magnetic field and thus the Lorentz forces, setting off a runaway process, which could start with smaller natural fluctuations and result in the stronger asymmetry found by Takahashi et al. (2019). In principle, this could also happen in other dynamo simulations but has not been reported so far. A possible reason could be related to the separate modeling of thermal and compositional buoyancy made by Takahashi et al. (2019). The only other Mercury dynamo models following this approach yield much less Mercury-like field geometries, but also assume somewhat different setups and parameters (Manglik et al. 2010). The main source of these discrepancies, however, is still unclear.

5.4.4 Towards More Realistic Dynamo Models for Mercury

The current numerical dynamo models strongly suggest that the outer part of Mercury’s core is stably stratified. This seems like the most likely explanation for the weakness and the strong degree of axisymmetry of the planet’s internal field. The large quadrupole to dipole ratio, often interpreted as an offset dipole, is another striking feature. It was already within the possible solutions at Mariner times and may thus have lasted at least 45 years.

At Earth, the dipole component changes rather slowly over millennia. Higher harmonics vary faster and the characteristic time scale for the quadrupole is about 200 years. Estimating the respective time scales for Mercury is difficult. A stable layer would filter out faster variations so that a direct observation is also problematic. However, time scale significantly faster than centuries seem unlikely. Whether the high quadrupole to dipole ratio represents only a snapshot thus remains unclear. Three possible reasons for a permanently high ratio have been discussed above, but all need to be further explored.

The BepiColombo mission will offer a better global coverage and separation of the different field contributions. It will thus allow constraining more details of the dynamo generated field. Spectra beyond the currently realistic degree $\ell =3$, perhaps 4, will help to more clearly distinguish between different dynamo models, dismissing some of the hypotheses. We may even be able to constrain the size of the inner core or the thickness of the stable layer based on magnetic observations. The tangent cylinder that touches the inner core boundary is an important dynamical boundary in rotating spherical shells. At Earth it leaves its trace in the magnetic field, which is weaker inside than at or just outside the tangent cylinder. However, this effect could only help to constrain the size of Mercury’s inner core if Mercury also harbors an dipole dominated dynamo. According to the line of arguments presented above, this seems unlikely.

The geomagnetic field also allows to deduce the top of the dynamo region. With the exception of the dipole, the magnetic spectrum seem to be roughly white at the core-mantle boundary. Simple upward continuation of the potential field then predict how the tilt of the spectrum changes with radius. The measured tilt can then be used to deduce the depth of the core-mantle boundary, which works amazingly well. At Mercury two problems may arise: 1) an additional stable stratified layer would further influence the tilt and 2) the spectrum determined by BepiColombo may not cover enough spherical harmonics to clearly identify a tilt. A comprehensive geophysical interpretation of magnetic and geodetic constraints on Mercury’s internal structure will be then fundamental to enhance our knowledge of its deep interior.

5.5 Long-Term Evolution of Mercury’s Interior

Geological and spectral measurements of Mercury’s surface are intimately related to the thermal evolution of the deep interior. Similarly, the thickness, composition, and mechanical properties of the crust and lithosphere, which can be inferred (directly or indirectly) by combining various sets of remote observations, ultimately result from the way heat has escaped from the mantle and core over billions of years. Understanding the long-term evolution of the interior is thus a fundamental step towards the interpretation of Mercury’s observational record.

Despite the uncertainties in the formation mechanism that led to Mercury’s peculiar interior structure (e.g., Asphaug and Reufer 2014; Lykawka and Ito 2017; Nittler et al. 2018)), the planet likely underwent a magma ocean phase early in its history. The compositional diversity of the northern hemisphere (Weider et al. 2015) has been interpreted as the result of partial melting of a chemically-heterogeneous mantle originating from the fractional crystallization of a magma ocean (e.g., Charlier et al. 2013). Furthermore, the so-called low-reflectance materials identified in heavily cratered regions (Murchie et al. 2015) are compatible with a carbon-rich composition (Peplowski et al. 2016) thought to originate from a buried layer of graphite (see Sect. 5.1), which formed as a flotation crust in a crystallizing magma ocean (Vander Kaaden and McCubbin 2015).

Following core formation and magma ocean crystallization, secular cooling of mantle and core is largely governed by internal heat production due to the long-lived radiogenic isotopes of K, U, Th, and heat loss at the surface, which is controlled by convection and conduction in the solid mantle. Due to the strong dependence of the mantle viscosity on temperature, the cold outermost layers of a planet are extremely stiff and tend to form an immobile “stagnant lid” (e.g., Solomatov 1995). In contrast to the Earth, where plate tectonics operates and the oceanic crust and lithosphere are an active part of the mantle convection system, Mercury has likely been characterized by a stagnant lid throughout its history.

The evolution of Mercury’s interior has been quantified with models of varying complexity, from one-dimensional, so-called ‘parameterized’ models (Hauck et al. 2004; Grott et al. 2011; Tosi et al. 2013; Knibbe and van Westrenen 2018; Hauck et al. 2018), to fully dynamic two- and three-dimensional models (Roberts and Barnouin 2012; Tosi et al. 2013; Michel et al. 2013; Padovan et al. 2017). The former are based on global energy balances of the mantle and core and supplemented by scaling laws to parameterize convective heat transfer. The latter rely instead on the numerical solution of the full set of conservation equations appropriate for highly viscous media such as the mantle. During the past decade, several models have been tailored to fulfill the different observations obtained by the MESSENGER mission, constraining Mercury’s thermal evolution (Roberts and Barnouin 2012; Michel et al. 2013; Tosi et al. 2013; Padovan et al. 2017; Knibbe and van Westrenen 2018; Hauck et al. 2018). Despite their discrepancies in the adopted method, scientific goals, and input observations, these models have led to a relatively consistent picture of Mercury’s long-term evolution.

The thickness and time of emplacement of the crust represent important pieces of information to guide the development of evolution models consistent with the observational record (Tosi and Padovan 2021). The thickness of the putative graphite flotation crust that formed upon magma ocean crystallization (i.e., the primary crust) is unknown, but can hardly exceed a few km (Vander Kaaden and McCubbin 2015). Long-lived mantle partial melting and subsequent volcanism are thought instead to be responsible for the formation of the bulk of the crust observed today (Denevi et al. 2013). Gravity and topography data of the MESSENGER mission have been used to infer Mercury’s crustal thickness (see Sect. 5.1). At any rate, the retrieved figures are broadly consistent with the results of thermal evolution models accounting for the effects of mantle melting and crust production (Tosi et al. 2013; Padovan et al. 2017; Hauck et al. 2018). Such models additionally predict that crust production was largely completed by $\sim 3.5-3$ Ga, in agreement with the age of the youngest large volcanic units (Byrne et al. 2016).

Since crust and mantle residuum have a larger volume than the primordial mantle, the silicate differentiation associated with mantle melting, melt migration and freezing is also accompanied by a net expansion of the planet (Kirk and Stevenson 1989). In addition, the phase of melting and crust production coincides with a period of mantle heating during which convective cooling is offset by radiogenic heat production. As a consequence, the resulting temperature increase also contributes to a net expansion of the planet, which likely characterized the first few hundred million years of evolution (e.g., Grott et al. 2011; Tosi et al. 2013; Hauck et al. 2018). Yet, this expansion phase is difficult to constrain since Mercury’s surface hardly presents traces of extensional tectonics. If present, these have likely been obliterated by the formation of the oldest volcanic units emplaced around 4 Ga ago (Marchi et al. 2013).

Mercury’s surface is rather dominated by compressive landforms associated with surface faulting – lobate scarps and wrinkle ridges – that have long been attributed to periods of global contraction (Solomon 1977). By mapping these tectonic features, displacement–length scaling properties of faults can then be used to estimate the amount of radial shortening experienced by Mercury (e.g, Watters et al. 2004, 2009; Di Achille et al. 2012; Byrne et al. 2014). Indeed, one of the key observables that has aided the development of Mercury’s evolution models is the amount of radial contraction accumulated by the planet throughout its history (see Sect. 5.3). Models consistent with the latest available observations tend to predict that the mantle is still convecting at present-day (Hauck et al. 2018). Nevertheless, even if mantle convection is active today, it is likely sluggish and its effects hardly recognizable in surface observables such as gravity and topography (Padovan et al. 2015; Tosi et al. 2015).

Upon considering solely the contribution of mantle and core cooling after the expansion phase, thermal evolution models have been constructed that satisfy the constraint of $\sim 7$ km of contraction (Hauck et al. 2018). Yet, the process of inner core freezing, being associated with a decrease in volume, also contributes to planetary contraction (Solomon and Chaiken 1976). To what extent core freezing is responsible for the observed global contraction of the planet strongly depends on the composition of the core. On the one hand, freezing of even a small core with a Fe-S composition could easily add several km of contraction, making it difficult to satisfy the available constraints (e.g., Grott et al. 2011). However, the constraint of 7 km or more of global contraction still needs to be evaluated in the context of models including freezing of a Fe-S core. On the other hand, if Si rather than S was alloyed with Fe in the core (Knibbe and van Westrenen 2018), which is likely because of the reducing conditions at which Mercury formed (e.g., Malavergne et al. 2010), the contribution of core freezing could be much smaller due to the small density difference between solid and liquid Fe-Si alloys (Kuwayama and Hirose 2004). The effects of possibly more realistic core compositions including multiple light elements (e.g., S, Si and C) are still difficult to quantify because of the lack of suitable equations of state for complex, non-binary mixtures. In combination with a robust estimate of the size of the inner core, which will be provided by the BepiColombo mission (Sect. 5.2), new experiments on core mixtures will be crucial to make progress in understanding the crystallization of Mercury’s core and its role on global contraction.

The timing and extent of crust production and planetary expansion and contraction discussed above strongly depend on the amount of radiogenic heat sources present in the mantle, which also affect the possibility that convection is presently ongoing (Michel et al. 2013). In absence of meteorites and samples from Mercury, the bulk concentration of heat-producing elements (HPE) cannot be determined directly and represents thus an important unknown upon modeling of the interior evolution. Nevertheless, surface abundances of K, U and Th have been measured by the MESSENGER gamma-ray spectrometer (Peplowski et al. 2011, 2012). Under the (rather strong) assumption that the observed concentrations are representative of the entire crust, the bulk concentration of HPE can be inferred considering that upon melting, these elements, being incompatible, will be partitioned into the crust according to a certain enrichment factor (Michel et al. 2013; Tosi et al. 2013; Knibbe and van Westrenen 2018; Hauck et al. 2018). Notably, models compatible with the available constraints suggest crustal enrichment factors between $\sim 2$ and 4.5, which ultimately imply bulk concentrations of HPEs quite similar to those of the Earth and Mars (Tosi et al. 2013; Hauck et al. 2018).

A few attempts have also been made to estimate the thickness of the elastic lithosphere (Melosh and McKinnon 1988; Nimmo and Watters 2004; James et al. 2016), and its evolution (Tosi et al. 2015). This is an important additional quantity as it contains clues on the surface heat flux and can thus help to reconstruct the thermal history of the planet. From gravity-to-topography ratios of the northern hemisphere, James et al. (2016) inferred an elastic thickness of $\sim 30$ km around 3.8 Ga, roughly in agreement with the work by Nimmo and Watters (2004). However, this low value is difficult to reconcile with independent estimates based on thermal history models and on the elastic response of the planet to the surface temperature distribution induced by the 3:2 resonance (Tosi et al. 2015). These results predict a significantly thicker lithosphere ($>70$ km at 3.8 Ga), in line with the earlier work by Melosh and McKinnon (1988). The reasons of this discrepancy are not fully clear (see for a thorough discussion the work by Phillips et al. (2018)), but could be related to a too simplistic treatment of the relation between lithosphere temperature and thickness, to uncertainties in the time at which Mercury acquired its resonance, or to the fact that the values inferred by Tosi et al. (2015) could be more representative for the present-day. More work both on modeling, e.g., by accounting for time-dependent viscoelastic deformation of the lithosphere, as well as on dating the surface units used for the gravity-topography analysis will be necessary to make further progress on this topic.

With numerous impact basins with diameters of several hundred km (Fassett 2016), Mercury’s surface shows clear traces of the tumultuous early phases of the solar system. Basin-forming impacts deposit a vast amount of energy into the interior that is dissipated into heat, which raises the question of the influence of these events on the overall evolution of the planet. From a global perspective, even the largest impacts such as the one that formed the Caloris basin have only short-lived effects on the evolution of the mantle and core (Roberts and Barnouin 2012). Such impacts locally modify the planform of mantle convection inducing the formation of hot upwellings beneath the impact site and in turn the production of partial melt, which ultimately erupts forming basin-filling melt sheets (Ghods and Arkani-Hamed 2007; Padovan et al. 2017). The volume and time of emplacement of the melt sheets within Mercury’s large basins such as Caloris or Rembrandt can be estimated from dedicated stratigraphic analyses (Ernst et al. 2015) and by dating of the corresponding surface units (Fassett 2016). Padovan et al. (2017) showed that the possibility to match these observations with numerical models crucially depends on the global thermal state of the mantle at the time of the impact, which indicates that regional-scale datasets related to specific impact sites can provide important additional constraints on the global evolution of Mercury.

The presence of a dynamo-generated magnetic field – present or past – on terrestrial bodies also poses fundamental constraints on the evolution of the interior (Tosi and Padovan 2021). However, interpreting Mercury’s present-day magnetic field and the ∼3.7-billion-years-old crustal magnetization inferred during MESSENGER low-altitude campaign (Johnson et al. 2015; Oliveira et al. 2019) in terms of global thermal evolution models remains challenging, largely because of uncertainties in the composition of Mercury’s core (as discussed above) and in the main mechanism(s) driving the dynamo (Wicht and Heyner 2014, and see Sect. 4.3.1). Evolution models generally predict an early phase during which a ‘thermal dynamo’ driven by rapid cooling of an initially super-heated core is possible (Tosi et al. 2013; Hauck et al. 2018). Dynamo action in later evolutionary stages must rely on compositionally driven convection, release of latent heat, or heat produced by radiogenic elements in the core.

6 Gravitational Theories and Heliophysics Experiment

The radio science instrumentation onboard the MPO spacecraft will also enable a precise determination of Mercury’s ephemeris through extremely accurate spacecraft ranging data (Iess et al. 2009). An enhanced knowledge of Mercury’s trajectory will allow testing gravitational theories, and characterizing the properties of the Sun’s interior. Because of Mercury’s proximity to the Sun, strong gravitational forces induce orbital perturbations that are not predicted by the Newtonian physics. The first measure of Mercury’s perihelion precession helped demonstrating Einstein’s theory of general relativity (GR) (Einstein 2019). The unprecedented accuracies of BepiColombo ranging data will be well-suited to constrain the dynamical evolution of Mercury’s orbital elements leading to a thorough investigation of GR space-time modeling in the weak–field approximation.

Planetary ephemerides in the Solar System are modeled through the formulation of Einstein’s field equations in terms of the parametrized post–Newtonian (PPN) parameters (Will 2018b). The main effect on Mercury’s heliocentric accelerations is given by the PPN parameters $\beta$ and $\gamma$ that quantify the nonlinearity in superposition of gravity and space–time curvature produced by a unit rest mass, respectively. These effects are not included in the Newtonian formulation (i.e., $\beta =\gamma =0$), but are fully accounted for in GR (i.e., $\beta =\gamma =1$). The PPN parameters $\alpha _{1}$ and $\alpha _{2}$ also provide measurable Mercury’s orbital perturbations by defining a gravitationally preferred frame (Milani et al. 2002; Imperi et al. 2018). Spatial isotropy would imply that both $\alpha _{1}$ and $\alpha _{2}$ are equal to 0.

Einstein’s theory of GR is based on the validity of the equivalence principle, which states the equality between gravitational and inertial masses. Laboratory tests with torsion–balance experiments provided precise measurements (i.e., ∼10⁻¹³; Wagner et al. 2012) of the weak equivalence principle (WEP) indicating that objects with different composition and structure fall with the same acceleration in a uniform gravitational field. To account for the self–gravitational energy of the observed body ($\varOmega _{B}$), the strong equivalence principle (SEP) was introduced as an extension of Galileo’s postulate (e.g., Bertotti and Grishchuk 1990)

$$ \frac{m^{G}}{m^{I}}=1+\eta \frac{\varOmega _{B}}{m^{I}c^{2}}, $$

(31)

where $m^{G}$ and $m^{I}$ are the gravitational and inertial masses, respectively, $c$ is the speed of light, and $\eta$ is the Nordtvedt parameter that is equal to zero if the SEP is valid. The accurate estimate of this effect requires massive objects, and, therefore, planetary and satellites ephemerides are adequate to measure possible SEP violations. Lunar laser ranging (LLR) over the past 40 years allowed to retrieve several solutions that are consistent with the SEP (Williams et al. 1976, 2004, 2012; Müller et al. 2014). The latest LLR solution led to an uncertainty of $\sigma _{\eta }\sim$ 3.0 ×10⁻⁴ (Müller et al. 2014), and further investigations are necessary to confirm the SEP with enhanced accuracies. A refined knowledge of Mercury’s ephemeris is a unique opportunity to estimate the Nordtvedt parameter in the weak–field approximation, since SEP violations would induce orbital perturbations as well as a redefinition of the Solar System Barycenter (SSB) (Milani et al. 2002).

Mercury’s orbit is also significantly perturbed by the Sun’s interior and dynamics. A better understanding of the Sun’s internal mass distribution is obtained by measuring the solar gravitational parameter ($GM_{ \odot }$) and oblateness ($J_{2_{\odot }}$). A long–duration observation of planetary ephemerides is fundamental for retrieving a precise estimate of the Solar System expansion due to the time variation of $GM_{\odot }$. Combined secular variations in the gravitational constant $G$ and the Sun’s mass $M_{\odot }$ would induce a few m changes in the planet’s relative distance to the Sun in few years (Smith et al. 2018). These effects are detectable through accurate ranging data over an extended period of time.

These helio– and fundamental physics objectives will be addressed by the BepiColombo radio science investigation. MORE highly accurate ranging data will strongly constrain Mercury’s ephemeris leading to unprecedented measurements of quantities related to theories of gravitation and Sun’s interior. Numerical simulations of MORE relativity experiment have been carried out to predict the achievable formal uncertainties of the adjusted PPN and heliophysics parameters (Milani et al. 2002; Imperi et al. 2018). These results were initially yielded through a method that is based on the recovery of Mercury and Earth initial states only (Milani et al. 2002). This estimation strategy relies on a priori assumptions on Earth’s initial conditions (i.e., rotation) and Sun’s position with respect to the SSB (i.e., rescaling) to solve the degeneracy in the inversion of the ephemeris solution (Milani et al. 2002; De Marchi et al. 2016; Schettino et al. 2018). Recent results of Einstein’s GR and solar system expansion through the analysis of the MESSENGER data indicated that discrepancies in the planetary ephemerides affect significantly the estimation of relativistic and heliophysics parameters (Genova et al. 2018). Therefore, BepiColombo data will be processed in a global solution of the Solar System planetary ephemerides with complementary datasets of spacecraft in orbit about other celestial bodies (Fienga et al. 2015). This approach allows improving the planetary ephemerides without the need of any a priori assumption. The latest numerical simulations of the BepiColombo relativity experiment demonstrated the benefits of this method, leading to the results reported in Table 4 (De Marchi and Cascioli 2020).

Table 4 GR and heliophysics results of the MESSENGER solution (Genova et al. 2018) and the BepiColombo latest numerical simulations (Imperi et al. 2018; De Marchi and Cascioli 2020). MESSENGER estimation reports the formal uncertainties (1–standard–deviation, $1-\sigma$) and errors due to discrepancies in the planetary ephemerides suggesting the limitations of the single planet integration only. BepiColombo expected results are in the range of accuracies presented by Imperi et al. (2018) and De Marchi and Cascioli (2020). The former uncertainties were obtained by assuming rotation and rescaling constraints, whereas De Marchi and Cascioli (2020) recovered a global estimation of the planetary ephemerides including other mission datasets (e.g., Mars Reconnaissance Orbiter) with conservative a priori information

The accurate estimation of the PPN parameter $\gamma$ with the BepiColombo relativity experiment stands out. The triple–link radio tracking system of the MORE investigation has been designed to obtain precise data in proximity of superior solar conjunctions leading to an accurate measurement of the radio signal bending and delay, which are proportional to $\gamma +1$ due to the curvature of space–time induced by the Sun (Shapiro 1964). The MORE estimate of $\gamma$ will be an order of magnitude more precise than the current value, which was retrieved by the radio science experiment of the Cassini mission (Bertotti et al. 2003). This PPN parameter also impacts on the planetary dynamical equations leading to perturbations that are highly correlated to $\beta$, $\alpha _{1}$, $\alpha _{2}$, and $J_{2_{\odot }}$ effects. To estimate this full set of quantities, the following a priori constraint was applied

$$ \eta =4\left (\beta -1\right )-\left (\gamma -1\right )-\alpha _{1}- \frac{2}{3}\alpha _{2} , $$

(32)

which is known as Nordtvedt equation (Nordtvedt 1970). This assumption helps especially to disentangle $\beta$ and $J_{2_{\odot }}$ contributions, yielding the accuracies shown in Table 4.

The combination of MORE and other radio science datasets will enable the estimation of $\dot{GM}_{\odot }/GM_{\odot }$. The time–varying solar gravitational parameter will be measured with an unprecedented accuracy, which is fundamental for monitoring the solar mass loss rate and testing possible time variations of the gravitational constant G. The comparison between the measured $\dot{GM}_{\odot }/GM_{\odot }$ and the theoretical computations of the Sun’s mass loss rate (Pinto et al. 2011; Pitjeva and Pitjev 2012) will significantly constrain the universal constant time variation.

Alternative theories of gravitation will also be tested by the MORE investigation. The joint analysis of Mars orbiters range data, which are assumed to have an accuracy of ∼1 m (as, e.g., the Mars Reconnaissance Orbiter; Zurek and Smrekar 2007; Zuber et al. 2007) and to entirely cover the 2025–2026 timeframe, and the BepiColombo range data, for example, will enable the refined estimate of the lower bound on the Compton wavelength $\lambda _{g}=1.1\times 10^{14}$. This result will significantly constrain the physical properties of the hypothetical massive carrier of the gravitational interaction (i.e., massive graviton) (Will 2018a; De Marchi and Cascioli 2020).

7 Conclusions

Improved measurements of the rotation, tides and the magnetic field of Mercury by the BepiColombo mission will enable independent constraints on Mercury’s core. However, the main progress is expected from a combined analysis of those datasets. Only interior models that satisfy all constraints can be a candidate for describing Mercury. Accurate rotational data as will be acquired by the MPO science instruments are required to better constrain the density profiles in the core and silicate envelope. This places essential bounds on the composition of the core and silicate envelope in addition to surface compositional data. The bulk composition of the core is intimately linked to the size and composition of the inner core. An indirect knowledge of the composition of the core can, therefore, be obtained by investigating the inner core, whose presence and properties will be studied by measuring the libration amplitude at 88 days, the orientation of the planet (i.e., pole coordinates $\alpha$ and $\delta$), and the tidal Love numbers $k_{2}$ and $h_{2}$.

Further constraints will be provided by the magnetic field that is intimately linked to the thermal state of the core and solidification processes. An accurate modeling of the internal dynamo will allow constraining the size of the inner core and the properties of a possible stably stratified layer below the core-mantle boundary. Current dynamo models for Mercury, however, are still limited by the available computational power, and different assumptions are required. Existing dynamo models are based on the assumptions of large viscosities to damp the small scale flows or of simplistic representation of the mixed convective driving by thermal and compositional buoyancy differences. A significant improvement in performances of the available resources will allow implementing more realistic models in the future, leading to new insights into Mercury’s deep interior.

Additional information on the inner core will be obtained by the observation of long-period librations, if the free librations are close to long-period librations or happen to be excited to an observable level. Since the size of the inner core depends strongly on the melting temperature of the core alloy, which itself depends on the core composition, an accurate characterization of iron alloy properties is required as a fundamental complementary observation to BepiColombo geodesy and geophysics data. Both laboratory studies and theoretical investigations can contribute to improving the amount and quality of data on the density, liquidus and solidus temperature, elastic properties, and transport properties of core alloys. Besides binary alloys with one light element, such studies should also consider ternary compositions as Fe–Si–S. Furthermore, other light elements as carbon (Shimoyama et al. 2016; Steenstra et al. 2017), oxygen (Tsuno et al. 2007), nickel (Ringwood 1959; Zhu et al. 2014), and hydrogen (Clesi et al. 2018) might also be present in small concentrations in Mercury’s core, and their possible effects on Mercury’s core need to be assessed.

Significant progresses in understanding Mercury’s interior and evolution have been made during the MESSENGER mission. Several questions are still open, and the combination of the BepiColombo observations with a new generation of models and laboratory data will provide crucial information on Mercury’s internal structure and history. Furthermore, the accurate determination of Mercury’s ephemeris through the analysis of the MORE data will allow addressing fundamental objectives regarding theories of gravitation, and Sun’s interior and evolution.