Up-to-fivefold reverberating waves through the Earth's center and distinctly anisotropic innermost inner core

The study addresses the long-standing challenge of probing the Earth’s centermost region to understand planetary formation, evolution, and the geodynamo. Traditional seismological probes (e.g., PKIKP travel times and normal modes) provide limited sampling of the innermost inner core (IMIC) because suitable earthquake–station geometries near antipodes are rare and normal modes have low sensitivity at the center. The research question is whether previously unobserved, direct-wave reverberations of compressional waves passing multiple times along Earth’s diameter (PKIKP multiples) exist in seismograms and, if so, whether their differential travel times can constrain anisotropy in the inner core—particularly the properties of a hypothesized IMIC distinct from the outer inner core (OIC). The purpose is to leverage modern, dense global networks and waveform stacking to reveal these weak signals, thereby improving sensitivity to the centermost ~650 km of the inner core and refining models of its anisotropy and structure.

Early work after discovery of the inner core focused on isotropic structure and the inner-core boundary. Since the 1980s, multiple studies have revealed inner-core anisotropy, initially proposing depth-independent cylindrical anisotropy inferred from PKIKP travel times and normal-mode splitting. Subsequent work identified hemispherical dichotomy and radial variations in anisotropy, with more recent models invoking complex spatial variations in P-wave anisotropy and attenuation. S-wave anisotropy has also been reported. Mineral physics debates whether hcp or bcc iron stabilizes under inner-core conditions; both can produce anisotropy patterns consistent with seismological observations depending on preferred crystal orientations. The IMIC was hypothesized as a central ball with anisotropic properties distinct from the OIC, with early estimates suggesting a ~300-km radius and a slow direction oblique to the Earth’s rotation axis (ERA). Subsequent analyses using ISC datasets, dedicated PKIKP picks at antipodal distances, and normal modes provided supporting evidence, but uncertainties persisted about the IMIC’s radius, transition to the OIC, and precise anisotropic strength and directions. Coda-correlation approaches emerged as a complementary tool with different sensitivity kernels to the deep Earth and were recently used to constrain inner-core anisotropy, though their complex kernels remain challenging to interpret. These studies motivated a search for exotic, direct-wave PKIKP multiples to more directly sample the centermost inner core.

