Superconductivity and strong anharmonicity in novel Nb-S phases

In this work we explore the phase diagram of the binary Nb-S system from ambient pressures up to 250 GPa using $ab$ $initio$ evolutionary crystal structure prediction. We find several new stable compositions and phases, especially in the high-pressure regime, and investigate their electronic, vibrational, and superconducting properties. Our calculations show that all materials, besides the low-pressure phases of pure sulfur, are metals with low electron-phonon coupling strengths and critical superconducting temperatures below 15 K. Furthermore, we investigate the effects of phonon anharmonicity on lattice dynamics, electron-phonon interactions, and superconductivity for the novel high-pressure phase of $\mathrm{Nb}_2\mathrm{S}$, demonstrating that the inclusion of anharmonicity stabilizes the lattice and enhances the electron-phonon interaction.


I. INTRODUCTION
In recent years, transition metal chalcogenides (TMC) have gained significant research interest based on their chemical and physical variety, as well as their tendency to create layered structures [1,2]. This constraint in geometry has been shown to be the source of many interesting phenomena related to electron-phonon (ep) processes such as charge density waves (CDW) and superconductivity (SC) [3][4][5][6]. For some TMC materials, Mottinsulating behaviour or potential topological SC has also been predicted [7][8][9], properties that could find future applications in high-performance electronics, as controllable Mott transitions allow very fast, precise and efficient field-effect transistors (FET) [10], and Majorana states in topological SC could smooth the way to fault-tolerant quantum computing [11,12].
Extensive studies have been conducted specifically on transition metal dichalcogenides (TMD). Featuring weakly bound van-der-Waals (vdW) layers, they offer a plethora of ways to manipulate and influence their physical properties, for example by varying the stacking order of the layers, intercalating other elements, doping via ionic liquids, etc. [13][14][15]. This leads to applications as cathodes for batteries [16], solid lubricants due to low friction between the layers [17,18], electrocatalysts for energy conversion [19,20], and makes them versatile candidates for 2D beyond-graphene materials [21][22][23].
The semiconducting phases with group-4 and group-6 transition metals are studied for their applications in photovoltaics and electronics [24][25][26]. Many TMDs undergo an indirect-to-direct band gap transition with decent carrier mobilities when bulk materials are exfoliated down to the monolayer limit [27][28][29][30], making them favourable as FET components in terms of efficiency as well as size [31][32][33][34].
Group-5 TMDs, by contrast, can be metallic by filling the lowest d band and exhibit conventional ep-mediated SC and CDW order. The coexistence and interaction of * christoph.heil@tugraz.at these two, at first glance mutually exclusive phases, is the focus of intensive research at the moment. Recent studies suggest that the ep coupling responsible for creating the CDW phase is strongly localized in k-space, allowing for the remaining ep interactions at other wave vectors to promote SC [35][36][37][38][39][40][41][42][43].
External pressure can have significant effects on many physical properties, as it strongly influences atomic bonding and hence electronic properties, lattice dynamics, and ep interactions, driving materials towards or away from lattice instabilities [44][45][46][47][48][49][50]. While many TMD phases are well studied up to higher pressures, including investigations of their electronic and vibrational properties, and phase diagrams of other specific TMC stoichiometries at lower pressures have been reported [2,3,[51][52][53], comprehensive and systematic studies on possible stable phases up to high pressures are still missing to the best of our knowledge.
With this work we want to contribute to fill this void by analysing the high-pressure phase diagram of Nb-S. We find several novel materials and phases, for which we determine their thermodynamic stabilities, as well as their electronic and vibrational behaviours, placing particular focus on ep interactions and associated physical phenomena such as SC.

II. METHODS AND COMPUTATIONAL DETAILS
We investigated the phase space using the USPEX package for evolutionary crystal structure prediction [54,55] in the pressure range from 0 to 250 GPa. For an overview, we started variable-composition runs in steps of 50 GPa, providing already known phases as additional seed structures, and further refined our search to find phase boundaries within an 5 GPa accuracy [56]. We employed density functional theory (DFT) via the Quantum ESPRESSO (QE) package [57] for all thermodynamically stable and metastable structures to relax the unit cells below a threshold of 10 −6 Ry/a 0 for all force components, and to calculate the electronic band structure and density of states (DOS) with an accuracy in the total energy be- Convex hulls of enthalpy of formation ∆H for compounds NbxSy over sulfur fraction y/(x + y) for different pressures. The circles mark stable structures, the crosses indicate a few selected metastable structures that are less than 60 meV above the convex hull. For Nb 2 S, we additionally employed maximally localized Wannier functions and the fully anisotropic Migdal-Eliashberg theory as implemented in the EPW package to calculate λ q,ν , nesting functions ζ q , and superconducting gaps ∆ k (T ) as a function of temperature T [68, 69].

III. PHASE DIAGRAM
Following the evolutionary crystal searches we performed further relaxation and enthalpy calculations on the best candidates and determined convex hulls in the pressure range from ambient to 250 GPa pressure, as shown in Fig. 1 for five different pressures. A full phase diagram of stable structures and their space groups is presented in Fig. 2, where we indicate in blue previously reported phases and novel phases in green.
In the studied pressure range, we reproduce the known phases of the Nb-S system in excellent agreement with literature [39,53,54,[70][71][72][73][74][75][76][77][78][79][80][81][82]: At ambient pressure the convex hull shows the stability of previously reported phases Nb 14 S 5 , Nb 3 S 4 , (2H-)NbS 2 , and NbS 3 , while Nb 21 S 8 and NbS observed in higher-temperature experi- Phase diagram for the stable crystal structures in the Nb-S system. The numbers above the bars indicate the space group, the horizontal fading the transition pressures, and the colours previously reported (blue) and novel phases (green). The hatched region of Nb2S denotes the dynamically unstable pressure range including anharmonic corrections. ments [83-86] appear as metastable, 4 meV and 29 meV above the convex hull tie-line, respectively. Increasing pressure, the formation of NbS becomes favourable and it represents the convex hulls' minimum above 30 GPa with a maximum enthalpy difference of about 1.8 eV/atom at 250 GPa.
Our structure searches also reveal several novel phases and materials of the Nb-S system, for which we provide all crystallographic details and schematics in the Supplemental Material (SM) [87]: i) NbS in the high-pressure CsCl-type phase (space group 221), ii) Nb 2 S 3 in different low-and high-pressure phases (space groups 160 and 139), appearing as "doubled" 3R-NbS 2 and high-pressure NbS 2 , respectively (see Fig. S2 [87]), iii) Nb 3 S 5 , which appears as 3R-NbS 2 intercalated with Nb atoms in octahedral coordination, resulting in the low-symmetry space group 1, iv) Nb 3 S, a simple cubic material (space group 223) where the S atoms are in icosahedral coordination with the surrounding 12 Nb atoms, and the icosahedra being stacked in a bcc superlattice fashion, v) three phases of Nb 2 S (space groups 42, 11, and 189) that are described in detail in Sec. IV.
Based on the types of appearing crystal structures, the phase diagram can be be separated roughly into low-and high-pressure regions between 50 and 100 GPa. The low-pressure phases are predominantly built by layered NbS units in the trigonal prismatic or octahedral geometry known from transition metal dichalcogenides [88, 89], resulting in many layered (meta)stable structures with combined stacking patterns known from 1T-, 1H-, 2H-and 3R-NbS 2 (see Figs. S1 and S2 [87]). On the other hand, the high-pressure sulfur-rich phases, i.e. y/(x + y) ≥ 0.5 with x being the number of Nb atoms per formula unit and y the number of S atoms , bringing certain Nb atoms in a distance to S atoms (dashed lines) that is only 10% longer than the drawn bonds. (c) High-pressure Nb2S (space group 189) with S atoms in tricapped trigonal prismatic coordination with Nb. The coordination of a single S atom is highlighted in the lower 3D view, as a guide to the eye for the complex coordination geometry, a full (yellow) and an empty (grey) trigonal prism are shown in the upper right 3D view. A 2 × 2 × 1 super cell is shown from z-direction on the left, where again the two prisms are indicated to highlight their hexagram orientation. They are shown for two S atoms inside the unit cell and one S atoms at the unit cell border, which are conceptually equal, but differ slightly in bond lengths. (d) High-pressure Nb2S3 with a 2:2 stacking (see SM [87]). A 2 × 1 × 2 super cell in y-direction is shown on top, the stacking and cubic coordination are emphasized in the bottom panels.
per formula unit, consist of different layered stackings of the CsCl-like simple cubic NbS unit (space group 221) as illustrated in Fig. 3d.
For the high-pressure stacking with an even ratio n : n, i.e. n simple cubic NbS units on top of each other and n units shifted by (a/2, b/2, 0) shown in Fig. 3d and Fig. S3 [87], we find the stoichiometric formula Nb n S n+1 and the tetragonal space group 139, while an odd ratio n : (n + 1) results in Nb 2n+1 S 2n+3 and space group 123. Stable structures are predominantly formed with an even ratio and n = 1, 2, namely NbS 2 and Nb 2 S 3 , respectively, and metastable structures can be found with even and odd ratios up to n = 5 ( Fig. S3 [87]).
As widely reported in literature, the formation of layers is without doubt a common feature of transition metal chalcogenides and we therefore want to emphasize on some highly three-dimensional structures here: i) Nb 3 S 4 [73,90], in a hexagonal crystal structure, where all Nb atoms are in octahedral coordination, forming a 3D mesh-like structure. ii) The so far unknown Nb 3 S, being the only Nb-rich member of the phase diagram at pressures above 235 GPa. We also find an isostructural material with inverted empirical formula, NbS 3 , that is metastable at high pressures. iii) The so far unknown high-pressure structure of Nb 2 S, where nine Nb atoms forming a 14-faced structure around the S atoms are three-dimensionally packed in a hexagonal lattice. The crystal structure of Nb 2 S, which we will discuss in detail in Sec. IV, as well as examples of the characteristic crystal geometries in the Nb-S system are shown in Fig. 3.
For all materials not extensively discussed in literature we carried out electronic, vibrational, and ep calculations, revealing that all phases except pure sulfur at low pressures are metallic and superconducting with T c < 15 K. To not deter from the main focus of this paper, we provide an extensive discussion of all materials and their properties in the SM [87] and in the following analyse in detail the novel Nb 2 S compound, for which we find imaginary phonon modes within DFPT at certain parts of the Brillouin zone (BZ) that are indicative of lattice instabilities. In the next section we will show that this putative instability is due to the disregard of anharmonic effects and that considering the full anharmonic potential, Nb 2 S is in fact stable over a wider pressure range exhibiting superconductivity and strong ep coupling.

