Pressure, as an external control method, can stabilize exotic compounds that defy the traditional understandings of stoichiometry, structural form, and properties in ambient conditions.1–4 One of the most representative cases of this phenomenon is the burgeoning field of high-pressure superconducting hydrides, which have remarkably high critical temperatures (T_c > 200 K).5–9 Within these compounds, H atoms or H and other atoms can form rigid covalent motifs, exhibiting diverse coordination configurations such as symmetric S–H–S bonds (twofold coordination) in H₃S, 1D nanotubes via threefold coordination in HfH₉,10 2D penta-graphene containing twofold and threefold coordination in HfH₁₀,11 and 3D cages with fourfold coordination in LaH₁₀ and CaH₆.12–14 These coordination patterns are much beyond the bonding limit of single-electron H atoms at ambient pressure, and their corresponding covalent bonds thus possess fewer electrons than a conventional shared electron pair, appearing to be unsaturated. Strikingly, these kinds of H frameworks have been demonstrated to be able to induce strong electron–phonon coupling (EPC).15 Thus far, the understanding of unsaturated-bonding-related superconductivity is quite limited. However, bonding unsaturation could provide a new approach to the design of superconductors.

Phosphorus (P) and sulfur (S), as adjacent elements in the Periodic Table, have moderate electronegativity as well as a strong ability to form rich P/S allotropes16–19 and P–S compounds.20–24 Interestingly, the majority of P–S compounds exhibit cage-like structures,21–24 and this bestows them with intriguing properties and gives them a broad range of potential applications.25,26 For instance, P₄S₃, which has strong resistance to volume changes during electrochemical reactions, is considered as an ideal anode material with high capacity and cycling stability for Li⁺/Na⁺/K⁺ storage, providing a reversible capacity of up to 1266 mAh g⁻¹ in Li–ion batteries.25 P₄S₁₀, having high chemical and electrochemical stability, can act as a protective layer for anodes, overcoming the corrosion and dendrite problems of Li metal anodes.26

In ambient conditions, P and S, which have 3s²3p³ and 3s²3p⁴ valence electrons, usually show three-/two-coordination in their covalent compounds, satisfying the well-known octet rules. However, in their high-pressure allotropes and compounds, their coordination numbers can be boosted to 8 and 6, respectively,17–19,27 resulting in the formation of diverse covalent frameworks with notable bonding unsaturation, contrasting with the classical covalent bonds formed by sharing an electron pair. This unique bonding means that P and S are two significant candidate elements for high-T_c superconductors.28–31

The first experimental validation of a superconducting hydride (T_c > 200 K) was conducted with H₃S, which has a record T_c of 203 K at 155 GPa,32 and this consists of nested SH₆ octahedral networks with each other.4 Moreover, monoclinic (C2/m) PH₃, with a P–H covalent framework and six-coordinated P, exhibits a maximum T_c of 83 K at 200 GPa in P-H compounds.30 Based on the difference in valence electrons, the S atoms in H₃S can be partially substituted with non-metallic atoms from the IV-, V-, VI-, or VII-A groups,33–36 and this can be used to modulate the position of the van Hoff singularity with respect to the Fermi level (E_F), raising the T_c value. In fact, the bonding unsaturation in doped H₃S can also be tuned. Interestingly, the T_c value has been found to increase significantly to 290 K at 250 GPa after replacing 7.5% of the S atoms with P atoms.34 To date, numerous binary and ternary P- and S-based compounds have been predicted to exhibit excellent superconductivity associated with their various covalent networks.37–41

Pressure-induced P–S compounds have also attracted great research interest due to the unique bonding behavior of P and S atoms. In particular, two pioneering works explored the exotic properties of P–S compounds at high pressure.42,43 Li et al. found that as the pressure increases, the stable phases (e.g., P₂S₇ and P₄S_x, x = 3–10) at 0 GPa gradually decompose, and above 8 GPa, the only predicted stable phase is R-3m PS₂; this was confirmed by experiments.42 More interestingly, R-3m PS₂, which is isostructural to 3R MoS₂, is a van der Waals layered compound with a predicted T_c of ∼11 K at 0 GPa. In addition, Liu et al. found several P–S compounds including P₃S, P₂S, PS, PS₂, and PS₃, showing abundant S- and P-centered polyhedrons (e.g., SP₆, SP₄, SP₃, PS₆, and PS₃).43 However, there is a difference regarding the stability of high-pressure PS₂ phases: C2/m and P-3m1 PS₂ were found to stable from 5 to 200 GPa in Liu et al.’s work,43 but the synthesized R-3m phase disappears.42 In summary, the veil over the nature and behavior of P–S compounds at high pressure has not yet been fully lifted.

