Seismic detection of a 600-km solid inner core in Mars

Data preprocess

To avoid pitfalls related to transient one-sided pulses, typically referred to as glitches⁵⁰, we first use two glitch-removal algorithms sequentially: the MPS, followed by the UCLA method proposed in ref. ⁵¹. Next, the UVW components are rotated to the ZNE components⁵². A Butterworth bandpass filter is then used to filter the glitch-reduced data to 0.2–0.8 Hz to avoid the 1-Hz tick noise, high-frequency idiosyncratic signals and long-period noises. Finally, we apply a time-domain polarization filter method⁵³, routinely used in planetary seismology^{7,10,20,54,55,56}, to suppress the non-linearly polarized noise and enhance the linearly polarized body wave phases.

Array analysis

Vespagram analysis^30,57 can efficiently distinguish different seismic phases based on their unique travel time and slowness. By stacking the traces shifted along a specified slowness, a given phase is enhanced. Vespagram analysis has been applied to the detection of core-reflected phase ScS on Mars⁷. Unlike typical receiver array analyses for earthquakes, which allow for direct stacking waveforms to preserve important information about the amplitude and polarity, the vespagram analysis here is based on a source array³⁰. Variations in marsquake properties, such as focal mechanism, magnitude and depth, combined with complicated propagation effects from possible three-dimensional structures, cause large differences in absolute amplitudes and waveforms among different events, making vespagram analysis more challenging.

To minimize absolute amplitude differences among different marsquakes, we normalize each polarization-filtered trace within the targeted time window as

where μ and σ are the mean value and standard deviation of waveform data X, respectively. In the subsequent section on determining the amplitude and polarity of PKKP and PKiKP, we also normalize the entire waveform record using the amplitude of the reference phase or by deconvolving the reference phase.

Both polarization-filtered waveforms from the vertical component and their envelopes (Supplementary Information section A3.2.2.1) are tested as inputs for vespagram analysis. To improve the slowness resolution, non-linear stacking methods are applied. For waveform data, we test both fourth-root vespagram and phase-weighted stack (PWS)³⁰, whereas envelope data use only the fourth-root vespagram. It is important to note that differences in waveform polarities and source durations among events may reduce the stacking energy even when they are aligned with a specified slowness. Vespagrams for synthetic P waveforms, assuming different source mechanisms and durations for different events, demonstrate that envelope data produces stronger and more focused energy than waveforms. Moreover, owing to the variation and incoherence in waveforms across different events, PWS cannot produce focused energy (Supplementary Fig. 8).

Moreover, FDPA (see section ‘Analysis for individual events’ for details) shows distinct differences in vertical (VRM) and horizontal (HRM) rectilinear motion between compressional and shear waves. To further enhance vertically polarized signals, we test multiplying the VRM–HRM to the polarization-filtered waveforms and conduct vespagram analysis subsequently (Supplementary Fig. 24).

Here we select 23 low-frequency marsquakes (Extended Data Table 1) and focus on using the fourth-root vespagram to detect core phases (Supplementary Information section A3.2.2). The traces are aligned with the direct P provided in the MQS²⁹. Therefore, the travel time and slowness derived from the vespagram analysis are relative to those of the direct P. The reference epicentral distance is set at 29°. The grid exhibiting the highest coherent energy within the targeted time-slowness window, and demonstrating the most consistent presence in the bootstrap tests, is then picked as our candidate phase.

Assessment of vespagram result robustness

To address uncertainties and assess the reliability of our vespagram results, bootstrap resampling test⁵⁸ is used. We conduct three types of bootstrap resampling tests.

In type I, we randomly select two-thirds of the events to generate a vespagram, minimizing the influence of events with anomalous amplitudes.

In type II, given that most events are concentrated at a distance range of 29°–32° (Extended Data Table 1 and Supplementary Fig. 7), which could bias the vespagram results, we randomly select half of the events from this distance range, along with events at other epicentral distances, to generate the vespagram.

