A First-Principles Exploration of Na x S y Binary Phases at 1 atm and Under Pressure

: Interest in Na-S compounds stems from their use in battery materials at 1 atm, as well as the potential for superconductivity under pressure. Evolutionary structure searches coupled with Density Functional Theory calculations were employed to predict stable and low-lying metastable phases of sodium poor and sodium rich sulﬁdes at 1 atm and within 100–200 GPa. At ambient pressures, four new stable or metastable phases with unbranched sulfur motifs were predicted: Na 2 S 3 with C 2/ c and Imm 2 symmetry, C 2-Na 2 S 5 and C 2-Na 2 S 8 . Van der Waals interactions were shown to affect the energy ordering of various polymorphs. At high pressure, several novel phases that contained a wide variety of zero-, one-, and two-dimensional sulfur motifs were predicted, and their electronic structures and bonding were analyzed. At 200 GPa, P 4/ mmm -Na 2 S 8 was predicted to become superconducting below 15.5 K, which is close to results previously obtained for the β -Po phase of elemental sulfur. The structures of the most stable M 3 S and M 4 S, M = Na, phases differed from those previously reported for compounds with M = H, Li, K. optB88-vdW, convex hull, lowest and C 2 phases 3 and 25 meV/atom These results suggest that other energetically competitive polymorphs, based upon unbranched S 2 − 5 units with either cis or trans geometries, could potentially be constructed, and that the computed energy ordering depends upon the method used to treat dispersion. PBE calculations showed that all three Na 2 S 5 polymorphs had indirect band gaps with values of 1.73, 1.47, and 1.30 eV for the (cid:101) , α , and C 2 phases, respectively (see the Supplementary Materials). Better estimates could be obtained using hybrid density functionals or GW, however the PBE results suggest that the conformation of the S 2 − 5 anion and the geometry of the cell both have an effect on the band gap.