Taking into account the peculiar bonding behavior of P and S atoms and the imperfect high-pressure P–S phase diagram, herein, we perform a comprehensive structural search with various stoichiometries P_xS_y (x/x + y = 0–1, 38 compositions) to explore as-yet untapped P–S compounds and their electronic properties. Three unknown P-rich phases, Cmcm P₅S, P-31m P₈S, and P-1 P₁₁S, were identified, in which both P and S atoms present sixfold octahedral coordination with unsaturated covalent bonds, inducing superconductivity. Furthermore, in P-rich P–S compounds and several MH₉ hydrides, a positive correlation was established between unsaturated bonding and the superconducting transition temperature, indicating the feasibility of modulating bonding saturation to improve superconductivity. In addition to the three proposed PS₂ phases, a new low-pressure P-3m1 PS₂ phase was found, representing a missing link in the phase-transition sequence with increasing pressure.

An intelligence-based particle swarm optimization algorithm using the CALYPSO software package has been demonstrated to be an important way to discover new structures.7 In this work, structures with various P_xS_y compositions (x = 1, y = 1/12, 1/11-1/3, 3/8, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 1, 5/4, 4/3, 3/2, 5/3, 7/4, 2, 9/4, 5/2, 8/3, 3, 7/2, 4, 5-12) were extensively searched with simulation cell sizes of 1–4 formula units at 25, 50, 100, and 200 GPa. Structural relaxation and electronic structure calculations were carried out within framework of density functional theory with the Perdew–Burke–Ernzerhof generalized gradient approximation,44 as implemented in the Vienna Ab initio Simulation Package.45 Projector augmented waves with 3s²3p³ and 3s²3p⁴ electrons were used as valences for P and S atoms, respectively.46 The reliability of our adopted pseudopotentials was confirmed by the perfect fitting of the Birch–Murnaghan equation of states with the full-potential linearized augmented-plane-wave method [Fig. S3(a)]. A plane-wave cutoff energy of 600 eV and a Monkhorst–Pack scheme with a k-point grid of 2π × 0.03 Å⁻¹ in the Brillouin zone were used to ensure an enthalpy convergence of less than 1 meV/atom. The electron localization function (ELF) is used to describe and visualize the chemical bonds of the stable compounds.47 The integrated crystal orbital Hamilton population (ICOHP) was calculated using the LOBSTER software package.48 Phonon-dispersion curves were calculated using the supercell finite-displacement method in the Phonopy software package.49

The EPC calculations were carried out within the framework of density-functional perturbation theory, as executed in the Quantum ESPRESSO package.50 We employed ultrasoft pseudopotentials, with 3s²3p³ and 3s²3p⁴ considered as the valence electrons for P and S, respectively. A kinetic-energy cutoff of 70 Ry and Gaussians of width 0.02 Ry were applied. To reliably calculate the EPC in metallic systems, we need to sample dense k-meshes for the electronic Brillouin zone integration and enough q-meshes to evaluate the average contributions from the phonon modes. Dependent on specific structures of the stable compounds, different sizes of k- and q-meshes were used: 15 × 15 × 9 k-meshes and 5 × 5 × 3 q-meshes for P₂S in the C2/m-I structure, 12 × 12 × 3 k-meshes and 3 × 3 × 3 q-meshes for P₃S in the C2/m structure, 15 × 15 × 6 k-meshes and 5 × 5 × 2 q-meshes for P₅S in the Cmcm structure, 8 × 8 × 12 k-meshes and 4 × 4 × 6 q-meshes for P₈S in the P-31m structure, and 12 × 12 × 9 k-meshes and 4 × 4 × 3 q-meshes for P₁₁S in the P-1 structure.

