Hidden pressure-stabilized lead reservoirs in Earth’s mantle

Domain: physics
Published: 2026-02-18
Editorial outcome: minor_revision

Abstract
The isotopic imbalance between the accessible silicate Earth and chondritic meteorites implies a missing reservoir of unradiogenic lead. While prior models invoked sequestration in the core, the mechanisms capable of isolating lead from its radiogenic parents within the silicate Earth remain uncertain. Here we show, using first-principles structural prediction and melting-point calculations, that high-pressure Pb-S phases could stably host unradiogenic Pb deep within the mantle. PbS crystallises early and remains solid to core-mantle boundary conditions, potentially forming a long-lived reservoir isolated from U and Th. Under sulfur-rich conditions, PbS may react to form polysulfide phases, PbS2 and PbS3, whose contrasting melting behaviours could promote both storage and limited upward transport of Pb. These findings provide a plausible physical mechanism for the preservation and episodic release of ancient lead, which could account for unradiogenic signatures in the upper mantle. The results suggest a link between mantle redox and sulfur cycling to the isotopic evolution of the silicate Earth.

Introduction
Lead isotopes systematics underpin our understanding of Earth’s differentiation, yet the accessible mantle and crust are markedly more radiogenic than chondritic meteorites—the so-called “missing lead paradox”. This disparity implies that a portion of Earth’s original, unradiogenic lead is sequestered in hidden reservoirs. Although partitioning into the metallic core has been proposed, the known partition coefficients of Pb between metal, sulfide, and silicate phases are insufficient to account for the deficit. Here we show that high-pressure Pb-S phases stable under lower-mantle conditions provide a plausible silicate-Earth mechanism for long-term Pb storage and periodic remobilisation.

Lead (Pb) has four stable isotopes: 204Pb, 206Pb, 207Pb, and 208Pb. Among these, 204Pb is unradiogenic and often referred to as “common lead”, whereas 206Pb, 207Pb, and 208Pb are radiogenic isotopes produced by the decay of 238U, 235U, and 232Th over varying half-lives. The parent elements U and Th are both highly incompatible in mantle melting, whereas Pb shows broadly similar or slightly lower incompatibility under oxidised conditions, but can behave less incompatibly in the presence of sulfide, owing to its chalcophile affinity. Their differential enrichment between the continental crust and mantle thus produces distinct lead isotopic signatures through time. Consequently, ratios of radiogenic to common Pb (206Pb/204Pb, 207Pb/204Pb, and 208Pb/204Pb) serve as powerful tracers of long-term silicate Earth differentiation, including crust-mantle separation, sediment recycling, and ore formation. However, the Th-U-Pb isotope system is associated with an enduring puzzle, the “missing Pb paradox”. This paradox arises from the observation that the combined Pb isotopic compositions of the continental crust and mantle do not match the expected composition of the BSE, assuming a chondritic starting composition. Specifically, the average Pb isotopic composition of the BSE appears excessively radiogenic, exhibiting elevated 206Pb/204Pb and 207Pb/204Pb ratios. The majority of rocks from the accessible Earth are plotted on or to the right of the geochron, indicating a more radiogenic composition than chondrites.

Given its importance, geoscientists have investigated the missing Pb paradox for decades. The current hypotheses proposed to resolve the missing Pb paradox can be grouped into two categories: (1) the presence of a hidden, unradiogenic Pb reservoir in the core or lower mantle, and (2) a non-chondritic bulk Earth composition, implying deviation from solar system building blocks. For example, a prominent hypothesis suggests that Pb was partitioned into the core through sulfide melt segregation during core formation. However, the rapid formation of the core and the metal-sulfide-silicate partition coefficients of Pb indicate that the core alone could not retain a significant amount of Pb, implying the existence of additional hidden reservoirs with highly unradiogenic Pb in the bulk silicate Earth (i.e., the mantle and crust). Crustal xenolith studies show that the lower crust does not contain enough unradiogenic Pb to balance other reservoirs. While some mantle reservoirs, such as Horoman mantle peridotites and sulfides in abyssal peridotites, contain highly unradiogenic Pb, resolving the Pb paradox requires a considerably larger reservoir, which remains largely unsampled by mantle-derived melts. Some models propose internally consistent Pb isotope evolution curves for the bulk silicate Earth based on specific assumptions about Earth’s history. Despite decades of study, the persistence of the Pb paradox underscores the limitations of our current understanding. In this study, we investigate the nature and mechanism of a hidden reservoir enriched in unradiogenic Pb, under the assumption of chondritic BSE.