IV. Nb2S -ANHARMONICITY AND SUPERCONDUCTIVITY
At low pressures, we find two metastable phases of Nb 2 S (see convex hulls in Fig. 1), namely an orthorhombic structure (space group 42) at ambient pressure that consists of vdW-coupled layers in trigonal prismatic geometry around sulfur atoms [91], shown in Fig. 3a, and above 20 GPa a monoclinic structure (space group 11) with a different trigonal prismatic arrangement exhibiting a denser clustering of Nb atoms around S atoms (indicated by dashed lines in Fig. 3b). The former lies 6 meV above the convex hull tie-line, and the latter 60 meV. Above 40 GPa, a hexagonal phase (space group 189) with nine atoms per unit cell becomes more favourable, but appears on the convex hull only above 60 GPa due to the presence of Nb 14 S 5 . Here, the main repetition unit is a non-uniform 14-faced triaugmented triangular prism around S atoms, formed by two shifted and rotated triangular prisms in hexagram geometry as illustrated in Fig. 3c. As indicated by the shaded areas in the supercell, the side lengths of the two auxiliary trigonal prisms (yellow and grey) are different, as well as the side lengths of the conceptually equal prisms around S atoms inside the unit cell and around S atoms at the unit cell border. Above 235 GPa the simple cubic Nb 3 S structure becomes thermodynamically more favourable and Nb 2 S leaves the convex hull.
In Fig. 4, we show the electronic band structure around the Fermi level of the high-pressure Nb 2 S phase along a high-symmetry path of the hexagonal BZ. We find highly anisotropic dispersions in all directions, supporting our structural observations that this is indeed a highly threedimensional material. The DOS together with the partial projections onto the Nb d and S p orbitals are reported in the same figure. We find the Fermi level on the lower shoulder of a peak in the DOS, a property that remains robust over the stable pressure range, as demonstrated in Fig. S4 [87]. In the energy range around and especially at the Fermi energy, the main contributions to the total DOS originate from the Nb d orbitals, a prominent feature already observed in group-5 TMDs [26] and also present in most other phases in the Nb-S system [87].
Vibrational properties of high-pressure Nb 2 S are shown in Fig. 5. The main contribution to the phonon DOS up to the peak below 60 meV stems from the heavier Nb atoms, and above that from the lighter S atoms, as indicated by the white and grey areas underneath the DOS curve in Fig. 5b. The S-dominated DOS range is separated by two phonon band gaps around 70 and 85 meV. In the harmonic phonon dispersion (blue lines in Fig. 5a) we find a single soft phonon mode exhibiting imaginary phonon frequencies in a BZ region around q 1 = (0, 0, 1/2) = A, clearly visible along the path Γ − A − L, and suggesting an apparent lattice instability. As demonstrated in Ref. [42] for the case of NbS 2 though, imaginary (harmonic) phonon frequencies are not a conclusive indication for lattice instabilities, but rather require a deeper quantum mechanical treatment. For that purpose, we calculated explicitly the corresponding adiabatic potential energy surface (APES) for those q-points that have imaginary phonon modes. This was done via a frozen phonon approach, for which we constructed appropriate supercells in order to fold a specific q-point back to Γ, where the DFPT solution provides a real phonon eigenvector and therefore directly the atomic displacements. In our case, where we find imaginary modes at q 1 = A and q 2 = 2/3 A, this leads to 1×1×2 and 1×1×3 supercells, respectively. We performed total energy calculations as a function of the phonon eigenvector amplitude, allowing us to construct the full anharmonic APES, for which we then solved the one-dimensional Schrödinger equation [92].
To calculate phonon dispersion relations and electronphonon interactions in the presence of anharmonicity, we retain the (harmonic) DFPT eigenmodes and use the anharmonically corrected frequencies, approximated via the difference between the first two eigenenergies of the fully anharmonic APES, to compute the interatomic force constants and dynamical matrices [93]. We also find that taking into account the full anharmonic potential leads to ground state probability densities centered at the highsymmetry structures, i.e. non-displaced atoms, and nonvanishing anharmonic phonon frequencies, and therefore to a stabilization of the Nb 2 S structure. The anharmonically corrected phonon dispersion is presented in Fig. 5a using the solid red line.
In Fig. 5b we report the Eliashberg spectral function α 2 F (ω) and cumulative ep coupling strength λ(ω), including anharmonic corrections. The Eliashberg function essentially follows the phonon DOS for energies above 15 meV, but exhibits a few additional peaks in the low- energy range corresponding to the flat dispersion of the anharmonic mode along the path L − M . In accordance with the sharp increases in the cumulative λ due to the peaks in α 2 F (ω), we find that the low-frequency, anharmonic mode contributes around 50% to the total ep coupling of about λ = 1.5. In order to elucidate the origin of this strong ep interaction, we calculated the mode-and wave vector resolved ep coupling strength λ q,ν . As expected, we find a large portion of the total ep coupling stemming from the anharmonic mode around A and close to L − M (Fig. 5c). Conversely, we observe no strong response at these regions of the BZ for the nesting function ζ q (Fig. 5d), indicating that the observed anharmonicity and softening of phonon modes is purely due to a strong, spatially localized ep coupling.
In Fig. 6a, we trace the anharmonic modes for the two irreducible wave vectors q 1 and q 2 over pressure and find an opposed behaviour, i.e., the frequency of the q 2 mode increases with pressure, while it decreases for q 1 . Extrapolating the energy of the q 1 mode, we find a hypothetical dynamic stability limit around 265 GPa. (We want to note at this point that Nb 3 S already becomes thermodynamically more favourable for pressure above 240 GPa.) An extrapolation of the results for q 2 shows that the energy of this mode will vanish for pressures below 130 GPa.
Decreasing the pressure, we also find that the BZ re-gions of imaginary harmonic modes increase, extending ever closer to and eventually including Γ and additionally appearing on the path Γ−L, as shown in Fig. S6 [87]. At 150 GPa and below, the ground state probability density of the soft mode at q 2 is no longer centered around the high-symmetry structure, but has its maxima in the minima of the double-well APES, thus indicating the transition to a different structure with lower symmetry.
Based on our phonon dispersions for various pressures, we do expect that for pressures between 175 GPa and 200 GPa, Nb 2 S would still be dynamically stable when anharmonicity would be included for all wave vectors. As the corresponding supercells would be very large and the calculations therefore computationally prohibitively expensive, however, we did not correct the additionally appearing imaginary modes below 200 GPa by a full calculation of the APES. Instead, we chose to approximately investigate the sensitivity of the results with respect to anharmonic effects by shifting the soft mode energies of the additional wave vectors "by hand": The choice of our q-grid results in two irreducible q-points in these regions, namely q 3 = 1/3 A and q 4 = (0, 1/3, 1/3) = 2/3 L, for which we set the modes' energies such as to achieve fully positive dispersions in the whole BZ at 175 GPa and 200 GPa, and no soft mode energy below 1 meV in order to avoid artificially high ep couplings.
This allows us to evaluate the superconducting properties over the full stable pressure range of Nb 2 S. In Fig. 6b, we show the total ep coupling strength λ, ω log , and T c as functions of pressure. The circles and solid lines correspond to calculations for which anharmonic effects are only considered for wave vectors q 1 and q 2 (all ep contributions from imaginary phonon modes at q 3 and q 4 are set to zero), while shifting the frequencies of q 3 and q 4 "by hand" as described in the previous paragraph leads to the triangles and dashed curves. The latter yield smoother and dome-shaped curves for T c , supporting our initial hypothesis that the APES associated with the phonon modes at 1/3 A and 2/3 L have indeed a strong anharmonic part.
We find T c values ranging between 12 and 14 K (see Tab. S1 in the SM for a comparison with the other materials in the Nb-S system [87]). λ increases dramatically when decreasing pressure due to increasing anharmonic effects, reaching values comparable to the current T c record holders in the high-pressure hydride class [94][95][96][97][98][99][100]. The characteristic phonon frequency ω log , however, is dominated by low-frequency modes of heavy Nb atoms, thus counteracting the high ep coupling, resulting in only a moderate critical temperature T c .
We further performed calculations in the framework of the fully anisotropic Migdal-Eliashberg theory, in order to obtain the superconducting gap as function of temperature, and find an isotropic and single-gap distribution of the superconducting gap, as detailed in Fig. S5 [87].