Introduction
At atmospheric pressures, the size of the metal atom is thought to be important in determining the alkali metal polysulfide stoichiometries that are stable. In the case of the lightest alkali, lithium, only Li 2 S and Li 2 S 2 have been made [1], whereas for the heavier metal atoms (M = K, Rb, Cs) the known phases include M 2 S n with n = 1-6 [2][3][4][5]. The solid state Na-S system is particularly fascinating, with manuscripts reporting the failed synthesis of previously published stoichiometries, or the synthesis of new polymorphs. Na 2 S, Na 2 S 2 , Na 2 S 4 and Na 2 S 5 are stable or metastable at atmospheric conditions [2,[6][7][8]. Two polymorphs of Na 2 S 2 are known: the α form with three formula units per cell that is stable below 170 • C, and the higher temperature β form with two formula units per cell [9]. Several studies have failed to synthesize or isolate Na 2 S 3 , with the reaction products yielding a mixture of Na 2 S 2 and Na 2 S 4 instead [6]. A novel synthesis route in liquid ammonia yielded Na 2 S 3 , but the product decomposed around 100 • C [10]. The crystal structure of α-Na 2 S 5 has been solved [11], and several metastable polymorphs with this stoichiometry have been detected via Raman spectroscopy [2] and X-ray diffraction [8], but their structures were not determined. This amazing structural versatility is in part due to the ability of sulfur to form anionic polysulfide chains, S 2− n , with various lengths. Only unbranched chains are found in the sodium polysulfides, but they may have different arrangements or conformations.
Research on the Na-S system has been motivated by Kummer and Weber's development of the battery containing these two elements at the Ford Motor Company in 1966 [12][13][14]. The advantages of the Na-S battery include the fact that it is made from inexpensive materials, has a long cycle lifetime,

Computational Details
Evolutionary structure searches were carried out to find stable and low-lying metastable crystals in the Na-S phase diagram: in the sulfur rich region the Na 2 S n stoichiometry with n = 2-6, 8, and in the metal rich region the Na n S stoichiometry with n = 2-4 were considered. The calculations were carried out using the open-source evolutionary algorithm (EA) XTALOPT [38,39] version 10 [40], wherein duplicate structures were removed from the breeding pool via the XTALCOMP algorithm [41], and random symmetric structures were created in the first generation using RANDSPG [42]. EA searches were performed on structures containing 1-4 formula units in the primitive cell at 0, 100, 150 and 200 GPa. Minimum interatomic distance constraints were chosen to generate the starting structures in each generation. The minimum distances between Na-Na, Na-S and S-S atoms were set to 1.86, 1.45, and 1.04 Å, respectively. To improve the efficiency of the CSP searches, and increase the size of the unit cell that could be considered, the evolutionary search was seeded with experimentally determined and theoretically predicted structures from the literature, when available. The most stable structures found in the high pressure EA searches were optimized between 100 and 200 GPa, and their relative enthalpies and equations of states (EOS) are provided in the Supplementary Materials.
Geometry optimizations and electronic structure calculations including the densities of states (DOS), band structures, electron localization functions (ELFs) and Bader charges were carried out using DFT as implemented in the Vienna Ab-Initio Simulation Package (VASP) [43,44]. The bonding of select phases was further analyzed by calculating the crystal orbital Hamilton populations (COHP) and the negative of the COHP integrated to the Fermi level (-iCOHP) using the LOBSTER package [45]. At all pressures, the gradient-corrected exchange and correlation functional of Perdew-Burke-Ernzerhof (PBE) [46] was employed. It has been shown that it is important to include van der Waals (vdW) interactions to obtain reasonable estimates for the volume of α-S at 1 atm [19]. Therefore, the most stable structures from the 0 GPa PBE EA searches were reoptimized with the optB88-vdW functional [47,48]. In the EA searches, we employed plane-wave basis set cutoff of 325-400 eV, and the k-point meshes were produced by the Γ-centered Monkhorst-Pack scheme with the number of divisions along each reciprocal lattice vector chosen so that the product of this number with the real lattice constant was 30 Å. This value was increased to 50 Å for precise optimizations. The atomic potentials were described using the projector augmented wave (PAW) method [49]. The S 3s 2 3p 4 electrons were treated explicitly in all of the calculations. At 0 GPa the Na 3s 1 valence electron configuration was used, and at higher pressures an Na 2s 2 2p 6 3s 1 valence configuration was employed. For precise optimizations, the plane-wave basis set cutoffs were increased to 700 eV at 0 GPa, and 1000 eV at high pressures.
To verify the dynamic stability phonon band structures were calculated via the supercell approach [50,51], wherein the dynamical matrices were calculated using the PHONOPY code [52]. The electron-phonon coupling (EPC) calculations were performed using the Quantum Espresso (QE) program [53]. The Na (2s 2 2p 6 3s 1 ) and S (3s 2 3p 4 ) pseudopotentials, obtained from the PSlibrary [54], were generated by the method of Trouiller-Martins [55] with the PBE functional, and an energy cutoff of 90 Ry was chosen. The Brillouin zone sampling scheme of Methfessel-Paxton [56] and 24 × 24 × 6 k-point grid and a 8 × 8 × 2 q-point grid were used for P4/mmm Na 2 S 8 at 200 GPa. The EPC parameter, λ, was calculated using a set of Gaussian broadenings with an increment of 0.02 Ry from 0.0 to 0.600 Ry. The broadening for which λ was converged to within 0.05 Ry was 0.10 Ry. T c was estimated using the Allen-Dynes modified McMillan equation [57] with a renormalized Coulomb potential, µ * , of 0.1.

Stable and Metastable Na-S Phases at Atmospheric Conditions
The enthalpies of formation, ∆H F , of the most stable Na-S phases found in our EA searches are plotted in Figure 1 as a function of pressure. The phases whose ∆H F lie on the convex hull are thermodynamically stable, while those whose ∆H F are not too far from the hull may be metastable stable provided their phonon modes are all real. At atmospheric pressures, calculations carried out with both the PBE and optB88-vdW functionals showed that the Na 2 S, Na 2 S 2 , Na 2 S 4 and Na 2 S 5 stoichiometries lay on the hull. The ∆H F of Na 2 S 3 , Na 2 S 6 and Na 2 S 8 were slightly above the hull, whereas those of Na 3 S and Na 4 S were further away from it. The inclusion of vdW interactions lowered the total ∆H F by no more than 80 meV/atom, but the stoichiometries that lay on the hull were the same as those found within PBE. The inclusion of the zero point energy, ZPE, increased the ∆H F by no more than 9 meV/atom, but it also did not affect the identities of the thermodynamically stable phases. The structural parameters and ∆H F of the stable and important metastable phases are provided in the Supplementary Materials. All of the structures identified in the EA search that were within 15 meV/atom of the lowest enthalpy geometry for a particular stoichiometry were examined, and none of them contained sulfur anions with branched chains. , open symbols lie above it. ∆H F was calculated using the enthalpies of the experimentally known structures: body-centered cubic (bcc, 0 GPa), face-centered cubic (fcc, 100 GPa) [58], cI16 (150 GPa) [59,60], and hp4 (200 GPa) [61] for Na, and α-S (0 GPa) [62]. Because the crystal structure of S above 83 GPa is still disputed [63,64], the β-Po phase [65] was employed between 100 and 200 GPa, as in recent studies [35].
Momida [19] and Wang et al. [21] studied the effect of vdW interactions on the unit cell volumes and formation enthalpies of the sodium polysulfides. The effects of dispersion were approximated by using Grimme's semi-empirical DFT-D2 method [66] in conjunction with the PBE functional. Both studies found that vdW interactions resulted in slightly smaller volumes for the Na-S compounds, and more negative formation enthalpies. However, the inclusion of dispersion was crucial so that reasonable estimates of the unit cell volume of elemental α-S, which is composed of molecular S 8 rings, could be calculated [19]. Herein, the optB88-vdW method [47,48], which employs a non-local correlation functional that approximately accounts for dispersion interactions, was used. It has been demonstrated that this functional is among those that provides the best agreement with experiment for the volumes and lattice constants of layered electroactive materials for Li-ion batteries [67], as well as a broad range of metallic, covalent and ionic solids [48]. Other choices that might provide even more accurate lattice parameters at ambient pressures include combining PBE [46], and PBEsol [68] or its improvements [69] with rVV10 [70,71], or SCAN+rVV10 [72]. A comparison of the calculated volumes per atom obtained via PBE and optB88-vdW are provided in the Supplementary Materials. For elemental sulfur, PBE yields a cell volume of 37.21 Å 3 /atom, which is 36.4% larger than the experimental volume of 25.76 Å 3 /atom [73]. The optB88-vdW functional yields a volume of 24.90 Å 3 /atom (c.f. 26.88 Å 3 /atom with PBE-D2 [19]), which is only 3.4% lower than experiment. In Na-S systems that contained at least 50 mole % sodium, the difference between the PBE and optB88-vdW volumes was <4%, otherwise the difference ranged from 5-16%, depending on the stoichiometry and polymorph. For Na 2 S, α-Na 2 S 2 , β-Na 2 S 2 , Na 2 S 4 , and Na 2 S 5 optB88-vdW yielded volumes of 22.77, 21.51, 21.90, 21.85, and 22.38 Å 3 /atom, respectively, which differs by <4% from Momida's PBE-D2 results [19]. The 0 K optB88-vdW ∆H F values for the most stable Na 2 S, Na 2 S 2 , Na 2 S 4 , and Na 2 S 5 polymorphs were calculated to be −1.17, −0.94, −0.67, and −0.58 eV/atom, which is in good agreement with the experimental ∆H 0 F values at 298.15 K of −1.26, −1.03, −0.71, and −0.61 eV/atom, respectively [74].
The anti-CaF 2 Na 2 S structure with Fm3m symmetry was the lowest point on the 0 GPa convex hull [19][20][21], and the Na 2 S 2 stoichiometry had the second most negative ∆H F . Our EA searches were seeded only with the known α and β-Na 2 S 2 polymorphs, but they also readily identified a higher energy γ-Na 2 S 2 phase that has recently been predicted [21]. All of these three polymorphs are comprised of Na + cations and S − 2 anions. In agreement with previous DFT calculations [19][20][21], the β polymorph was computed to have the lowest ∆H F , followed by the α, and the γ configurations. The PBE/optB88-vdW differences in energy between the α and β structures were comparable to the difference between the β and γ structures, 4/8 meV/atom and 5/9 meV/atom, respectively.
In addition to Na 2 S and Na 2 S 2 , the Na 2 S 4 and Na 2 S 5 stoichiometries also lay on the PBE and optB88-vdW convex hulls. The EA search was seeded with the known I42d-Na 2 S 4 structure [75], which contains an unbranched S 2− 4 chain whose S-S bond angle and dihedral angle were computed to be 111.3 • and 96.7 • , within PBE, respectively. No other polymorphs with comparable energies were found.
The EA search was also seeded with the α-Na 2 S 5 [11] and -Na 2 S 5 [20] polymorphs illustrated in Figure 2a,b. In the α form, the unbranched S 2− 5 anion adopts a bent (cis) configuration, whereas, in the form, it is stretched (trans), as in K 2 S 5 [76], Rb 2 S 5 [77] and Cs 2 S 5 [78]. PBE and PBE-D2 calculations have shown that neither the α nor the phases lay on the convex hull [19,20], but CSP has found another currently unsynthesized phase with stretched S 2− 5 anions that was thermodynamically stable [20]. The coordinates of this phase were not provided in Ref. [20], however it appears to be different from the lowest enthalpy C2 symmetry Na 2 S 5 phase from our EA searches shown in Figure 2c. In the structure of Mali and co-workers [20], neighboring S 2− 5 chains point in opposite directions (similar to what is observed in the phase along the b-lattice vector), whereas in C2-Na 2 S 5 they face the same direction. C2-Na 2 S 5 lay on the PBE convex hull, and it's enthalpy was 4 and 5 meV/atom lower than the and α polymorphs, respectively. The order of stability was not affected by the ZPE contributions. On the other hand, within optB88-vdW, the phase, which lay on the convex hull, had the lowest enthalpy with the α and C2 phases being 3 and 25 meV/atom higher, respectively. These results suggest that other energetically competitive polymorphs, based upon unbranched S 2− 5 units with either cis or trans geometries, could potentially be constructed, and that the computed energy ordering depends upon the method used to treat dispersion. PBE calculations showed that all three Na 2 S 5 polymorphs had indirect band gaps with values of 1.73, 1.47, and 1.30 eV for the , α, and C2 phases, respectively (see the Supplementary Materials). Better estimates could be obtained using hybrid density functionals or GW, however the PBE results suggest that the conformation of the S 2− 5 anion and the geometry of the cell both have an effect on the band gap.  [11] and (b) -Na 2 S 5 [20] polymorphs, as well as the newly predicted (c) C2-Na 2 S 5 structure. Sodium atoms are colored blue, and sulfur atoms are yellow.
Besides the previously mentioned thermodynamically stable stoichiometries, several sulfur-rich containing phases, Na 2 S n (n = 3, 6, 8), were found to be low-energy metastable species, as confirmed by phonon calculations. For the Na 2 S 3 stoichiometry, the Cmme, C2/c-I, C2/c-II, and Imm2 polymorphs illustrated in Figure 3 lay 7/22, 3/1, 13/9, and 17/40 meV/atom above the convex hull within PBE/optB88-vdW calculations. CSP previously predicted the Cmme [20] and C2/c-I [21] structures, used as seeds in our EA searches, whereas the C2/c-II and Imm2 phases discovered here have not been reported before. All four of these polymorphs contained V-shaped S 2− 3 anions with S-S bond lengths of 2.087-2.107 Å and bond angles of 106.00-111.21 • . The main difference between them was the relative arrangement of the S 2− 3 motifs. In Imm2 all of the V-shaped anions pointed in the same direction, whereas in C2/c-II those in a single layer pointed in the same direction, but those in the adjacent layer were rotated by 180 • . In Cmme and C2/c-I adjacent rows of anions in the same layer pointed in opposite directions. In Cmme the apex of the Vs in one layer were located directly behind those in an adjacent layer, but rotated by 180 • . In C2/c-I, the Vs in adjacent layers also faced opposite directions. All four polymorphs had indirect band gaps, with the PBE value for C2/c-I, 1.06 eV, being about 0.5 eV smaller than for the other three. The closeness of the ∆H F of these four polymorphs to the convex hull, and their dynamic stability suggests that they may be synthesizable. Experimentalists have not yet been able to synthesize a persistent Na 2 S 3 compound, yielding a mixture of Na 2 S 4 and Na 2 S 5 either directly [6,20,79] or after disproportionation near 100 • C [10], suggesting that the kinetic barriers towards decomposition may be low.
Seed structures were not employed in EA searches carried out on the Na 2 S 6 stoichiometry, and the two nearly isoenthalpic C2 and P1 symmetry structures illustrated in Figure 4a,b were found. Whereas the former was comprised of unbranched S 2− 6 units, the latter contained V-shaped trimers. These two dynamically stable phases lay only 2/3 and 5/16 meV/atom above the PBE/optB88-vdW convex hulls, respectively. The inclusion of the ZPE decreased the PBE difference in energy to 1 meV/atom. Mali and co-workers found the same two structures in their CSP searches [20]. The calculated phonon spectrum of P1-Na 2 S 6 contained two bands above 500 cm −1 (565 and 543 cm −1 at Γ), which is higher than in any of the other 0 GPa phase considered herein. The S-S bond distances in P1-Na 2 S 6 were shorter than in any of the predicted Na 2 S 3 phases, 1.972 and 2.048 Å, giving rise to the increased frequency. The decreased distance is likely a result of the smaller formal charge on the trimer, S − 3 in P1-Na 2 S 6 vs. S 2− 3 in the Na 2 S 3 polymorphs in Figure 3. Molecular calculations using the ADF program package [80,81] with a triple-zeta basis set with polarization functions (TZP), and the PBE functional yielded S-S bond lengths of 2.037 and 2.166 Å in the monoanion and dianion, respectively, which is in reasonable agreement with the periodic calculations. Both phases were computed to have indirect PBE band gaps, 1.75 eV in the C2 and 0.72 eV in the P1 phases, respectively.  [20], and (b) C2/c-I Na 2 S 3 [21] phases, as well as the newly found (c) C2/c-II Na 2 S 3 , and (d) Imm2-Na 2 S 3 geometries. Sodium atoms are colored blue, whereas sulfur atoms in adjacent layers are yellow and orange. Finally, EA searches were carried out on the Na 2 S 8 stoichiometry, which has not been considered in other studies before. The dynamically stable C2 symmetry structure shown in Figure 4c was found to have the lowest energy and it lay 10/30 meV/atom above the PBE/optB88-vdW convex hulls. This phase is comprised of an unbranched S 2− 8 stretched, helical-like chain. Its PBE band gap was indirect and measured 1.56 eV.
Two sodium rich metastable phases, whose structures can be derived from the known 1 atm Fm3m-Na 2 S configuration, which has one formula unit per primitive cell and is shown in Figure 5a, were found. The P3m1-Na 3 S and Pmn2 1 -Na 4 S phases illustrated in Figure 5b,c lay 69/71 and 67/70 meV/atom above the PBE/optB88-vdW convex hull, respectively. In the antifluorite Na 2 S structure the S 2− anions are arranged on a face centered cubic (fcc) lattice, and all of the tetrahedral sites are filled with Na + cations, with the Na-Na nearest neighbor distances measuring 3.267 Å. The P3m1-Na 3 S phase can be derived from a three-formula unit supercell of the antifluorite structure where one sulfur atom has been removed from the fcc lattice. This results in a distortion from perfect cubic symmetry, with Na-Na-Na angles measuring 89.3 and 89.9 • , and Na-Na nearest neighbor distances measuring 3.288 and 3.242 Å. Another dynamically stable R3m symmetry phase that was 0.7/1.5 meV/atom less stable within PBE/optB88-vdW was also found in the structure searches. Its coordinates and computed properties are provided in the Supplementary Materials. The Pmn2 1 -Na 4 S structure contained four formula units. It can be obtained by inserting layers of sodium atoms into a supercell of the antifluorite structure, which results in a minimal structural distortion in the anti-CaF 2 layers, with Na-Na distances along the a and c lattice vectors measuring 3.344 and 3.267 Å, respectively. The PBE band gap in Fm3m-Na 2 S was computed to be 2.47 eV, as expected for an ionic phase. As shown in the Supplementary Materials, the sodium rich P3m1-Na 3 S and Pmn2 1 -Na 4 S compounds are metals. Plots of the ELF, given in Figure 6, illustrate that regions where the bonding is metallic, indicated by an ELF value of ∼0.5, are localized along the two-dimensional layers of sodium atoms. The DOS of both structures between −2.5 eV to the Fermi level has step-like features, as would be expected for a two-dimensional electron gas [82]. The ELF plots also show spherical regions with high ELF values that encompass the sulfur atoms, suggestive of an S 2− oxidation state, and regions with a high ELF value between adjacent Na layers within Pmn2 1 -Na 4 S . The latter are labeled "Es" in Figure 6d, since they resemble interstitial electrons that are paired in anion-like species, which have been computed for a number of high-pressure electrides exhibiting ionic bonding [83,84].

Na-S System at High Pressure
Because dispersion interactions are unlikely to be important under pressure [85], calculations at 100, 150 and 200 GPa were carried out with the PBE functional. Figure 1 illustrates that only Na 3 S and Na 2 S were thermodynamically stable at all of these pressures. Na 2 S 2 lay on the 100 and 150 GPa convex hulls, whereas Na 4 S and Na 2 S 4 lay on the 200 GPa convex hull. The inclusion of the ZPE corrections did not affect the identity of the thermodynamically stable stoichiometries at 100 and 150 GPa, but it added Na 2 S 2 and Na 2 S 3 to the 200 GPa hull, and pushed Na 2 S 4 away from it. Because phonon calculations (see the Supplementary Materials) confirmed that all of these phases were dynamically stable, and high-pressure experiments can lead to the formation of metastable species [86]; the structures of all of these phases are discussed below.
The Na 2 S stoichiometry was the lowest point on the convex hull among all other binary phases between 100-200 GPa. The species seeded into the EA, including the experimentally known Ni 2 In structure (P6 3 /mmc) [26] and a previously predicted anti-AlB 2 structure (P6/mmm) [31] were found to be the most stable at 100-150 GPa and 200 GPa, respectively. DFT calculations have predicted that Na 2 S is an insulator up to 300 GPa, with PBE band gaps of ∼3.50, 3.25, and 1.70 eV at 100, 150 and 200 GPa, respectively [31].
The evolutionary search for Na 2 S 2 at 100-200 GPa found a Pbam symmetry structure, shown in Figure 7a, which contained two formula units in its primitive cell. It lay on the convex hull at 100 and 150 GPa, and was only 5 meV/atom away from the hull at 200 GPa. Its shortest S-S distances were no longer than 2.040 Å at all of the pressures studied, which is comparable to the measured S-S bond length within α-S at ambient conditions, 2.055 Å [73]. A plot of the ELF, shown in the Supplementary Materials, confirmed the presence of layers of S 2− 2 dimers, arranged in a herringbone type fashion. The Na 2 S 3 stoichiometry lay above the convex hull at all pressures considered, however the distance from the hull decreased from 65 to 25 meV/atom from 100 to 200 GPa. The C2/c phase shown in Figure 7b, which was metastable between 100 and 150 GPa, is comprised of one-dimensional branched tertiary sulfur chains with S-S distances of 2.180 Å along the chain, and 2.030 Å along the branch. By 150 GPa, the I4/mmm structure illustrated in Figure 8a, was enthalpically preferred. It is comprised of two layers of square sulfur nets, with S-S distances of 2.092 Å at 200 GPa, separated by sodium nets. Na 2 S 4 was thermodynamically stable at 200 GPa, but it lay 55 and 3 meV/atom from the 100 and 150 GPa hulls. The most stable structure between 100-110 GPa, shown in Figure 7c, had C2/c symmetry and contained layers of hexagonal puckered corner-sharing chair cyclohexane-like rings with S-S distances of 2.105 and 2.232 Å. The Bader charges on the 4-and 2-coordinated sulfur atoms were −0.16e and −0.60e, respectively. Above 110 GPa, the I4/mmm structure, shown in Figure 8b, was preferred. Similar to I4/mmm-Na 2 S 3 , it was comprised of two sets of two-dimensional square sulfur layers with S-S distances of 2.092 Å at 200 GPa separated by a single square net of sodium atoms.
The stoichiometries with an even larger ratio of sulfur to sodium lay between 35-85 meV/atom from the non-ZPE corrected convex hull. The structure dubbed C2/m-I Na 2 S 5 was the most stable at 100 GPa, but by 130 GPa another phase with the same symmetry, C2/m-II Na 2 S 5 , became preferred (see Figure 7d,e). At 100 GPa, the sulfur sublattice in the former contained zigzag chains with S-S distances of 2.076 Å, as well as one-dimensional corner sharing squares with S-S distances of 2.194 Å running along the b-lattice vector. The higher pressure structure was comprised of these same one-dimensional square nets fused at the edges with S-S distances ranging from 2.029 to 2.070 Å at 200 GPa.
At 100 GPa, the lowest enthalpy Na 2 S 6 structure adopted C2 symmetry (see Figure 7f), and it lay 60 meV/atom above the convex hull. It was comprised of sulfur atoms as well as dimers and V-shaped trimers, as shown in the ELF plots in the Supplementary Materials, with S-S bond lengths of 2.062 Å, and 1.998-2.180 Å, respectively. The I4/mmm phase illustrated in Figure 8c lay 35 and 45 meV/atom above the hulls at 150 and 200 GPa, respectively. It resembled the high-pressure Na 2 S 4 and Na 2 S 3 phases, except that its sulfur layers contained both edge-sharing cubes and square nets with S-S distances of 1.916-2.111 Å and 2.111 Å, respectively at 200 GPa. Evolutionary searches carried out at 100, 150, and 200 GPa for Na 2 S 8 discovered the metastable P4/mmm phase shown in Figure 8d. This species was reminiscent of I4/mmm-Na 2 S 6 with layers comprised of fused cubes and square nets with S-S distances of 1.903-2.124 Å and 2.124 Å, respectively, at 200 GPa. ELF plots for P4/mmm-Na 2 S 8 are provided in the Supplementary Materials because they are representative of the ELFs obtained for phases that contained square and/or cube-like sulfur nets. In all of the structures studied, the ELF plots suggested that the S-S bonds within the two-dimensional sheets were weaker than the bonds between atoms in different sheets. For example, in P4/mmm-Na 2 S 8 the S-S distances along the aand b-directions were 2.124 Å and the negative of the crystal orbital Hamilton populations integrated to the Fermi level (-iCOHPs) were 2.5 eV. For the S-S bonds that were oriented along the c-direction these values became 1.903 Å and 5.5 eV, which is comparable with a bond length of 2.168 Å and -iCOHP of 4.8 eV in the S 2− 2 dimer in α-Na 2 S 2 at 0 GPa. The pressure induced polymerization and bond weakening observed in the two-dimensional sulfur layers has been seen in other systems before. For example, it has been shown that the Cl 2 molecules present in XeCl at 40 GPa undergo polymerization to form one-dimensional zigzag chains by 60 GPa with the -iCOHPs between nearest neighbor atoms computed to be 3.81 eV and 2.16 eV, respectively [85]. By 100 GPa, however, there is no evidence of Cl-Cl bond formation with computed ELF and -iCOHPs characteristic for monoatomic Cl.
Within the 100-200 GPa pressure range, the most stable Na 3 S stoichiometry unearthed contained two formula units in its primitive cell, and it possessed P6 3 /mmc symmetry. This structure, illustrated in Figure 9a, consisted of an ABAB. . . stacked triangular net where each sulfur atom was [6+6] coordinated with Na-S distances of 2.216 and 2.286 Å at 200 GPa. We also optimized Na 3 S stoichiometries that were analogous to the R3m and Im3m structures proposed for superconducting H 3 S [37]; as shown in the Supplementary Materials, their enthalpies were at least 600 meV/atom and 1000 meV/atom higher, respectively, between 100-200 GPa. Moreover, sodium analogs of the previously predicted Pm3m, I4/mmm and Fm3m-Li 3 S [35] and Pm3m, I4/mmm-K 3 S [36] phases were relaxed, and they were found to be at least 160 meV/atom less stable than P6 3 /mmc-Na 3 S between 100 and 200 GPa. Figure 8. Predicted crystal structures of Na 2 S n (n = 3, 4, 6, 8) phases at 200 GPa whose sulfur sublattices were comprised of two-dimensional square or cube nets, and whose sodium sublattices were comprised of square nets: (a) I4/mmm-Na 2 S 3 ; (b) I4/mmm-Na 2 S 4 ; (c) I4/mmm-Na 2 S 6 ; and (d) P4/mmm-Na 2 S 8 . Sodium atoms are colored blue, and sulfur atoms are yellow. The Cmcm symmetry Na 4 S phase, shown in Figure 9b, was also found to be the most stable in the whole pressure range of 100-200 GPa. At 200 GPa, it lay on the convex hull, but at 100 and 150 GPa, it was 35 and 10 meV/atom above the hull. It can be viewed as a cut triangular two-dimensional lattice that is joined to another such lattice via a square net. Individual sheets are stacked in an ABAB fashion. Geometry optimizations of the sodium analog of Cmcm-K 4 S [36] were between 25 and 50 meV/atom less stable than Cmcm-Na 4 S, whereas the enthalpies of the sodium analogs of R3m-Li 4 S [35] and K 4 S [36] were at least 200 meV/atom higher than the most stable structure found here.

Electronic Structure and Superconductivity under Pressure
The PBE band structures and electronic DOS plots of the high-pressure dynamically stable phases, provided in the Supplementary Materials, showed that within their range of stability they were all metallic. Under pressure elemental sodium transforms from the ambient pressure bcc structure to an fcc phase at 65 GPa [58], followed by the semi-metallic cI16 structure at 103 GPa [87]. Near the minimum of the melting temperature a number of complex crystal structures, with unit cells as large as 512 atoms, have been identified experimentally within a narrow temperature/pressure range [59]. Remarkably, as first predicted by theory [88], sodium undergoes a metal to insulator (MIT) transition by 200 GPa, which results from the overlap of core electrons and the concomitant localization of p − d hybridized valence electrons, which render this phase an electride [61]. Theoretical calculations have found that the maximum electron-phonon-coupling parameter, λ, for sodium occurs for the cI16 phase near 140 GPa [89]. Estimates of the T c within the Allen-Dynes modified McMillan equation and µ * = 0.13 found a maximum T c of 1.2 K, leading to the conclusion that superconductivity in sodium prior to the MIT is weak or non-existent. In addition, because DFT calculations showed that most of the metal rich binary lithium and potassium sulfides had no or low T c [35,36], we thought it unlikely that the P6 3 /mmc-Na 3 S and Cmcm-Na 4 S phases found herein would be good superconductors.
In contrast, compressed sulfur is among the elemental phases with the highest T c , reaching up to 17 K in the rhombohedral β-Po geometry near 160 GPa [90]. The β-Po structure can be viewed as a simple cubic (sc) lattice compressed along the body diagonal, and both can be described using the same unit cell with rhombohedral angles of 90 • (sc) and 104 • (β-Po) [91]. This made us wonder if the metastable I4/mmm-Na 2 S 6 or P4/mmm-Na 2 S 8 phases, which both contained sc-like sulfur layers, could potentially be candidates for conventional superconductivity? Calculations were carried out on the latter since it had a higher sulfur content, and a larger normalized density of states at the Fermi level. Figure 10 shows the computed phonon band structure, Eliashberg spectral function, and the electron-phonon integral for P4/mmm-Na 2 S 8 . Because of the similar masses of the sodium and sulfur atoms, both atom types contribute to the phonon modes across the frequency spectrum. The EPC was calculated to be λ = 0.79 and the logarithmic average frequency ω log = 344.8 K, resulting in a T c of 15.5 K using a µ * = 0.1. For comparison, LDA calculations on the β-Po structure at 200 GPa obtained λ = 0.78 and T c = 19.2 K using the Allen-Dynes approximation with µ * = 0.11 [91], and T c ∼16 K via the multiband approach within the superconducting density functional theory (SCDFT) formalism [92]. Thus, the propensity for superconductivity in P4/mmm-Na 2 S 8 under pressure appears to be similar to that of elemental sulfur.

Conclusions
Evolutionary structure searches coupled with first-principles calculations have been employed to explore the phase diagram of sodium-rich and sodium-poor sulfides at 1 atm and within 100-200 GPa. At ambient pressure, the Na 2 S n , n = 1, 2, 4, 5, stoichiometries were thermodynamically stable, whereas n = 3, 6, 8 were low lying metastable species. In addition to identifying experimentally known or previously predicted species, we also found the novel C2/c-II and Imm2-Na 2 S 3 phases, which differed in the orientation of the S 2− 3 anions, as well as C2-Na 2 S 5 , which contained unbranched trans S 2− 5 chains, and C2-Na 2 S 8 , which was comprised of S 2− 8 helical-like chains. Even though the inclusion of van der Waals interactions did not have much of an impact on the optimized volumes of these phases, in some cases it did have an effect on the computed energy ordering of different polymorphs. Two dynamically stable sodium rich phases, P3m1-Na 3 S and Pmn2 1 -Na 4 S, whose enthalpies of formation lay ∼70 meV/atom above the convex hull were also identified. Their structures were related to Fm3m-Na 2 S, and they had intriguing electronic structures whose metallicity was derived from two-dimensional sodium sheets.
With the exception of the Na 2 S stoichiometry, the high-pressure Na-S phase diagram has not been studied before. The stable and metastable sulfur-rich phases contained sulfur motifs that included atoms, dimers, V-shaped trimers, branched polymeric and zigzag chains, as well as two-dimensional corner-linked cyclohexane-like rings, and edge/corner-linked squares. At 200 GPa, the most stable Na 2 S n , n = 3, 4, 6, 8, phases were comprised of sodium nets and two-dimensional fused cube or square sulfur nets, whose electronic structure was interrogated. The Allen-Dynes modified McMillan formula was used to predict the superconducting critical temperature, T c , of Na 2 S 8 at 200 GPa, and it was found to be comparable to previous theoretical estimates for the β-Po sulfur phase, 15.5 K for the polysulfide vs. 16-19 K for elemental sulfur [91,92]. At 150-200 GPa, P6 3 /mmc-Na 3 S and Cmcm-Na 4 S were found to be thermodynamically stable, and their structures differed from those proposed earlier for H 3

Conflicts of Interest:
The authors declare no conflict of interest.