The strong chalcophile affinity of Pb promotes its migration and sulfide mineralization, effectively separating it from the parent isotopes of radiogenic Pb and thereby influencing the isotopic composition of the silicate Earth. Galena (PbS), the primary and economically significant Pb-S mineral, crystallizes in the rock-salt structure (B1, Fm3m) and commonly forms in hydrothermal veins from Pb-rich fluids flowing through fractures in the Earth’s crust. A small amount of anglesite (PbSO4) is also found in the oxidized zone of the Pb ore bodies. While Pb-bearing minerals in the deep Earth are less well characterized than those in the crust, mass balance estimates indicate that ~60% of Earth’s Pb resides in the mantle, based on abundances in the bulk silicate Earth and continental crust. Evidence for mantle Pb includes sulfur-bearing xenoliths and hydrothermal fluids. However, specific Pb-hosting mineral phases in the deep Earth have yet to be conclusively identified.

Given the reducing conditions of the deep Earth and the typical oxidation states of Pb and S, PbS and PbS2 are likely candidates for stable Pb minerals at high pressure and high temperature conditions. Experimental studies have shown that PbS undergoes two pressure-induced transitions: from the ambient pressure B1 phase to an intermediate orthorhombic phase at 2.2 GPa and subsequently to a B2 (Pm3m) phase at 21.5 GPa. The intermediate orthorhombic phase can crystallize in either the B33 (Cmcm) or B27 (Pnma) structure, depending on experimental conditions. High-pressure experiments up to 7 GPa have also identified two PbS2 crystal structures, SnS2 (P3m1) and CuAl2 (I4/mcm), which were synthesized through the reaction of PbS and excess sulfur. However, the viability of these compounds at temperatures and pressures relevant to deep Earth conditions remains uncertain.

In light of these results, we conducted an extensive crystal structure search to map the high-pressure phase diagram of Pb-S compounds. In addition to the known PbS, we identified two previously unreported Pb sulfide phases, PbS2 and PbS3, which feature unique polysulfur anionic units and exhibit thermodynamic stability across a broad range of pressure and temperature conditions. The high melting resistance of PbS suggests that common Pb could have been sequestered during early Earth differentiation and retained to the present day, effectively isolating Pb from its parent isotopes U and Th. In sulfur-rich regions, PbS may further react with elemental sulfur to form polysulfides, with PbS2 remaining stable in the upper mantle, while the low-melting PbS3 may facilitate limited upward migration of unradiogenic Pb. These findings could provide a fundamental mechanism for a hidden Pb reservoir.

Results
Stability and structures of Pb1−xSx
The relative thermodynamic stability of Pb-S compounds with various sulfur content was calculated by the equation and shown in the convex hull: ΔH(Pb1−xSx) = H(Pb1−xSx) - (1-x)H(Pb solid) - xH(S solid) (0 < x < 1). Elemental hcp-Pb, S-I, S-II, S-III, and S-IV were chosen as reference structures within their respective pressure ranges. Among the predicted compounds, PbS remains stable at all predicted pressures and exhibits the lowest formation enthalpy. At 10 GPa, PbS2 becomes thermodynamically stable but decomposes at higher pressures into PbS and PbS3. The pressure-composition phase diagram of the predicted compounds is shown in Fig. 1c. Within the studied pressure range (0–150 GPa), PbS undergoes two phase transitions, first from the B1 structure to B33 at 8.5 GPa, and then to the B2 structure at 25.3 GPa. B33 exhibits the lowest enthalpy in this study. The B1-Cmcm-B2 transition sequence reflects lattice-shear-driven symmetry lowering and structural reconfiguration, consistent with previous first-principles transition-path analyses showing that Cmcm acts as a mechanically stable intermediate between the two high-symmetry phases under compression. The calculated B1 to B33 transition occurs slightly higher than the experimental observations. The higher transition pressure obtained here reflects 0 K static enthalpy calculations without kinetic or stress-path effects, whereas experimental transitions at ~2 GPa are sensitive to non-hydrostatic stress and temperature. The high-pressure phase of PbS2, with space group I4/m, is energetically preferable between 4.5 and 14.5 GPa, while PbS3 adopts a P421m structure and remains stable from 11.5 to 42.7 GPa. Extensive structure searches over PbxSy compounds (x=1–4, y=1–4) did not identify any additional thermodynamically stable stoichiometries beyond PbS, PbS2, and PbS3 within the investigated pressure range. No imaginary phonon modes were found at the given pressures, corroborating the dynamical stabilities of those lead polysulfides.

The previously characterized I4/mcm-PbS2 phase is dynamically unstable at 10 GPa, as evidenced by imaginary phonon modes at the Γ point. However, phonon softening reveals a more stable I4/m structure, which involves a rotational distortion of the S2 dimers. Both phases were readily reproduced in the same structural prediction run at 10 GPa. The I4/m-PbS2 phase consists of alternating layers of S2 dimers and PbS8 square antiprisms oriented perpendicular to the (0 0 1) direction. In adjacent layers, the dumbbell-shaped S2 units are arranged perpendicular to one another. With increasing pressure, the S2 units corotate in the same direction, and the angle (θ) between the dimers and the unit cell axis gradually increases, driving a structural transition from the I4/mcm (θ=45∘) to I4/m. This continuous transition implies a second-order phase transition of PbS2, which causes the X-ray diffraction (XRD) characteristic peaks to shift imperceptibly. The subtle effect requires further experimental investigation. The structure of PbS3, which is commonly found in divalent metal trichalcogenides such as CaS3, CaO3, BaS3, SrS3, and MoS3, consists of isolated V-shaped S3 units and face-sharing PbS8 polyhedra. The trisulfur unit comprises a central S1 atom and two terminal S2 atoms. Pb atoms, located at the center and the edges of the structure, form square-planar bonds with four adjacent S1 atoms and cubic coordination with eight S2 atoms. Electron localization function (ELF) analysis reveals a strong electron concentration on the S-S bonds within the polysulfur units, demonstrating homonuclear covalent bonding, a typical characteristic of many high-pressure metal polysulfides. Compression increases orbital overlap within S-S units, widening bonding-antibonding splitting and lowering total energy, thereby stabilising covalent polysulfur units.