In type III, we randomly shift each trace within a range of −10 s and 10 s before stacking to mitigate possible uncertainties in the epicentral distance and picked P-arrival time.

For all tests, we repeat the random resampling process 200 times and calculate the mean vespagram along with the corresponding 95% confidence levels.

Contamination from glitches could potentially introduce artefacts⁵⁰, even with the application of glitch-removal processes. To evaluate possible effects of glitches, we conduct additional vespagram analysis on the raw data, and the data contain only identified glitch signals (Supplementary Fig. 23). The results indicate that the energy of the core phases identified on the vespagram is not attributable to these glitches.

Analysis for individual events

Guided by the results from the vespagram analysis, we apply complementary approaches to identify core phases for individual events, including time-domain envelopes (filter bank)⁵⁶ and polarized analysis³³.

In the filter bank approach, we first apply a zero-phase-shift bandpass filter to the vertical velocity trace of each event, using a bandwidth half an octave wide around each central frequency, ranging from 1/16 Hz to 2 Hz. Next, we use a time-domain polarization filter to construct a filter bank and compute envelopes in 5-s-long time windows. Given the uncertainty in marsquake locations, to better identify the arrivals, we select the time-domain envelope peaks within a ±15 s time window around the predictions from the vespagram analysis, focusing on a frequency band from 1/4.8 Hz to 1/1.2 Hz, for which clear energy packages are consistently observed.

We further use the FDPA method^7,20,59. We start by computing the S-transform of three-component waveforms and a 3 × 3 cross-spectral covariance matrix using 90% overlapping time windows, with duration varying inversely with frequency. The relative sizes of the eigenvalues of this covariance matrix indicate the degree of polarization (DOP) of the particle motion, ranging from 0 and 1, whereas the complex-valued components of the eigenvectors describe the particle motion ellipsoid. A polarized wave results in a notably larger eigenvalue and a DOP close to 1. The orientation of the semi-major axis provides information on the inclination, with the back-azimuth derived from the projection of the axis onto the horizontal plane and the deviation from the vertical plane determining the inclination. We then identify seismic arrivals with rectilinear particle motion dominantly polarized in vertical (VRM) and horizontal directions (HRM). Given the nearly vertical incidence of the core phases, a strong peak in the VRM–HRM plot would confirm their presence (Extended Data Fig. 2).

Travel time picks and measurement uncertainties

For the vespagram results, we estimate the uncertainties in the slowness and travel time of the candidate phases using the bootstrap resampling test type I (refs. ^60,61). In each vespagram from the resampling tests, grids with energy exceeding 85% of the peak are assigned a value of 1, whereas others receive a value of 0. We then sum the values at each grid and calculate their percentage occurrence among all tests (Fig. 2b,e), showing a concentrated distribution at the slowness and travel time of the candidate core phases. To evaluate the effect of the threshold selection, we also test two alternative levels: 50% and 70% (Supplementary Fig. 16). Although lower thresholds yield higher percentage occurrences, they also produce more dispersed energy distributions. Accordingly, we adopt the 85% threshold as an optimal balance between robustness and precision. Focusing on a window around these phases, we fit the distribution separately along the travel time and slowness axes with Gaussian functions (Fig. 2e) to derive mean values for each. The uncertainties are represented by the 1σ values obtained from these analyses.

In the analysis for individual events, the RMS value of envelope peaks from 1/4.8 Hz to 1/1.2 Hz in the filter bank is used to determine the arrival time of our candidate phase (Extended Data Fig. 2a). The corresponding standard deviation is computed to estimate the travel time uncertainty for subsequent inversions.

Amplitude and polarity of PKKP and PKiKP

To better capture the relative amplitude between different phases, we first normalize the time window for a reference phase (for example, P or PKiKP) using equation (1). Other phases are then normalized using the same ratio as the reference phase. Both stackings produce similar results, which shows that a stacked amplitude ratio of around 0.5 between PKKP and PKiKP (A_PKKP/A_PKiKP) (Supplementary Fig. 28).