To explore the phase diagram of the P–S system, an extensive structure search was performed by focusing on the 38 stoichiometries. The relative stabilities of the P_xS_y compounds were determined using their formation enthalpies with respect to the solid phases of P27 and S18,51 at the considered pressures (Fig. S1). We not only reproduced the theoretically proposed42,43C2/m P₃S, C2/m-I P₂S, C2/m-II P₂S, P2₁/m PS, R-3m PS₂, C2/m PS₂, P-3m1 PS₂, and C2/c PS₃ structures at high pressures, but also identified three hitherto unknown P-rich compounds: Cmcm P₅S, P-31m P₈S, and P-1 P₁₁S [Fig. 1(a)], as well as a new P-3m1 PS₂ phase (labeled as P-3m1-I, Fig. S2).

The optimized lattice parameters of R-3m PS₂ (a = 3.165 Å and c = 15.572 Å) were found to be in good agreement with the experimentally obtained values (a = 3.149 Å, and c = 15.501 Å) at 14 GPa;42 additionally, the use of different exchange–correlation functionals did not affect the thermodynamic stabilities of the phases proposed here [Fig. S3(b)], confirming the reliability of our calculations. Both C2/m and C2/c PS₃ were predicted by Liu et al. to be thermodynamically stable.43 However, in our work, we found that only the C2/c phase is thermodynamically stable. In other words, C2/m PS₃ is unstable against decomposition into S and PS₂ in our considered pressure range of 25–200 GPa, as shown in Fig. S2(b). All the predicted structures are dynamically stable due to the absence of any imaginary frequencies in the phonon dispersion curves (Fig. S4). The stabilizing pressure ranges of the identified phases are illustrated in Fig. 1(b), and this provides useful information for experiments.

There are obvious differences in the stable stoichiometries and structural characters between high and ambient pressures. Most of the P–S compounds are S-rich stoichiometries with the structural features of a cage-like molecular crystal in P₄S_x (x = 5, 7, 9, and 10) and infinite wrinkled chains in P₂S₇ at ambient pressure.20–24 Increasing the S content in the compounds induces a gradual transformation of coordination configuration from triangular cone to tetrahedron for the P atoms, and this is accompanied by an increase of the coordination number (Fig. S5). In sharp contrast, at high pressure, P-rich stoichiometries are stabilized (P_xS, x = 5, 8, and 11). Notably, these P-rich compounds exhibit similar P–S frameworks, including the reported P₂S and P₃S, which can be viewed as derivatives of S-doped simple cubic (sc) P [Figs. 2(a)–2(c)]. In other words, both the P and S atoms in these structures maintain a six-coordinated octahedron configuration independent of the P/S ratio, but the S-centered octahedra demonstrate different coordination atoms, such as four P atoms and two S atoms in P₂S, or mixed four P and two S atoms as well as six P atoms in P₃S, as well as six P atoms in P_xS (x = 5, 8, and 11) at 100 GPa (Fig. S6).

The P/S ratio can be used to effectively modulate the orientation and distribution of S-centered octahedra. Specifically, in P₂S, P₃S, and P₁₁S, S-centered octahedra are distributed horizontally with S atoms in the same layer, leading to the formation of alternating P layers perpendicular to the horizontal plane [Figs. S6(a), S6(b), and S6(f)]. In P₅S and P₈S, the central S atoms are alternately distributed in different directions, inhibiting the formation of the P layer [Figs. S6(c)–S6(e)]. As a result, different connections appear between the S-centered octahedra, including mixed edge/vertex sharing in P₂S and P₃S, vertex sharing in P₅S, and P₂ pair sharing in P₈S and P₁₁S. In the P_xS (x = 5, 8, and 11) compounds, the P–P and P–S bond lengths are in the range 2.12–2.17 Å at 100 GPa, which is close to the P–P (2.20 Å) and P–S (2.09 Å) bond lengths in black phosphorus and P₄S₃ at ambient pressure,52,53 indicating that covalent P–P/S bonds are responsible for the stabilization, as demonstrated by the ELF [Fig. 4(a) and Fig. S7].