Electronic properties of PbS2 and PbS3
Bader charge analysis quantifies electron transfer from Pb to S in these compounds, showing that Pb atoms have similar charge values in all three compounds, which are nearly twice those of the S atoms in PbS and the S2 atoms in PbS2 and PbS3. Bader charge analysis indicates partial electron transfer from Pb to S, consistent with predominantly ionic Pb-S interactions combined with significant covalency. While absolute Bader charges (≈+1 on Pb) are smaller than formal oxidation states, this behaviour is typical for heavy p-block elements and reflects charge delocalisation rather than reduced valence. Formal oxidation states (+II for Pb, -I for S in polysulfide units) are inferred from electron counting within the S2(1-) and S3(2-) anionic motifs, not from absolute Bader charge magnitudes. The near-zero Bader charge on the central S1 atom in PbS3 reflects its symmetric bonding within the S3 unit, analogous to neutral central atoms in molecular trisulfide anions. Consequently, the polysulfide units in these Pb-S compounds can be regarded as isolated S2(n-) anions. Atomic orbital energy level analysis shows that high pressure does not alter the energy sequence of valence electron orbitals for Pb and S atoms (Pb-6p > S-2p > Pb-6s). Thereby, Pb is limited to a maximum oxidation state of +II in sulfides, which prevents the formation of higher oxidation states and more complex stoichiometries. Band structures and projected density of states (PDOS) of PbS2 and PbS3 reveal that both are semiconducting at given pressures. S-p orbitals dominate near the Fermi level, while the Pb-p orbitals are mainly localized at the bottom of the conduction band, suggesting that Pb-p electrons are fully transferred to the polysulfur units. Pressure significantly shortens interatomic distances and strengthens atomic interactions, which leads to band broadening and causes the band gaps of PbS2 and PbS3 to generally decrease with increasing pressure. However, the band gap of PbS3 slightly increases at lower pressures (<15 GPa) as the conduction bands near the Fermi level shift to higher energies. Although band broadening generally reduces band gaps with pressure, competing shifts in conduction-band minima lead to a non-monotonic response in PbS3 at low pressures. Additionally, PbS3 undergoes a transition from an indirect to a direct band gap and then back to an indirect band gap. At 15 GPa, the highest occupied band at the Z points becomes more dispersive, shifting the valence band top from M to Z and resulting in a direct band gap. At higher pressures, the lowest unoccupied band at Z flattens, shifting the conduction band bottom to A, which completes the transition back to an indirect band gap.

P-T diagrams of Pb1−xSx
To evaluate the geological viability of these Pb-S compounds, we analyzed their energetic and structural stability at high temperatures. Phase boundaries for decompositions and phase transitions were determined using Gibbs free energy calculations within the quasi-harmonic approximation. Zero-point energy slightly modifies these boundary pressures at 0 K. Above 14.5 GPa, PbS2 decomposes into PbS and PbS3, with a minor temperature effect on this threshold pressure. At ambient pressure, PbS2 is energetically favorable above 400 K but decomposes into PbS and S at lower temperatures. The decomposition pressure of PbS3 into PbS2 and S increases from 43.7 GPa at 0 K to 50 GPa at 2485 K. Conversely, below 10.2 GPa at 0 K, it decomposes into PbS and S, with this pressure rising to 14.2 GPa at 3500 K. For PbS, the B1 phase transitions to the intermediate B33 phase above 670 K at ambient pressure, which furthur transforms into the B2 phase at pressures above 25.2 GPa. High temperatures accelerate this transition, lowering the transition pressure to 17 GPa at 4000 K. Furthermore, the melting temperatures were calculated using the NVE ensemble via Z-method to assess the physical states of Pb sulfides under geological conditions. The melting points of PbS3 are slightly below the geotherm, suggesting that it primarily exists in a molten state within the Earth. In contrast, PbS2 has melting points around or above the geotherm, indicating its potential stability in crystalline form within the upper mantle and crust. Melting point calculations for the B33 and B2 phases show that PbS has significantly higher melting temperature than geotherm. Even under extreme core-mantle boundary conditions (135 GPa and 5000 K), mean-squared displacement (MSD) calculations using the NVT ensemble show no atomic diffusion in PbS, confirming its crystalline stability throughout the Earth’s crust and mantle.