For the polarity, we first adjust the traces to ensure all P waveforms have positive polarity. However, the polarity of PKiKP can vary, either positive or negative, depending on the focal mechanism. Therefore, this approach does not yield focused energy. We then test by setting the polarity of all targeted PKiKP arrivals to positive, and analyse the corresponding vespagram for PKKP. Nevertheless, this adjustment does not markedly improve the results because of possible source complexity and complicated propagation effects (Supplementary Fig. 29).

Source complexity and propagation effects at shallow depths can be partially removed by deconvolving the P (using −2 s to 8 s of the P-arrival time here) or PKiKP (−5 s to 5 s) waveforms. Vespagrams from deconvolved waveforms show more focused energy than those from the original waveforms (Supplementary Fig. 30). Notably, deconvolving PKiKP shows that the vespagram of PKKP has opposite polarity to PKiKP, with an A_PKKP/A_PKiKP of about 0.5, consistent with the individual event measurement (Extended Data Fig. 3 and Supplementary Fig. 31). However, we exercise caution about whether the selected P or PKiKP waveforms accurately represent the source. Consequently, extracting amplitude and polarity information for marsquakes with low signal-to-noise ratios and unknown source mechanisms remains challenging. Moreover, using envelope data proves to be more effective and stable for vespagram analysis in detecting core phases.

We also measure the A_PKKP/A_PKiKP for event S0235b, which has concurrent and clear observations of both phases. In the bandpass-filtered data (Extended Data Fig. 7a), the amplitudes of PKiKP and PKKP are determined as the maximum amplitude in the target windows. We estimate the amplitude uncertainties as the maximum amplitude in the noise window, defined as 5–15 s before the main phases. This allows us to calculate the A_PKKP/A_PKiKP, along with its uncertainty, as shown in the blue-shaded region in Extended Data Fig. 7c.

Seismic inversion

To solve the inverse problem of determining core seismic velocity profile, we use the probabilistic approach⁶², with a solution given by

where k is a normalization constant, f(m) is the prior model parameter probability distribution, L(m) is the likelihood function, which is a measure of the similarity between the data and the predictions from model m, and σ(m) is the posterior probability distribution.

Assuming that data noise is uncorrelated and described by a Laplace distribution (L₁-norm), the likelihood function takes the form

where φ is the misfit function. The general expression for the misfit is

where d_obs and d_cal denote the vectors of observed and synthetic travel times relative to the P, with N expressing the total number of travel times, which comprises all possible combinations of phases with respect to P: ΔT_PKiKP–T_P, ΔT_PKKP–T_P, ΔT_{P′P′r_ab}–T_P, ΔT_P′P′n–T_P and σ is the travel time uncertainty obtained from either vespagram or individual event analysis.

We use two strategies to define the travel time dataset and the associated σ. The first strategy, referred to as M_vesp, uses travel time picks of core phases and their uncertainties from the vespagram as d_obs and σ, respectively (Supplementary Table 1 and Supplementary Fig. 2). In this case, each core phase has only one travel time pick and the corresponding d_cal represents the theoretical travel time at a single epicentral distance of 29°. Alternatively, the second strategy, M_pick, uses travel times relative to P for selected individual events at their epicentre distances as d_obs (Supplementary Tables 3–6).

As the slowness of PKiKP offers additional constraints on IC size, we also incorporate the slowness into the M_vesp by modifying the misfit function φ as

where p_{PKiKP_obs} and σ_p represent the observed slowness and its corresponding uncertainty derived from the vespagram analysis in Fig. 2e, whereas p_{PKiKP_cal} is the calculated slowness for a select model; w denotes the weight balancing the contributions of travel time and slowness in the misfit function. By incorporating the slowness of PKiKP (Extended Data Table 2, case 9), the inverted IC radius is consistent with those obtained by inverting only travel times, but with a more concentrated posterior distribution.