As noted in the introduction to this paper, the previously reported stable phases of PS₂ are discrepant under pressure. Aside from the R-3m phase being stable in the range 8–20 GPa,42 the C2/m and P-3m1 phases are found to be stable in the ranges 5–178 and 178–200 GPa,43 respectively. Here, we not only reproduce the stable R-3m, C2/m, and P-3m1 phases, but we also identify a new P-3m1-I phase, which stabilizes from 36.3 to 87.6 GPa; this is an intermediate phase between the R-3m and C2/m phases. Notably, the discovery of the P-3m1-I phase settles the contradiction regarding the PS₂ phase diagram and provides a plausible description of the phase-transition process of PS₂ with pressure.

For R-3m PS₂ with an ABC stacking order, pressure makes the sublayers slide along ab plane into AA stacking, accompanying the transition from the R-3m phase to the P-3m1-I phase (Fig. S8). With further compression, the sublayers in the P-3m1-I phase move away from each other along the a axis until transitioning into C2/m symmetry with staggered sublayers in the vertical direction [Figs. 2(d) and 2(e)]. Then, the PS₂ sublayers in the C2/m phase move toward each other along the a axis with increasing pressure. Finally, AA stacking P-3m1-II PS₂ is obtained above 168.3 GPa [Figs. 2(e) and 2(f)], in which both S and P atoms form a covalent octahedron coordination (e.g., S-centered P₃S₃ and PS₆). Overall, the lattice parameters and distances between S atoms of adjacent sublayers progressively decrease with pressure (Fig. S9). Meanwhile, the charge is redistributed and the interlayered S–S interaction is enhanced, as supported by the ELF (Fig. S10), resulting in a transformation of PS₂ from a van der Waals layered structure to a 3D covalent framework.

One prominent structural characteristic is the formation of covalent octahedron coordination for P and S atoms in P-rich P–S compounds at high pressure. This is mainly attributed to the head-to-head overlap of three p-orbital wave functions [Fig. 3(a)], as verified below. Here, we follow the principle from simplicity to complexity. First, sc P is selected due to the identical orthometric basis vector between the crystal structure and the covalent bond orientation. The calculated ICOHP exhibits dominant bonding interactions between the p_x and p_x, p_y and p_y, and p_z and p_z orbitals along the three bonding directions (Fig. S11). Meanwhile, P 3s electrons are involved in bonding with the 3p electrons of corresponding orientation. A similar bonding pattern is also demonstrated in Im-3m H₃S (Fig. S12). Similarly, P_xS compounds follow the bonding pattern of P₈S (Figs. S13 and S14) by excluding the basis-vector difference between the atomic p orbitals and the covalent bond orientations.

The corresponding covalent bonds for P atoms in these P-rich compounds have unsaturated bonding due to there being five valence electrons from each P atom for six bonds. Here, we propose a standard for evaluating the unsaturation of covalent bonds in these compounds. Each S or P atom has 6 or 5 valence electrons, respectively. On average, a single S atom can provide one electron per covalent bond in a six-coordinated octahedron, whereas a P atom can provide 5/6 of an electron. Therefore, the number of electrons is (1 + 1)e⁻ in each S–S bond, (1 + 5/6)e⁻ in each P–S bond, and (5/6 + 5/6)e⁻ in each P–P bond; the S–S bonds are saturated (e.g., a two-centered two-electron bond), and the other two are unsaturated.

For ease of measuring the covalent unsaturation of these three kinds of bonds, the saturation of each S–S bond is normalized to 1—meaning that the saturations of P–S and P–P bonds are 0.917 and 0.833, respectively—and the bonding unsaturation is defined to be the difference between 1 and the saturation, e.g., 0, 0.083, and 0.167 for S–S, P–S, and P–P bonds, respectively. Correspondingly, the average bonding unsaturation of these five P-rich P–S compounds is 0.106 for P₂S, 0.125 for P₃S, 0.139 for P₅S, 0.148 for P₈S, and 0.153 for P₁₁S [Fig. 3(b)], the detailed calculations can be found in supplementary material. It is apparent that the average bonding unsaturation increases monotonously with increasing P content.