V. CONCLUSIONS
In this work we studied the binary phase diagram of the Nb-S system and investigated ep coupling and superconducting properties of its phases up to 250 GPa, using fully ab initio methods. Our results are in excellent agreement with literature for previously reported phases and we find several new materials, revealing an intricate phase diagram at all pressures. Apart from crystal structures consisting mainly of layers with distinct high-and low-pressure building blocks, we also find phases with particularly strong three-dimensional character, and except for the low-pressure phases of pure sulfur, all investigated structures are superconducting metals with T c 's below 15 K.
We focused in particular on the highly threedimensional, high-pressure phase of Nb 2 S that exhibits an apparent lattice instability. By taking into account the full anharmonic potential energy surface, however, we show that this instability is suppressed, demonstrating that imaginary harmonic phonon frequencies are not a conclusive indicator for lattice instabilities, but require a more sophisticated quantum-mechanical treatment. By ruling out Fermi surface nesting, the origin of the strong anharmonicity is found to be in most parts due to an ep interaction particularly strong for a single phonon mode and certain wave vector regions of the Brillouin zone.
Employing the fully anisotropic Migdal-Eliashberg theory, we find that the superconducting gap function is fairly isotropic with only a single gap. We further investigate the influence of the anharmonic corrections on ep coupling and superconductivity, as well as their dependence on pressure, and find a dome-shaped behaviour for T c with a maximum of 14.4 K at 225 GPa.
Our research reveals a much more complex phase diagram for the Nb-S system than indicated by available literature, featuring not only the characteristic, layered vdW structures, but also highly three-dimensional phases at low and high pressures, some of which exhibit interesting anharmonic phonon behaviour and strong ep interaction. Shedding new light on the Nb-S system, our findings encourage to revisit familiar systems with the promise of discovering novel and fascinating phases.