Finally, to sample the posterior distribution (equation (2)), we use the Metropolis sampling algorithm⁶². This algorithm randomly samples the model space, ensuring more frequent sampling of models that align with previous information while simultaneously fitting the data.

Seismic model parameterization

Given that our target seismic phases are primarily related to the Martian core, a crucial consideration is whether to invert the Martian mantle structure simultaneously. As PKiKP and PKKP travel with nearly identical paths in the mantle (Fig. 1b), their differential travel time are primarily affected by the core structure, with only minimal impact from the one-dimensional mantle structure variations (see Supplementary Information section A4.1 for more details). This supports our decision to focus solely on inverting core velocity structure with ΔT_PKKP–P and ΔT_PKiKP–P. By opting for differential travel times over absolute travel times, we also mitigate the effects of the uncertainties of the mantle model. Furthermore, the mantle near the antipode sampled by the P′P′ could potentially be heterogeneous^{59,63,64,65,66}. Therefore, we refrain from amalgamating travel times obtained in previous studies^{7,20,31,32,56,59} and our new measurements to invert the entire Martian velocity structure. Instead, we invert the P-wave velocity (V_P) structure of the core based on different mantle models and assess how different models can affect the inverted core models. Seven representative mantle models: the mean, lower and upper bound V_P models from the AK_subset models⁷; the two inverted V_P models incorporating SKS travel times²⁰; and the two models with MSL^31,32 are selected to perform inversion separately. These mantle models encompass a range of models produced in previous studies.

The Martian core P-wave velocity structure is then parameterized with five model parameters (Extended Data Fig. 4): the OC radius (R_OC), the V_P at the CMB (V_P__CMB), the V_P of the OC side at the ICB (V_P__OC__ICB), the IC radius (R_IC), and the V_P jump at ICB (δV_P__ICB). Here, the model spaces of the OC are based on the results in ref. ²⁰. In the case of the model with an MSL overlying the liquid OC (Extended Data Fig. 4c), the thickness of MSL, the mean velocity gradient in MSL and the V_P jump at the bottom of MSL (the real CMB) are set according to the two models, MSL_ETH³¹ and MSL_IPGP³². To simplify the travel time calculation using TauP⁶⁷, we treat the MSL as a part of the liquid OC. Thus, our identified PKKP corresponds to PKKP_MSL in this case and the inverted R_OC contains the thickness of MSL.

We further assume that the mean velocity gradients in the OC and IC are the same. The core velocity at a certain depth r is interpreted with a smooth function as

Here, r₁ and r₂ are the depths of the CMB and ICB, respectively. \({V}_{{{\rm{P}}}_{1}}\) and \({V}_{{{\rm{P}}}_{2}}\) denote the corresponding velocities at these interfaces.

The inverted core velocity profiles, based on different mantle models and strategies, are listed in Extended Data Table 2.

Composition of the Martian core

The presence of H substantially lowers the melting point of Fe, making it difficult to form a solid IC⁴¹. Therefore, in our density and velocity models, we considered only O, C and S as light elements in the Martian core. We examined two distinct core composition models. In the first model, when the C content exceeds 4 wt%, a C-rich IC may precipitate^42,43. Under this scenario, the O and S contents range from 0 to 6 wt% and 12 to 16 wt%, respectively²¹. The second model assumes an O-enriched core. A solid O-rich IC will precipitate when the oxygen content exceeds about 6.7 wt% (ref. ²¹). Based on the seismologically inferred radius of the IC, we estimate the O content in the Martian core to be 6.7–9.0 wt%, corresponding to C and S concentrations of approximately 0–4 wt% and 12–16 wt%, respectively²¹. Both models account for the partitioning of O and C between the Martian OC and IC. The density and velocity profiles of solid Fe₃C and FeO under high pressure and temperature conditions are as follows^68,69:

The influence of O, C and S on the OC velocity and density are calculated as follows⁵⁴:

where Cont._O, Cont._C and Cont._S are mole per cent of O, C and S used in our modelling.