Data retrieval and preprocessing: Broadband seismograms were obtained from IRIS, ORFEUS, GFZ (GEOFON), ETHZ, and INGV. Stations with continuous 120-minute recordings from event origin time were used. Instrument responses were removed to obtain velocity seismograms, resampled to 10 sps (N=72,000 samples per trace). Processing used obspyDMT and ObsPy. Global stack construction: For each large earthquake (Mw ≥ 6.0) over the last decade, records were bandpass filtered 10–100 s (zero-phase, three corners) and grouped into 1° epicentral distance bins. A median-amplitude outlier rejection discarded traces with peak amplitudes >5× bin median. Remaining traces were linearly stacked per bin to build 2D lapse-time vs distance images (global stacks). To enhance visibility of very late arrivals (e.g., PKIKP4, PKIKP5), stacks were multiplied by a common polynomial of elapsed time f(t)=(t×10^1)^4; visualization used linear interpolation. Identification of PKIKP multiples: Flat, near-zero slowness features consistent with steep-incidence phases reverberating along Earth’s diameter were identified and labeled PKIKP (I), PKIKP2 (I2), PKIKP3 (I3), PKIKP4 (I4), and PKIKP5 (I5), indicating 1–5 passages through the diameter. Multiple examples were found in single-event stacks; three- and fourfold multiples are commonly observed; fivefold multiples appear clearly in a few events (e.g., Mw 7.9 Solomon Islands, 22 Jan 2017). Regional array stacking and measurement: For differential time measurements, regional dense arrays (USArray continental and Alaska branches; AlpArray; European networks) were used. For these arrays, a shorter period band 7–13 s (sometimes 1–10 s for specific cases) was applied to improve SNR given more coherent mantle paths. Individual traces were time-corrected for ellipticity and mantle heterogeneity using ak135 predictions plus 3D mantle models (DETOX-P3 as primary; MIT-P08 and LLNL-G3Dv3 for sensitivity tests). Traces were then stacked across the array to enhance weak higher-order multiples. Differential residuals: For podal geometry (epicentral distance <50°), residuals between PKIKP4 and PKIKP2 were measured; for antipodal geometry (>155°), residuals between PKIKP3 and PKIKP were measured. Differential residual Δt = (PKIKP(n+2)−PKIKP(n))obs − (PKIKP(n+2)−PKIKP(n))pred were obtained via cross-correlation of stacked pairs. Bootstrap with 5000 resamples of stations quantified uncertainties. Residuals were converted to fractional velocity perturbations Δv/v using theoretical inner-core segment travel-time differences for the relevant PKIKP legs. Sampling direction parameterization: A representative sampling angle ξ′ relative to the ERA was defined for each pair: for PKIKP4−PKIKP2, ξ′ combines forward and backward PKIKP2 legs; for PKIKP3−PKIKP, ξ′ equals the PKIKP angle. Distance constraints ensured limited deviation of ξ′ from individual leg angles (<25° for PKIKP/PKIKP2; <12.5° for PKIKP3/PKIKP4). Mean and standard deviation of cos²ξ′ across array elements quantified direction and its uncertainty. Anisotropy modeling: Cylindrical transverse isotropy was assumed with Δv/v = γ + ε cos²ξ + α sin²ξ cos²ξ (Creager, 1992). Orthogonal Distance Regression (ODR) estimated anisotropy parameters accounting for uncertainties in both cos²ξ′ and Δv/v. Monte Carlo simulations propagated correlated parameter uncertainties. Two model configurations were tested: (1) bulk IC model fitted directly to data; (2) two-layer IC with a fixed outer shell (OIC) from Stephenson et al. (2021) and an IMIC of radius H=650 km inverted for distinct anisotropy. Sensitivity tests with alternate mantle models (MIT-P08, LLNL-G3Dv3) and no mantle correction assessed robustness. Finite-frequency 2D spectral element experiments illustrated feasibility and sensitivity of 7–13 s PKIKP multiples to a 650-km-radius IMIC.

• Discovery of up-to-fivefold PKIKP reverberations: Global single-event stacks reveal PKIKP, PKIKP2, PKIKP3, PKIKP4, and in a few cases PKIKP5 arrivals—multiples beyond two passages were previously unreported in direct seismograms.
• Robust measurement set: Using regional arrays and differential timing of PKIKP4−PKIKP2 (podal) and PKIKP3−PKIKP (antipodal), 16 differential travel-time measurements were obtained across events in Alaska, the mainland U.S., and Europe.
• Strong IMIC anisotropy with slow direction at mid-latitudes: Inferred cylindrically anisotropic models indicate the slowest P-wave direction at ~48–50° from the Earth’s rotation axis in the centermost inner core.
• Quantified anisotropy strengths: • Bulk IC fit (DETOX-P3 corrected): P-wave speeds are ~2.8% faster along the polar direction and ~1.7% faster along the equatorial plane relative to the slowest direction (~48°), implying a pronounced anisotropic pattern in the center-driven model. • Two-layer IC with IMIC (H=650 km): Within the IMIC, P-wave speeds are ~4.0% faster along the polar direction and ~3.4% faster along the equatorial plane relative to the slowest direction (~48°). The OIC exhibits weaker anisotropy with slow direction near the equatorial plane.
• Parameter estimates indicate slow directions deviating from the equatorial plane (e.g., σ≈8.3 for bulk model; σ≈−15.6 for IMIC case, both implying non-equatorial slow directions).
• Robustness: Similar anisotropy patterns (slow direction at oblique angles ~48–50°) are recovered using different mantle corrections (DETOX-P3, MIT-P08, LLNL-G3Dv3) and even without mantle correction, supporting the stability of the inferred depth-dependent anisotropy.
• Practical advantages: Differential timing reduces source location and path heterogeneity effects; repeated passages sample the centermost ~650 km multiple times, amplifying sensitivity to the IMIC.
• Geophysical implications: Results strengthen evidence for a distinctly anisotropic IMIC transitioning to a weakly anisotropic OIC; consistent with independent PKIKP, normal-mode, and coda-correlation studies; compatible with both hcp and bcc iron textures under inner-core conditions.