Methods
The crystal structures of PbxSy compounds (x=1–4, y=1–4) were investigated using the CALYPSO particle swarm optimization method and its corresponding software, at pressures of 0, 10, 30, 50, 100, and 150 GPa, covering the mantle pressure region from surface conditions to the core-mantle boundary. For each Pb-S stoichiometry, primitive cells containing 1, 2, 3, 4, and 6 formula units were considered, generating approximately 1500 candidate structures per cell. Geometric optimization, electronic structure, and phonon calculations were performed within the framework of Density-Functional Theory using the Perdew-Burke-Ernzerhof exchange-correlation functional in the Generalized Gradient Approximation, as implemented in the Vienna Ab initio Simulation Package. The projector augmented-wave method was used to describe the electron-ion interaction, with the 5d10 6s2 6p2 and 3s2 3p4 configurations considered as valence electrons for Pb and S, respectively. A plane wave kinetic energy cutoff of 900 eV and Monkhorst-Pack k-point meshes with a spacing of 2π×0.03 Å−1 were used to ensure convergence of total energies (~1 meV/atom) and residual forces (~10−3 eV/Å). Electronic charge transfer calculations were analyzed using the Bader charge analysis. Phonon dispersion relations were computed via the finite displacement method within the PHONOPY code, to assess the dynamical stability of each structure. Unless otherwise stated, enthalpies refer to static 0 K DFT including the pV term. Zero-point energies and finite-temperature vibrational contributions were incorporated only in Gibbs free-energy calculations using the quasi-harmonic approximation: G(T,p) = minV [U(V) + pV + Fphonon(T,V)], where the minV [functions of V] indicates finding the unique minimum of expression in the brackets by varying the volume, U(V) is the static DFT energy as a function of volume, pV is the pressure-volume term, and Fphonon(T, V) is the vibrational free energy obtained from phonon calculations. Ab initio molecular dynamics (AIMD) simulations were performed using both NVE and NVT ensembles for distinct purposes. Melting point isochores were determined from NVE simulations using the Z method, with temperature controlled by the Andersen thermostat. Simulations employed supercells containing 72, 128, 108, and 96 atoms for Cmcm-PbS, Pm3m-PbS, PbS2, and PbS3, respectively, with Γ point sampling, a 1 fs timestep, and 2000 to 5000 simulation steps for each temperature point. NVT simulations were calculated to evaluate atomic mobility and sublattice diffusion using a 128-atom supercell and Γ point sampling, and were run for at least 7000 steps with a 1 fs timestep using the Nosé-Hoover thermostat.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors systematically calculated the high-pressure phases of PbxS1-x, giving some important information about this\nbinary compound, such as the phase diagram, electronic properties, phonon dispersion, Bader charge and even the reduced\npartition function ratios of Pb. They used CALYPSO to predicted structure at high pressure, which is widely used in the\ncrystal structure prediction and proven correct. Besides, the characters of crystal were calculated using VASP, the most\nfamous theoretical calculation in materials. The sets for DFT calculation are high enough to ensure the results correct, such\nas the plane wave kinetic energy cutoff is set 900 eV, which is also presented in the phonon dispersion calculations where\nis no image frequency. However, I have some confusions:\n1. Geometric optimization, electronic structure, and phonon calculations is ok, using NVE ensemble to calculate the melting\npoint is also easy to understand though the time step is slight large and the simulation time is not so long. However, I do not\nfind the calculation details of NVT MD simulations shown in the Fig S10. Besides, there should be more supplementary\ninformation about the calculation methods and processes, because there used many different methods and predicted many\nproperties.\n2. I am confused the Bader charge calculation in the Fig 3a. The results show charge of Pb is near +1, and the S (PbS2), S2 (PbS3) are near -0.5, why do you say the charge of Pb is +2 and the S (PbS2), S2 (PbS3) are -1? Which is nearly double the\ncalculation results. The author said that the Bader analysis tends to underestimate absolute electron transfer, however, the\nerror should not be so large as the author justified. If the author is right, why do you think the charge of S1 (PbS3) is 0?\n3. For the ∆H calculation, I found no clear explanation about the calculation details. I do not know if the ∆H in this paper is\nstatic or considering the temperature effects. Because the enthalpy (H) has two parts, one is the internal energy (U), the other\nis the volume work (pV). For the theoretical calculation of ∆H, there also should consider the zero-point energy (ε) and the\nphonon contribution if the temperature effect is also considered. Through the paper and supplementary materials, I do not\nfind the description about these.\n4. The author’s calculation show that the PbS is stable until CMB and its melting curve is high than geothermal curve during\nthe whole mantle. However, under the initial CMB condition of earth magma ocean (60 Gpa 4000 K), the PbS is melting,\nwhich might be the mechanism of Pb into core. And the authors also think there should be some melts remain in the mantle\nand not participate in mantle convection during the future evolution of the earth. Thus, they think there should be a reservoir\nof Pb in the lower mantle. I think this deduction exists some troubles. One is the effect of density of PbS melt compared with\nthe surroundings. Because the density dominates the fate of the melts. The other is the immersion properties of the PbS in\nthe surrounding silicate rock. The authors think during the earth evolution, surface Pb is transport by the PbS2 melts.\nHowever, if the melt is immersion (contact angle is less than 90 degree) on the surface of silicate rock, it is ok to forming\nconcentrated melts chamber. If the melt is not immersion (contact angle is larger than 90 degree), these melts should not\nconcentrate together to forming the Pb reservoir.", "reviewer_2": "Review on the paper of Hidden pressure-stabilized lead reservoirs in Earth’s mantle\nIn the present studies, the authors of Siyu Liu et al. using first-principles structural prediction and melting-point calculations,\nthat novel high-pressure Pb–S phases can stably host unradiogenic Pb deep within the mantle. I read through this\nmanuscript very conscientiously. As a whole, the manuscript is well written, extremely well organized, is easy to read and\nwell-illustrated. The data support the conclusion of this study. Thus, I would strongly recommend its publication in Nature\nCommunications after minor revisions. Some indispensable revisions and improvements need to be conducted as follows:\n[1] In PbS, B1-B2 transition type is common. However, why such cubic structures have chosen the inter-mediate transition\npath like B1-Cmcm-B2 under high pressure regimes?\n[2] According to the Ref. 32, Cmcm phase takes place at just 2 GPa. But, in your experimental conditions, the Cmcm phase\nexist 8 GPa. Please justify the possible reasons for this critical pressure enhancement.\n[3] “As pressure increases, the energy gap between the bonding and antibonding states of these S-S covalent bonds\nwidens, which improves the bond strength and structural stability.” Authors can clearly explain how during the covalent\nbonds widens, the bond strength and structural stability can be enhanced.\n[4] “Pressure significantly shortens interatomic distances and strengthens atomic interactions, which leads to band\nbroadening and causes the band gaps of PbS2 and PbS3 to generally decrease with increasing pressure”. The above-\nmentioned statement is quite opposite trend for the band gap energy with respect to pressure increases. The authors are\nrequested make a clarification about this band gap changes with respect to pressure in PbS, PbS2 and PbS3.\n[5] Can PbS exhibits different stoichiometry phases under high pressure-temperature similar to NaCl\n(10.1126/science.1244989)?"}, "author_response": "Point by point response to reviewers’ comments\nBelow, we provide a point-by-point response to all reviewer comments, reproduced\nverbatim and followed by our responses in italics. All changes have been implemented\nin the revised manuscript and Supplementary Information, with tracked changes\nhighlighted.\nReviewer(s)' Comments to Author:\nReviewer #1:\nRecommendation: The authors systematically calculated the high-pressure phases of\nPb S , giving some important information about this binary compound, such as the\nx 1-x\nphase diagram, electronic properties, phonon dispersion, Bader charge and even the\nreduced partition function ratios of Pb. They used CALYPSO to predicted structure at\nhigh pressure, which is widely used in the crystal structure prediction and proven\ncorrect. Besides, the characters of crystal were calculated using VASP, the most\nfamous theoretical calculation in materials. The sets for DFT calculation are high enough to\nensure the results correct, such as the plane wave kinetic energy cutoI is set 900 eV,\nwhich is also presented in the phonon dispersion calculations where is no image\nfrequency.\nHowever, I have some confusions:\nResponse:\nWe appreciate the constructive suggestions provided and have carefully revised the\nmanuscript to address all the comments.\n1. Geometric optimization, electronic structure, and phonon calculations is ok, using\nNVE ensemble to calculate the melting point is also easy to understand though the time\nstep is slight large and the simulation time is not so long. However, I do not find the\ncalculation details of NVT MD simulations shown in the Fig S10. Besides, there should\nbe more supplementary information about the calculation methods and processes,\nbecause there used many diIerent methods and predicted many properties.\nResponse:\nWe agree and thank the reviewer for pointing out this lack of clarity. In the original\nsubmission, the distinction between NVE and NVT ab initio molecular dynamics (AIMD)\nsimulations were not made explicit.\nWe have now clarified this by:\n1. Explicitly stated that NVT simulations were used solely to evaluate atomic\nmobility and sublattice diWusion whereas NVE simulations were used for\nmelting-point determination via the Z-method.\nThe revised statement is “Furthermore, the melting temperatures were\ncalculated using the NVE ensemble via Z-method to assess the physical states of\nPb sulfides under geological conditions (Fig. S9). … Even under extreme core-\nmantle boundary conditions (135 GPa and 5000 K), mean-square displacement\n(MSD) calculations using the NVT ensemble show no atomic diIusion in PbS,\nconfirming its crystalline stability throughout the Earth’s crust and mantle (Fig.\nS10).” (Page 4, last paragraph)\n2. Separated NVE and NVT AIMD simulations and provided full details of the NVT\nsimulations (ensemble, thermostat, timestep, simulation length, supercell size)\nin the Methods section of main text.\nThe added text is “Ab initio molecular dynamics (AIMD) simulations were\nperformed using both NVE and NVT ensembles for distinct purposes. Melting\npoint isochores were determined from NVE simulations using the Z method47,\nwith temperature controlled by the Andersen thermostat48. Simulations\nemployed supercells containing 72, 128, 108, and 96 atoms for Cmcm-PbS,\nPm3m-PbS, PbS , and PbS , respectively, with Γ point sampling, a 1 fs time step,\n2 3\nand 2000 to 5000 simulation steps for each temperature point. NVT simulations\nwere calculated to evaluate atomic mobility and sublattice diIusion using a 128-\natom supercell and Γ point sampling, and were run for at least 7000 steps with a\n1 fs time step using the Nos𝑒́-Hoover thermostat49.” (Page 3, first paragraph)\n3. Added a new Supplementary Section (“Computational workflow”), which\nprovides a step-by-step overview of all computational methods employed,\nincluding structure prediction, static DFT, phonons, quasi-harmonic free\nenergies, NVE melting simulations, and NVT diWusion analyses.\nThese additions improve transparency and reproducibility.\n2. I am confused the Bader charge calculation in the Fig 3a. The results show charge of\nPb is near +1, and the S (PbS2), S2 (PbS3) are near -0.5, why do you say the charge of Pb\nis +2 and the S (PbS2), S2 (PbS3) are -1? Which is nearly double the calculation results.\nThe author said that the Bader analysis tends to underestimate absolute electron\ntransfer, however, the error should not be so large as the author justified. If the author is\nright, why do you think the charge of S1 (PbS3) is 0?\nResponse:\nWe thank the reviewer for highlighting this important conceptual issue. We agree that\nthe original wording could be misinterpreted as equating Bader charges directly with\nformal oxidation states, which was not our intention.\nWe have therefore revised the corresponding sentences to explicitly distinguish\nbetween Bader charges and formal oxidation states, which are fundamentally diWerent\nquantities. Bader charges represent a real-space partitioning of continuous electron\ndensity and are known to underestimate formal valence, particularly for heavy p-block\nelements with significant covalency and charge delocalization. We have removed the\nprevious phrasing that could be read as assigning formal oxidation states from Bader\ncharges. The revised text now reports Bader charges as qualitative indicators of charge\nredistribution, while formal valence is discussed separately via electron counting of\npolysulfide units.\nIn the revised manuscript we now:\n1. Avoid equating absolute Bader charge values with oxidation states.\n2. Clarify that formal oxidation states are inferred independently from electron-\ncounting arguments applied to the S 2- and S 2- polysulfide units, not from the\n2 3\nnumerical magnitude of Bader charges.\n3. Explicitly explain that the near-zero Bader charge on the central S1 atom in PbS\n3\nreflects its symmetric bonding environment within the trisulfide unit, analogous\nto neutral central atoms in molecular polysulfides.\nThe relevant paragraph in the Results section has been rewritten to remove ambiguity\nand prevent conceptual misunderstanding.\nThe added statement is “Bader charge analysis indicates partial electron transfer from\nPb to S, consistent with predominantly ionic Pb–S interactions combined with\nsignificant covalency. While absolute Bader charges (≈+1 on Pb) are smaller than\nformal oxidation states, this behaviour is typical for heavy p-block elements and reflects\ncharge delocalization rather than reduced valence. Formal oxidation states (+II for Pb, -I\nfor S in polysulfide units) are inferred from electron counting within the S 2− and S 2−\n2 3\nanionic motifs, not from absolute Bader charge magnitudes. The near-zero Bader charge\non the central S1 atom in PbS3 reflects its symmetric bonding within the S unit,\n3\nanalogous to neutral central atoms in molecular trisulfide anions.” (Page 4, first\nparagraph)\n3. For the ∆H calculation, I found no clear explanation about the calculation details. I do\nnot know if the ∆H in this paper is static or considering the temperature eIects. Because\nthe enthalpy (H) has two parts, one is the internal energy (U), the other is the\nvolume work (pV). For the theoretical calculation of ∆H, there also should consider the\nzero-point energy (ε) and the phonon contribution if the temperature eIect is also\nconsidered. Through the paper and supplementary materials, I do not find the\ndescription about these.\nResponse:\nWe agree that this distinction was insuWiciently explicit in the original manuscript.\nWe have now clarified that:\n• Formation enthalpies (ΔH) shown in convex-hull diagrams refer to static 0 K DFT\nenthalpies including the pV term only.\n• Zero-point energy and finite-temperature vibrational contributions are included\nexclusively in Gibbs free-energy calculations, which were performed using the\nquasi-harmonic approximation and used to construct the P-T phase diagrams.\nThis clarification has been added to the Methods section, together with an explicit\nexpression for the Gibbs free energy and a description of where each term is applied.\nThe added text is “Unless otherwise stated, enthalpies refer to static 0 K DFT including\nthe pV term. Zero-point energies and finite-temperature vibrational contributions were\nincorporated only in Gibbs free-energy calculations using the quasi-harmonic\napproximation: 𝐺(𝑇,𝑝) = 𝑚𝑖𝑛[𝑈(𝑉)+𝑝𝑉 +𝐹 (𝑇,𝑉)], where the\n\"#$%$%\n!\n𝑚𝑖𝑛[𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛𝑠 𝑜𝑓 𝑉] indicates finding the unique minimum of expression in the\n!\nbrackets by varying the volume, 𝑈(𝑉) is the static DFT energy as a function of volume,\n𝑝𝑉 is the pressure–volume term, and 𝐹 (𝑇,𝑉) is the vibrational free energy\n\"#$%$%\nobtained from phonon calculations46.” (Page 3 first paragraph)\n4. The author’s calculation show that the PbS is stable until CMB and its melting curve is\nhigh than geothermal curve during the whole mantle. However, under the initial CMB\ncondition of earth magma ocean (60 GPa 4000 K), the PbS is melting, which\nmight be the mechanism of Pb into core. And the authors also think there should be\nsome melts remain in the mantle and not participate in mantle convection during the\nfuture evolution of the earth. Thus, they think there should be a reservoir of Pb in the\nlower mantle. I think this deduction exists some troubles. One is the eIect of density\nof PbS melt compared with the surroundings. Because the density dominates the fate\nof the melts. The other is the immersion properties of the PbS in the surrounding silicate\nrock. The authors think during the earth evolution, surface Pb is transport by the PbS2\nmelts. However, if the melt is immersion (contact angle is less than 90 degree) on the\nsurface of silicate rock, it is ok to forming concentrated melts chamber. If the melt is\nnot immersion (contact angle is larger than 90 degree), these melts should not\nconcentrate together to forming the Pb reservoir.\nResponse:\nWe thank the reviewer for this important point of melt buoyancy and wetting angle,\nwhich indeed control whether a minor phase can form an interconnected melt network\nand segregate to form a physically “gathered” reservoir. We realize that our wording may\nhave unintentionally implied such a scenario, and we have clarified this in the revised\nmanuscript.\nFirst, we would like to emphasize that we do not state or require that Pb forms a distinct,\npooled melt reservoir in the lower mantle (e.g., a concentrated melt chamber). Our\nmain inference is mineralogical rather than mechanical: because the calculated\nstability field and melting curve of PbS (and related Pb–S phases) are high relative to the\nmantle geotherm over much of the mantle, Pb hosted in these phases can persist as a\nrefractory component during subsequent mantle evolution and does not readily enter\nsilicate melts during decompression melting. In this sense, the “reservoir” is simply Pb\nretained in stable high-melting phases (potentially dispersed, grain-scale, or\nheterogeneously distributed), and it need not be physically segregated or “gathered” to\nexplain reduced Pb extraction.\nSecond, regarding the reviewer’s point that PbS may melt under very early CMB magma-\nocean conditions, we agree that this could facilitate Pb partitioning toward the core\nduring the magma-ocean stage. However, our argument about later Earth evolution\nconcerns the persistence of Pb in high-melting sulfide/sulfosalt phases after magma-\nocean crystallization and during subsequent mantle processing. The existence of\ntransient melting at the CMB in the earliest stage does not require that later Pb must\nremain as an interconnected melt phase; Pb can be stored in solid phases that are\nstable and refractory under subsolidus mantle conditions.\nTo address this concern, we have narrowed and clarified the scope of our conclusions.\nSpecifically, we have:\n1. Added an explicit statement in the Discussion acknowledging that we do not\nrequire melt segregation and therefore do not rely on assumptions about melt\ndensity or wetting angles.\n2. Revised our wording throughout the manuscript to emphasise thermodynamic\nfeasibility rather than inevitability.\n3. Replaced deterministic language (e.g., “will form a reservoir”, “provides the\nmechanism”) with probabilistic language (“could form a reservoir”, “provides a\nplausible mechanism”), and the changes have been highlighted in blue.\nThe revised text now makes clear that our results demonstrate that pressure-stabilised\nPb–S phases can exist and remain stable under mantle conditions, thereby oWering a\nplausible mineralogical mechanism for Pb sequestration, without claiming to uniquely\ndetermine melt segregation eWiciency or exclude other mechanisms.\nThe revised statement is “The long-term fate of sulfide melts in the mantle depends on\ntheir density contrast with surrounding silicates and their wetting behaviour along grain\nboundaries, parameters not explicitly calculated here. … Our results therefore\ndemonstrate the thermodynamic feasibility of deep Pb sequestration in sulfide phases,\nrather than uniquely determining melt segregation eIiciency.” (Page 5 second\nparagraph)\nReviewer #2:\nRecommendation: Review on the paper of Hidden pressure-stabilized lead reservoirs in\nEarth’s mantle\nIn the present studies, the authors of Siyu Liu et al. using first-principles structural\nprediction and melting-point calculations, that novel high-pressure Pb–S phases can\nstably host unradiogenic Pb deep within the mantle. I read through this manuscript very\nconscientiously. As a whole, the manuscript is well written, extremely well organized, is\neasy to read and well-illustrated. The data support the conclusion of this study. Thus, I\nwould strongly recommend its publication in Nature Communications after minor\nrevisions. Some indispensable revisions and improvements need to be conducted as\nfollows:\nResponse:\nWe sincerely thank the reviewer for the positive evaluation of our manuscript and for\nrecognizing the significance of our findings. We have revised several statements in the\nmanuscript to enhance its clarity and readability.\n[1] In PbS, B1-B2 transition type is common. However, why such cubic structures have\nchosen the inter-mediate transition path like B1-Cmcm-B2 under high pressure\nregimes?\nResponse:\nWe have clarified that the B1–Cmcm–B2 sequence reflects a symmetry-lowering shear\ndistortion pathway that is mechanically favourable under compression. This behaviour\nis consistent with prior experimental and theoretical studies of PbS and related\nchalcogenides. A brief explanatory sentence and appropriate citations have been\nadded to the Results section.\nThe added discussion is “The B1–Cmcm–B2 transition sequence reflects lattice-shear-\ndriven symmetry lowering and structural reconfiguration, consistent with previous first-\nprinciples transition-path analyses showing that Cmcm acts as a mechanically stable\nintermediate between the two high-symmetry phases under compression36.” (Page 3\nsecond paragraph)\n[2] According to the Ref. 32, Cmcm phase takes place at just 2 GPa. But, in your\nexperimental conditions, the Cmcm phase exist 8 GPa. Please justify the possible\nreasons for this critical pressure enhancement.\nResponse:\nWe have added a clarification noting that our transition pressures are derived from\nstatic 0 K enthalpy calculations under ideal hydrostatic conditions. Experimental\ntransition pressures are sensitive to temperature, stress state, and kinetic eWects,\nwhich can lower observed transition pressures. This explanation has been added to the\nmanuscript.\nThe added clarification is “The calculated B1 to B33 transition occurs slightly higher\nthan the experimental observations56. The higher transition pressure obtained here\nreflects 0 K static enthalpy calculations without kinetic or stress-path eIects, whereas\nexperimental transitions at ∼2 GPa are sensitive to non-hydrostatic stress and\ntemperature.” (Page 3 second paragraph)\n[3] “As pressure increases, the energy gap between the bonding and antibonding states\nof these S-S covalent bonds widens, which improves the bond strength and structural\nstability.” Authors can clearly explain how during the covalent bonds widens, the bond\nstrength and structural stability can be enhanced.\nResponse:\nWe have added a concise explanation stating that increased orbital overlap under\ncompression widens the bonding–antibonding splitting, lowering the total electronic\nenergy and stabilising the covalent S–S units. This clarification has been incorporated\ninto the Results section.\nThe added explanation is “Compression increases orbital overlap within S–S units,\nwidening bonding–antibonding splitting and lowering total energy, thereby stabilising\ncovalent polysulfur units65.” (Page 3 second paragraph)\n[4] “Pressure significantly shortens interatomic distances and strengthens atomic\ninteractions, which leads to band broadening and causes the band gaps of PbS and\n3\nis quite opposite trend for the band gap energy with respect to pressure increases. The\nauthors are requested make a clarification about this band gap changes with respect to\npressure in PbS, PbS and PbS .\n2 3\nResponse:\nWe have reorganised this section to make clear that although band broadening\ngenerally reduces band gaps with pressure, competing shifts in band extrema lead to\nnon-monotonic behaviour in PbS at low pressures. The explanation has been moved\n3\nearlier in the paragraph to avoid confusion.\nThe revised explanation is “Although band broadening generally reduces band gaps with\npressure, competing shifts in conduction-band minima lead to a non-monotonic\nresponse in PbS at low pressures.” (Page 4, first paragraph)\n3\n[5] Can PbS exhibits diIerent stoichiometry phases under high pressure-temperature\nsimilar to NaCl (10.1126/science.1244989)?\nResponse:\nWe have clarified that extensive structure searches up to Pb S (x=1-4, y=1-4)\nx y\nstoichiometries did not reveal additional thermodynamically stable phases beyond\nPbS, PbS , and PbS within the investigated pressure range. This statement has been\n2 3\nadded to the Results section.\nThe added statement is “Extensive structure searches over Pb S compounds (x = 1–4, y\nx y\n= 1–4) did not identify any additional thermodynamically stable stoichiometries beyond\nPbS, PbS , and PbS within the investigated pressure range.” (Page 3 second paragraph)\n2 3"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thanks for your explanation about my confusions, I think the MS is complete now. Overall, this is a systematic work for the fate of Pb in the earth. Although this is just a possible mechanism, I think this work will spark a deep discussion about Pb."}, "author_response": "We thank the reviewer for the positive feedback."}]

Decision letter