The calculated projected densities of states (PDOS) illustrate that P₅S, P₈S, and P₁₁S are metallic on account of the large electronic contributions at E_F [Fig. 4(b) and Fig. S15]. The considerable orbital overlap between P and S atoms (3p–3p, and 3s–3s) also supports their strong covalent interaction. With respect to P/S 3s states, the 3p states have a much larger DOS at E_F, dominating the metallicity of these P–S compounds. Notably, their DOS amplitude at E_F (NEF) is consistent with the average bonding unsaturation. i.e., the larger the average bonding unsaturation, the higher the PDOS amplitude [Fig. 3(d)].

On account of the identical coordination modes and distinct bonding unsaturation values caused by different P/S ratios in P-rich compounds (e.g., P₂S, P₃S, P₅S, P₈S, and P₁₁S), they provide an ideal model for the study of bonding-unsaturation-related superconductivity. To accurately describe this relationship, we calculated their T_c values at an identical pressure (e.g., 100 GPa) by employing the Allen–Dynes-modified McMillan equation.54 Here, the T_c values of the previously proposed C2/m-I P₂S and C2/m P₃S structures are recalculated for better comparability. The resulting T_c values (2.63 K for C2/m-I P₂S, 3.06 K for C2/m P₃S) are in accordance with the reported values (Fig. S17).43 Interestingly, the value of T_c increases monotonously with the increasing average bonding unsaturation from P₂S to P₁₁S, indicating an evident positive correlation between bonding unsaturation and T_c values [Fig. 3(c)]. Furthermore, their λ values present the same evolution trend with T_c values; the trend in ω_log is opposite, manifesting a decisive effect of λ on the variation of T_c value.

The increasing λ value can be principally ascribed to the elevated NEF value, especially for P-derived DOS contributions at E_F [Fig. 3(d)]. This is related to the increase of the number of P–P bonds, which have greater bonding unsaturation than P–S bonds. However, the variation of the highest phonon frequency of P_xS compounds is inverse to the average bonding unsaturation, which may be derived from the lengthening of P–S and P–P bonds [Fig. 3(e)]. The lowering of the phonon frequency leads to a reduction in ω_log. However, the decreasing phonon frequency ω is also favorable to raising λ according to its definition formula λ=2∫dωα2F(ω)ω.55,56 Overall, a higher bonding unsaturation favors the superconductivity of P_xS compounds.

To further examine the effect of bonding unsaturation on superconductivity, we inspected reported superconducting TM hydrides based on these criteria: they must have the same H cage, different electrons donated by TM atoms that affect bonding saturation, and pressure differences within 30 GPa to minimize the effects of pressure. Among the three well-known TMH_x clathrates (LaH₁₀-type with H₃₂ cage, YH₉-type with H₂₉ cage, and CaH₆-type with H₂₄ cage),57 only YH₉-type clathrates can meet these requirements. Here, we assume that the smaller the number of TM valence electrons, the higher the average bonding unsaturation of the H cage. This is also confirmed by the calculated charge-transfer amount (Fig. S16).

As shown in Fig. 3(f) and Fig. S16, with increasing average bonding unsaturation, the T_c values of TMH₉ (M = Y, Th, and Pr) present an increasing trend.57–59 This is also similar to the T_c evolution of H₃S after incorporating a small number of Si/C/P atoms, i.e., doped H₃S with more electron-deficient atoms than the S atom can elevate the bonding unsaturation and the corresponding T_c values (e.g., T_c = 283 K at 230 GPa, 220 K at 225 GPa, and 212 K at 200 GPa for 4% Si, 7% C, and 6.25% P-doped H₃S with average bonding unsaturations of 0.257, 0.262, and 0.255, respectively, higher than H₃S).33,60,61 These results provide further support for the relationship between bonding unsaturation and superconductivity. It is noteworthy that our work opens the door to greater understanding of bonding-unsaturation-dependent superconductivity; further efforts are needed for a deeper understanding.