ACKNOWLEDGMENTS
This work was supported by the Austrian Science Fund (FWF) Project No. P 32144-N36, the dCluster of the Graz University of Technology, and the VSC4 of the Vienna University of Technology.
[58] We performed extensive convergence tests on k-point spacing, kinetic energy cutoff and smearing value for the Brillouin zone integration. Converged values were obtained with |∆k| = 0.01 ·2π/Å, a smearing of 0.01 Ry, and an energy cutoff of 125 Ry. We used a convergence threshold of 10 −10 Ry for the total energy in all electronic self-consistence calculations. The thresholds for the relaxation calculations were 10 −7 Ry in total energy and 10 −6 Ry/a0 for all force components. [67] In order to estimate the effect of µ * on Tc, we have calculated Tc for Nb2S at various pressures using the Allen-Dynes McMillan equation for several values µ * up to 0.14, resulting in a maximal total variation of around 2 K (see Fig. 6). In addition, we have also repeated our Migdal-Eliashberg calculations for Nb2S at 225 GPa pressure with µ * = 0.14, leading to a change in Tc also by ∼ 2 K (see Fig. S5  [92] The approach adopted here includes by construction phonon-phonon interactions within the same harmonic eigenmode, but it neglects those between different eigenmodes. This coupling could be included by considering multi-dimensional APES, for which the atomic displacements are expanded in terms of multiple phonon eigenmodes. To do this for all phonon modes is computationally prohibitive, which is why we tested the coupling between the two lowest phonon modes at A and 2/3 A by calculating the two-dimensional APES, as shown in Fig. S27 and Fig. S28

I. THERMODYNAMIC STABILITY
We determined the thermodynamically stable Nb x S y phases (π x,y ) at a fixed external pressure p via their enthalpy of formation ∆H(π x,y ) = H(π x,y ) − xH(π 1,0 ) + yH(π 0,1 ) with respect to the boundary phases of niobium π 1,0 and sulfur π 0,1 , where is the enthalpy per atom, E tot the total energy, N the number of atoms in the unit cell, and V the unit-cell volume.
The sets of phases that are additionally stable against decomposition into neighbouring phases form so-called convex hulls for a fixed pressure, as shown in Fig. 1 of the main text.

II. COMMENTS ON PREVIOUSLY REPORTED PHASES
• Sulfur is an element with many known allotropes. The most stable phases in our calculations are the well known α-S (S 8 rings in space group 70) at ambient pressure, S 6 rings from 5 GPa to 20 GPa (space group 148), spiral tetragonal chains in space group 142 from 20 GPa to 70 GPa, and the β-Po phase (space group 166) above 70 GPa [1][2][3]. We also found and included structures with trigonal spiral chains in one and three orientations [2], as well as other orientations of S 8 , but they did not appear as thermodynamically stable phases in our calculations. The same applies to an incommensurate base-centered orthorhombic (bco) S, whose stability has not been confirmed by theoretical calculations [3][4][5]. We want to note at this point that we found some structures in the transition range from space group 142 to 166 (60 GPa to 75 GPa) formally lower in enthalpy, but with an enthalpy difference smaller than our accuracy of 5 meV. These structures show some similarity to the puckered layers proposed for the bco S modelling [6], yet with lower symmetry.
• We find all known low-pressure phases of NbS 2 (1T-, 1H-, 2H-and 3R) and especially confirm the extensively studied 2H structure as the stable phase at ambient pressure [7,8]. We also confirm the high-pressure phase (space group 139), although we find a slightly higher transition pressure as predicted in Ref. [9] (between 30 GPa and 35 GPa instead of 26 GPa) that can be explained by taking into account vdW as shown in Ref. [10], where a transition pressure of 38 GPa is predicted. We further find imaginary phonon frequencies at q = X below 100 GPa to 150 GPa, which were not predicted in Ref. [9].
• NbS 3 is known in several phases usually (historically) labelled I to V and "HP" (for high pressure, referring to a <10 GPa range). I, IV and V are ZrSe 3 -like structures (space group 11 with 8 atoms in the unit cell) with different distortions, where V is the most symmetric, I is slightly distorted and the distortion in IV reaches over two unit cells. While NbS 3 -I is commonly considered to be the predominant structure formed in experiments, our calculations yield the larger NbSe 3 -like structure "NbS 3 -HP" (24 atoms in the unit cell) to be thermodynamically stable up to 30 GPa (labelled "NbS 3 a" in the phase diagram in Fig. 2 of the main text). From 30 GPa to 55 GPa, NbS 3 -V is more favourable (labelled "NbS 3 b"), for which we confirm imaginary frequencies in different BZ regions and strong ep couplings λ > 2 [11][12][13][14].

III. LOW-PRESSURE GEOMETRIES
The low-pressure structures in the Nb-S system are dominated by trigonal prismatic or octahedral coordination geometries possible for both element types. The coordination number for both geometries is six, stemming from two atomic triangles that are either aligned (prism) or rotated by 60°(octahedron). Figure S1 shows the coordination geometries of Nb and S atoms in the low-pressure NbS phase (space group 194). In many metastable phases, we find the coordination geometries in combination with stacking patterns known from other NbS 2 phases. For example, we report a stable phase for Nb 2 S 3 formed in that fashion from 3R-NbS 2 by doubling the individual vdW layers as shown in Fig. S2. Figure S2. Low pressure Nb2S3 (lower) with trigonal prisms stacked in the pattern from metastable 3R-NbS2 (upper). Both structures belong to space group 160.

IV. HIGH-PRESSURE STACKING
The high-pressure structures in the Nb-S system are dominated by diagonal stackings of simple cubic high-pressure NbS units, forming structures of space group 139 and 123 for even and odd stacking ratios, respectively. Figure S3 illustrates the stacking patterns. Figure S3. Stacking patterns of the simple cubic NbS (center) with even n : n ratios on the left and odd n : (n + 1) ratios on the right. Even structures include NbS2 (n = 1), Nb2S3 (n = 2), and Nb3S4 (n = 3). Odd structures include Nb3S5 (n = 1) and Nb5S7 (n = 2).
Even ratios n : n result in structures with stoichiometric formula Nb n S n+1 , and odd ratios n : (n + 1) in structures Nb 2n+1 S 2n+3 . We found many metastable structures up to n = 5, but only the even stacking for n = 1 and n = 2 resulted in the thermodynamically stable phases of NbS 2 and Nb 2 S 3 . For the sake of curiosity, we constructed structures with different even stackings, i.e. n : (n + 2), resulting in the same n : n formula, but with lower enthalpies of formation.
V. Nb2S Figure S4 shows the electronic DOS in Nb 2 S for pressures from 175 GPa to 250 GPa, emphasizing that the slope around the Fermi energy E F does not change significantly with pressure. Figure S4. Electronic DOS of Nb2S around the Fermi Energy EF for different pressures. Figure S5 shows the results from fully anisotropic Migdal-Eliashberg calculations for Nb 2 S at 225 GPa, revealing a single isotropic superconducting gap closing at ∼25 K. As reported before, the fully anisotropic, mode and wave vector resolved Migdal-Eliashberg calculations lead, in general, to higher T c 's [7,15,16] than the averaged, isotropic Allen-Dynes McMillan formula [17]. Figure S5. Energy distribution of the superconducting gap ∆ k of Nb2S at 225 GPa as a function of temperature, calculated within the anisotropic Migdal-Eliashberg theory as implemented in the EPW code [18], using µ * = 0.1 (black) and µ * = 0.14 (grey). Figure S6 shows the full phonon dispersion curves for Nb 2 S in the studied pressure range, including the anharmonic corrections for q 1,2 , where the opposed behaviour of the modes at A and between Γ − A can be appreciated.

VI. ELECTRON-PHONON PROPERTIES OF STUDIED PHASES
The electron-phonon and superconducting properties of considered, metallic phases are summarized in Tab. S1. Table S1. Electron-phonon and superconducting properties of considered, metallic phases. The asterisk * indicates phases with one or more imaginary frequencies in the harmonic approximation, who's contribution to the ep coupling has been set to zero, the dagger symbol † indicates metastability, and NF is the DOS per spin and per formula unit at the Fermi level. The reported values were obtained with a broadening value of 0.025 Ry.

VII. CRYSTAL STRUCTURE DETAILS
In the following Tables S2-S21 we provide crystal structure details for all considered phases.    In the following figures we report the electronic properties for all studied phases, where the partial DOS for Nb d orbitals is shown in green, Nb s in turquoise, Nb p in blue, and the S p orbitals in red and S s in ochre. We performed phonon calculations for phases for which we could not find corresponding data in literature.         Figure S27. Square moduli of the first six wavefunctions, together with their eigenenergies (with respect to the lowest state) calculated for the 2-dimensional APES of Nb2S at 225 GPa, taking into account the coupling between the first two phonon modes at A. The first four excited states belong to the first phonon mode and the fifth to the second phonon mode. The approximated eigenfrequencies for the two modes (∆ ω1 and ∆ ω5) are in good agreement with the uncoupled, 1-dimensional APES calculations (3.8 and 21.6 meV vs. 3.1 and 20.4 meV, respectively, 1D solutions not shown explicitly). Figure S28. Square moduli of the first six wavefunctions, together with their eigenenergies (with respect to the lowest state) calculated for the 2-dimensional APES of Nb2S at 225 GPa, taking into account the coupling between the first two phonon modes at 2/3 A. The first three and the fifth excited states belong to the first phonon mode and the fourth to the second phonon mode. The approximated eigenfrequencies for the two modes (∆ ω1 and ∆ ω4) are in good agreement with the uncoupled, 1-dimensional APES calculations (5.4 and 22.8 meV vs. 3.5 and 24.9 meV, respectively, 1D solutions not shown explicitly).