The findings directly address the central question: whether reverberating PKIKP multiples exist and can be exploited to probe inner-core anisotropy. By identifying and measuring up to fivefold direct-wave reverberations, the study demonstrates a new, practical class of observations with clear sensitivity kernels along PKIKP ray paths. Differential travel-time analysis of carefully selected podal and antipodal pairs mitigates earthquake location biases and reduces mantle heterogeneity impacts, enabling precise constraints on anisotropy in the centermost inner core. The inferred anisotropy, with the slowest direction at ~48–50° to the ERA and stronger anisotropy in the IMIC than in the OIC, is consistent with and extends independent lines of evidence (absolute PKIKP, dedicated antipodal picks, normal modes, and coda-correlation). This convergence strengthens the case for a distinct IMIC and suggests a change in inner-core growth regime or deformation history, potentially linked to thermal/compositional evolution and geodynamo behavior. Practically, the approach opens a new north–south sampling style and complements limited existing geometries (e.g., South Sandwich–Alaska paths) by accessing the very center with podal configurations. The results are relevant to mineral physics debates (hcp vs bcc iron textures) and geodynamic models explaining depth-varying anisotropy, including evolving convection strength, preferential crystallization, and chemical stratification at the ICB.

By leveraging modern global and regional seismic networks, the study reveals and exploits previously unreported direct-wave PKIKP multiples (up to fivefold) to probe the Earth’s centermost inner core. Differential timing of podal and antipodal multiple pairs yields robust constraints on cylindrical P-wave anisotropy: a distinctly anisotropic IMIC (~650-km radius) with the slowest direction at ~48–50° to the rotation axis, transitioning to a weakly anisotropic OIC with slow directions near the equator. These findings consolidate multiple independent seismological lines of evidence for the IMIC and imply a significant shift in inner-core growth and deformation history. Future research should focus on characterizing the IMIC–OIC transition depth and nature, expanding higher-frequency observations with smaller-aperture networks, improving mantle heterogeneity corrections, and integrating mineral physics and geodynamic models to link observed anisotropy patterns to inner-core composition and dynamics.

• Data and geometry constraints: Only a subset of large earthquakes with favorable source mechanisms and geometries produce clear higher-order PKIKP multiples; fivefold reverberations are rare. Access to some large arrays (e.g., ChinArray) was limited, and introducing new global sampling directions will be challenging without future large-scale deployments.
• Frequency and SNR limitations: Long-period filtering (10–100 s globally; 7–13 s regionally) is needed to mitigate heterogeneity-induced misalignments; higher-order multiples have low SNR and require stacking, limiting event-wise resolution.
• Modeling assumptions: Cylindrical transverse isotropy is assumed; anisotropy is inferred via ODR fits with uncertainties. Near-center measurements have minimal depth resolution, so data fits alone cannot decisively distinguish between bulk-IC and two-layer IC models without assuming an OIC model.
• Heterogeneity corrections: Results depend on ellipticity and mantle corrections; although robustness was tested across multiple models (DETOX-P3, MIT-P08, LLNL-G3Dv3), residual mantle effects remain a source of uncertainty.
• Kernel simplifications: While direct-wave sensitivity is well localized along ray paths, finite-frequency effects and 3D structure could introduce complexities not fully captured by the adopted parameterization.