To explore the origin of superconductivity in these materials, we take P₈S as an example due to it having the highest structural symmetry and the widest stabilizing pressure range among the P-rich compounds. The projected electronic band structure of P₈S reveals obvious P–S hybridization across the whole Brillouin zone [Fig. 4(c)]. Moreover, several bands crossing E_F demonstrates steep dispersion along the Γ–A, A–H, and L–H directions, as well as a flat feature around the Κ point [Fig. 4(c)]. The calculated EPC parameter (λ) of P₈S is 0.47 at 100 GPa [Fig. 4(d)], and this is comparable to the values of 0.52 for C2/m S₃PH at 150 GPa,37 0.59 for Cmcm PH₃ at 100 GPa,31 and 0.42 for Pnma H₄S₃ at 140 GPa,62 and it is slightly lower than the value of 0.64 for Cmc2₁ SP₂H at 150 GPa.37 Based on the Eliashberg spectral function α²F(ω) and the phonon density of states (PHDOS), we can establish that the phonons derived from P and S atoms show evident coupling across the whole phonon frequency range due to their similar mass and strong covalent interactions; the phonons in the frequency range 10.0–18.9 THz make the main contribution (59%) to the total λ, and these correspond to two large S- and P-dominated PHDOS peaks [Fig. 4(d)]. In particular, above 14.5 THz, there appears to be a rapid increase in λ, resulting in a large contribution of 40%; this is also illustrated in the phonon dispersion curves [Fig. 4(e)]. Such a strong EPC is thus mainly associated with P-atom vibrations and local phonon softening at the A point. Based on the Allen–Dynes-modified McMillan equation,54 the T_c value of P₈S is estimated to be 4.15 K at 100 GPa with a typical Coulomb pseudopotential parameter of μ* = 0.1.63,64 Combining this with the electronic band structure, we can draw the conclusion that the superconductivity of P₈S principally comes from the coupling between the P/S 3p electrons and the P-/S-derived phonons.

To obtain insight into the pressure dependence of the superconductivity of P₈S, we investigated its superconductivity in the pressure range 55–125 GPa. The variation of T_c and the two key parameters (λ and ω_log) with pressure are displayed in Fig. 4(f). As the pressure decreases, T_c gradually increases, reaching a maximum of 6.0 K at 55 GPa. Below 100 GPa, the variation of λ is linearly related to T_c; this is in sharp contrast to ω_log, indicating that λ dominates the evolution of T_c. Above 100 GPa, the value of λ increases slightly, whereas ω_log decreases significantly, and this is responsible for the evolution of T_c. Overall, the pressure dependence of T_c is the result of competition between λ and ω_log.

In summary, we performed an extensive first-principles structural search of P–S systems under high pressure and identified three stable P-rich compounds (Cmcm P₅S, P-31m P₈S, and P-1 P₁₁S) and a new low-pressure PS₂ phase with P-3m1 symmetry, which acts as an intermediate phase between the R-3m and C2/m phases, completing the understanding of the phase-transition process of PS₂. All the P-rich P_xS (x > 1) compounds have a common covalent framework, in which P and S atoms form an octahedral coordination via head-to-head hybridization of P/S 3p orbitals and unsaturated covalent bonds. More interestingly, the variation in bonding unsaturation induced by different P contents is nearly linearly related to the NEF value, the EPC parameter λ, and the T_c value. A similar relationship is also observed in MH₉-type clathrates and Si/C/P-doped H₃S. Furthermore, the superconducting origin of this kind of compound is closely associated with the coupling of 3p electrons with the middle-/high-frequency phonons of P/S atoms. Our results indicate that modulating the bonding unsaturation in covalent frameworks could be a useful strategy for improving superconductivity.

Acknowledgment. The authors acknowledge funding from the Natural Science Foundation of China (Grant Nos. 21873017 and 21573037), the Natural Science Foundation of Hebei Province (Grant Nos. A2019203507 and B2021203030), the Postdoctoral Science Foundation of China (Grant No. 2013M541283), and the Natural Science Foundation of Jilin Province (Grant No. 20190201231JC). This work was carried out at the National Supercomputer Center in Tianjin, and the calculations were performed on TianHe-1 (A).