First-principles high pressure structure searching, longitudinal-transverse mode coupling and absence of simple cubic phase in sulfur

We use first-principles density functional theory to conduct an extensive structure search using the AIRSS package for elemental sulfur in the range 50–550 GPa. We then obtain the low-temperature phase diagram of sulfur in the same pressure range, including vibrational effects through the harmonic approximation. We do not find any structures lower in energy than those already reported in experiment, although the phase diagram below 100 GPa is found to be crowded with structures separated by only a few meV. We report the transition sequence I 4 1 / acd → P 1 ¯ → ICM → C 2 / m → R 3 ¯ m → I m 3 ¯ m and obtain accurate pressures for each transition, although we find the second-order C 2 / m → R 3 ¯ m transition particularly difficult to define. Contrary to previous first-principles works (Pavel et al; Rudin and Liu 1999 Phys. Rev. Lett. 83 3049–52), we do not reproduce a trigonal → simple cubic transition at either the static lattice or harmonic level. We also undertake a detailed analysis of the incommensurately modulated (ICM) phase of sulfur phase using a commensurate approximant found in the structure search. We find that the modulation amplitude is zero above 96 GPa; some 40 GPa below the experimentally reported transition to the unmodulated phase. We find that the body-centred atoms in the relaxed ICM approximant are, in addition to the dominant transverse modulation (which is a frozen-in optical phonon mode), slightly displaced longitudinally in the b-direction. We subsequently discover that this (small) longitudinal modulation is coupled to the transverse mode, and hence report previously unnoticed weak-mode coupling between transverse and longitudinal optical phonons in the ICM phase.

R m Im m 2 3 3 and obtain accurate pressures for each transition, although we find the second-order C m R m 2 3 transition particularly difficult to define. Contrary to previous first-principles works (Pavel et al; Rudin and Liu 1999 Phys. Rev. Lett. 83 3049-52), we do not reproduce a trigonal  simple cubic transition at either the static lattice or harmonic level. We also undertake a detailed analysis of the incommensurately modulated (ICM) phase of sulfur phase using a commensurate approximant found in the structure search. We find that the modulation amplitude is zero above 96GPa; some 40 GPa below the experimentally reported transition to the unmodulated phase. We find that the body-centred atoms in the relaxed ICM approximant are, in addition to the dominant transverse modulation (which is a frozen-in optical phonon mode), slightly displaced longitudinally in the b-direction. We subsequently discover that this (small) longitudinal modulation is coupled to the transverse mode, and hence report previously unnoticed weak-mode coupling between transverse and longitudinal optical phonons in the ICM phase.

Introduction
The low-temperature phase diagram of elemental sulfur from ambient pressure up to around 200 GPa is now reasonably well understood from experiment. Upon increasing pressure from ambient, sulfur undergoes an orthorhombic (Fddd)  trigonal (P3 21 2 ) transition at around 2 GPa [1,2] (with the Fddd structure being the well-known ambient S-8 ring configuration [3]). At around 36 GPa, the trigonal P3 21 2 structure undergoes a phase transition to a body-centered tetragonal structure with space group I acd 4 1 , consisting of helical square chains that run along the tetragonal c-axis.
Upon further pressure increase, sulfur is experimentally observed to undergo a transition to an incommensurately modulated (ICM) superconducting phase at 83 GPa, which can be viewed as a distorted body-centered monoclinic structure (space group C2/m). This phase is reported to exist, with a modulation amplitude that decreases monotonically with pressure, up to 135 GPa, at which point sulfur adopts the trigonal β-Po structure with space groupR m 3 .
In-between the squared-chain tetragonal and ICM phases, reports of an unknown phase of sulfur have repeatedly surfaced [4,5], with Hejny et al [5] being the first to propose this structure as triclinic.
At higher pressures still, which remain experimentally unexplored, previous first-principles works [6,7] have predicted an¯ R m Pm m 3 3 (simple-cubic) transition in sulfur at around 300GPa. This transition is not, however, observed in the neighbouring Chalcogens Selenium and Tellurium, which are experimentally observed [1,2,8] to instead transform to a body-centered cubic phase.

Structure search
For our structure searching we used the AIRSS package [9] with the density functional theory electronic structure code CASTEP [10]. The PBE functional was used, with a plane-wave cutoff of 500eV and k-point spacing of 0.05 Å −1 . The on-the-fly generated (OTFG) pseudopotentials from Castep were used. The cutoff radius of these pseudopotentials was 0.95 Å, which is much lower than half the interatomic separation at the lower pressures investigated in this paper; at the highest pressures (≈550 GPa), this value approaches half the interatomic separation, and we do not go above this regime.
Searching was carried out at 50 GPa and then 100-500 GPa at 100 GPa intervals. The random input structures built by AIRSS were constrained to randomly contain between 1 and 20 sulfur atoms, and to satisfy from 1 to 48 symmetry operations. The cells were initialised with a volume appropriate to the pressure at which the searching was conducted [11,12]. The simplest AIRSS approach would be to sample structures quasirandomly, but in practice we steer the searches towards more realistic structures which maintain structural diversity. This is achieved by imposing constraints on symmetry, interatomic distances, coordination numbers, stoichiometry, dimensionality, and structural units.
The 50 GPa search yielded ≈13500 structures, and the searches at 100-500 GPa yielded ≈3000 structures per search; a larger number of structures were generated at 50 GPa than at the higher pressures, in anticipation of sulfur's crowded and complicated phase diagram at lower pressures.
The lowest-Enthalpy 10% of structures from each search were then selected as energetically competitive phases, and similar structures were then grouped together using the CRYAN tool bundled with AIRSS, which identifies nearly-identical structures and unites them.
The raw results of our structure search, subjected to the procedure above, are shown in table 1. These raw results were calculated with computationally inexpensive parameters in order to sample the structure space efficiently, and they do not include phonon contributions; therefore, the enthalpy rankings shown in table 1 are different from the final results presented later, which use fully converged parameters (see next section). Nonetheless, these tables provide for a good appreciation of the most energetically competitive structures in configuration space at the pressures concerned and show how often each structure turns up in the search (related to the size of that structure's basin of attraction).

Gibbs free energies
At non-zero temperature and pressure a physical system at equilibrium minimises its Gibbs free energy G(T, p): , , tot elec phonon at various volumes using an applied pressure, and then fit the results to the Vinet equation of state. The Gibbs free energy is then calculated using equation (1). This approach also fully accounts for pressure contributions from (harmonic) phonons. The structure with the lowest Gibbs free energy at a given temperature and pressure will be the most thermodynamically favourable phase (although energetic barriers in configuration space could make these states difficult to access in practice). Using a k-point spacing of 0.010 Å −1 , plane-wave cutoff of 600 eV, a smearing temperature of 300 K and a 4×4×4phonon q-point grid, our Gibbs free energies are shown to be converged to ±0.5 meV below 100 GPa, and around ±1 meV above this (see supplementary material available online at stacks.iop.org/NJP/22/023020/ mmedia).
We first consider the Gibbs free energy landscape in the region 100-550 GPa, and then separately below 100 GPa. We use units of energy per atom throughout. 4. Results and discussion: 100-550 GPa Figure 1 shows Gibbs free energy curves for a number of structures, neglecting zero-point energy contributions, in the pressure range 100-550 GPa, most of which came from the AIRSS searches conducted at 100 GPa and above (a few structures from the 50 GPa search are also included). It can be seen that, at the static lattice level, the trigonalR m 3 (β-Po) structure is the most favourable up until 500±20 GPa, where a first-order transition to theĪm m 3 BCC phase occurs. We do not reproduce a transition to the simple-cubicPm m 3 phase at any pressure, however the relatively close approach (≈10 meV apart) of the two curves at 380 GPa invites us to consider whether zero-point energy contributions could induce a phase transition. Figure 2 shows the Gibbs free energies, both with and without zero-point energies, of the trigonal and simple cubic phases at their point of closest approach. It can be seen that even after the inclusion of zero-point energies, a transition still does not occur, and in fact the effect of this is to increase the stability of theR m 3 phase relative to the SC phase (now by ≈35 meV), owing to the latter's higher ZPE.
The lack of a transition to the simple cubic phase is in disagreement 1 with Gavryushkin et al [6] and Rudin [7], who give pressures of 333 and 280 GPa, respectively, for this supposed transition. However, the absence of this transition would make the phase diagram analogous to those of Selenium and Tellurium, which transform directly from theR m 3 phase to the BCC phase under pressure, without an intervening transition to the SC phase. Unlike sulfur, which has not yet been observed beyond the trigonal phase, this direct¯ R m 3 BCC transition has actually been observed for Selenium and Tellurium in experiments [1,2], as a much lower pressure is required to induce the transition in these other Chalcogens. Figure 3 shows that theR m 3 phase remains the most thermomechanically favourable until ≈500 GPa, where sulfur undergoes a first order phase transition to the BCCĪm m 3 phase. Inclusion of zero-point energies reduces this pressure to 476±3 GPa, due to the greater ZPE of theR m 3 phase relative to theĪm m 3 phase.

5.
Results and discussion: < 100 GPa Figure 4 shows the Gibbs free energy of various sulfur phases below 100 GPa, neglecting and including zeropoint energies, respectively. The structure labelled ICM is an approximant to the experimentally observed incommensurately modulated phase of sulfur (see figure 5), the unmodulated form of which has C2/m symmetry. This approximant was actually found in the AIRSS search-only once-and is the P2 1 m structure ranked seventh in the 50 GPa search in table 1. Figures 5-9 illustrate some of the sulfur phases below 100 GPa. The structure designated C2/m-A is the (unmodulated) monoclinic structure given by Degtyareva et al, geometry optimised for each pressure. The structure labelled C2/m-B is a different structure, also with C2/m symmetry, consisting of zig-zag sulfur chains. The Cm andP 1 structures are very similar in appearance to the C2/m-B structure, also consisting of zig-zag chains, but are slightly distorted and have lower symmetry. TheP 1 variant of this structure has particularly low energy, and is in fact, in a relatively narrow pressure range (see below), the ground state of sulfur.
It can be seen that inclusion of zero-point energies serves only to change the pressure boundaries of each transition by a few GPa. Having accounted for (harmonic) vibrations, we obtain the transition sequencē   I acd P 4 1 1 ICM with increasing pressure, at 56.1±0.1 GPa and 69.2±0.5 GPa respectively. This transition sequence is in good agreement with experiment 2 [5]. The pressure boundaries for each transition depend very weakly on temperature, with each transition pressure being lowered by around 1 GPa when increasing temperature from 0 K to 300 K.
The static-lattice Gibbs free energy curves in figure 4 seem to indicate, contrary to experiment [1,2,5], that the incommensurately modulated phase survives only to ≈85 GPa, at which it undergoes a transition to either its  Rudin [7] uses LDA, so it may be tempting to conclude that this disagreement arises from the particular functional used. Gavryushkin et al [6], however, also use PBE like ourselves, and still report an¯ R m 3 SC transition. 2  unmodulated form (monoclinic C2/m) or a trigonalR m 3 phase (these two structures are so similar that we were unable to energetically separate the structures anywhere in the Gibbs free energy landscape; see section 5.1).
Inclusion of zero-point energies increases the stability range of the ICM phase, as it is found that the C2/m and R m 3 structures exhibit imaginary phonon modes below ≈96 GPa, which distort them both into the ICM phase. However, these modes become real above 96 GPa, indicating that the modulation now raises the energy, and thus implying that the stability range for the ICM phase within the harmonic approximation is 69.2-96 GPa. It should be noted that since the energetic differences between the ICM, C2/m andR m 3 phases in this region are beyond our convergence limits, the imaginary  real transition of the phonon modes, described in the next section, provides a much clearer indication of phase stability. Experimentally, the ICM phase is observed to survive up to 135 GPa .
The non-smooth joining of the Gibbs free energy curves on the right of figure 4 is an artefact of using a commensurate approximant; the truly incommensurate structure would have a slightly lower energy, and the curve would join the unmodulated (C2/m-A) curve smoothly at ≈96 GPa. This is also confirmed by the width of the imaginary phonon mode; between 93 and 95 GPa, the mode is still imaginary, however its width is such that the q-vector corresponding to a four-cell long modulation now has a positive, not negative, energy. The energies are normalised to the trigonalR m 3 structure, and the fitting error is shown as a shaded region for each curve. A square indicates that beyond this pressure, the structure spontaneously collapsed into another symmetry upon geometry optimisation, and a circle indicates that no data was taken beyond this point. The absolute energies were converged to ±1 meV, although the fitting error in a few structures, such as the / I mmm 4 case, is larger than this by up to a factor of five. Nonetheless, this is still more than sufficient at this scale to determine the energy separations.      . TriclinicP 1 structure, viewed along the a axis at 60 GPa. This structure is the ground state of sulfur between 56 and 69 GPa, consisting of zig-zag chains stacked on top of one another in the a-axis direction, and is visually identical at this scale to the C2/m-B structure. It can be seen that the bonds running diagonally bottom-left to top-right are slightly longer than those running top-left to bottom right. Each sulfur atom is fourfold coordinated; in addition to the bonds shown here, each S atom is also bonded to an S atom into the page, and one out-of the page.

Similarity of unmodulated C2/m andR 3m structures
The unmodulated C2/m structure (C2/m-A), which we find the ICM phase changes into above 96 GPa, is very nearly equivalent to the trigonalR m 3 structure. This can most easily be seen by noting that a one-atom rhombohedral primitive cell can be chosen for the C2/m structure, in which two of the lattice lengths are equivalent and the other is not; this is a distortion of a perfect rhombohedron, in which all three lengths are equivalent.
In their ab-initio paper, Gavryushkin et al [6] present that, up to a pressure of ≈250GPa, the inequivalent rhombohedral length differs from the other two lattice parameters by ≈0.06 Å, which would certainly make the C2/m structure noticeably distinct from theR m 3 structure. Above 250 GPa, Gavryushkin et al report that the C2/ m structure suddenly adopts theR m 3 symmetry, thus giving a clear C m R m 2 3 transition. Gavryushkin et al obtain these structures by geometry optimisation of the experimental structure given by Degtyareva et al in [1,2].
We have found however, that if the unmodulated C2/m structure is geometry optimised to a high tolerance, it always transforms to a structure with a negligably small deviation from perfect trigonal symmetry, even when no assumptions are made in the geometry optimisation about the symmetry of the structure (which can bias the optimisation). For example, at 150 GPa, where Gavryushkin et al present the rhombohedral lattice parameters for the C2/m structure as a=b=2.27 Å, c=2.21 Å (a difference of 0.06 Å), we obtain a=b=2.156, c=2.154 (a difference of 0.002 Å); at higher pressures, the deviation is smaller still; reducing to 0.0001 Å at 180 GPa. Furthermore, the abrupt change in lattice parameters reported by Gavryushkin et al is not observed in experiment [1,2].
We are able to partially recover 3 some of the results of Gavryushkin et al if we use their k-point spacing of 0.035 Å −1 , which we have found to be too coarse to converge to anything better than 20 meV at these pressures Figure 8. The monoclinic Cm structure at 60 GPa. This structure is closely related to theP 1 and C2/m-B structures via a distortion. Unlike in theP 1 and C2/m-B cases, the distortion here is slightly larger and is enough to create new (inter-chain) bonds; some of the sulfur atoms are now fivefold coordinated, as opposed to being all fourfold coordinated. Figure 9. The V-shaped corner units that connect the zig-zag chains, shown here for the C2/m-B structure at 60 GPa, and tilted slightly into the page. The C2/m-B, Cm andP 1 structures all feature this corner unit, with the latter two structures being symmetrylowering distortions of the first. TheP 1 structure, for example, is obtained by lengthening the shortest bond (0.205 nm bond above) by ≈1.5%, whilst keeping the other bonds fixed. These distortions are small enough that the C2/m-B andP 1 structures are visually identical, although the Cm structure is distorted just enough for new bonds to be created (see figures 7 and 8). 3 We can reproduce the non-negligable deviation from trigonal symmetry of the C2/m structure, but still cannot reproduce the 'abrupt change in lattice parameters'.  (we used 0.010 Å −1 ). We also note that Gavryushkin et al used a geometry force tolerance of 1×10 −3 eVÅ −1 , which we find also allows the geometry optimiser to relax to a structure much further from perfect trigonal symmetry, giving rise to the apparent large distinction between the C2/m andR m 3 structures (we used a tolerance of 1×10 −5 eVÅ −1 ).
We are therefore, through comparing lattice parameters or differences in Gibbs free energies, unable to provide a transition pressure for the unmodulated C m R m 2 3 transition due to the near complete equivalence of the energetic and structural properties of these phases. However, a similar problem arises experimentally in detecting this transition; Luo et al [4] were able to determine the transition by observing a discontinuity in the first derivative of the volume as a function of pressure, and fitting to a different equation of state either side of this discontinuity.
We are able to replicate this method with success; the fitting error in V(p) using two Birch-Murnaghan fits is reduced than when using just one, indicating that there are two different phases either side of the discontinuity in the derivative of V(p) (to wit: a two-phase fit gives a smaller error than a one-phase fit). For this second order transition, we obtain a pressure of 194.5±12 GPa using this approach, as shown in figure 10. Experimentally, Luo et al and Degtyareva et al obtain 162 and 153 GPa respectively. For comparison, in their ab-initio paper, Gavryushkin et al obtain a pressure of 250-270 GPa for this transition, although they determine this using the aforementioned abrupt change in lattice parameters that we do not observe.
It should be remembered that, even at pressures well below this 'transition' within the stability range of the unmodulated C2/m structure, we find the deviation from trigonal symmetry to be extremely small.
The C m R m 2 3 transition is also observed in Selenium and Tellurium [13], where the strongly second order nature of the phase transition again makes the small discontinuity in the gradient of V(p) the clearest way in which to detect the transition.  6. Phonon dispersion curves for C2/m phase Figure 11 shows harmonic phonon dispersion curves in the [010] direction of the two-atom monoclinic cell (space group C2/m), and figure 12 shows the amplitude and q-vector dependence as a function of pressure. The phonon dispersion curves were obtained by constructing a primitive cell of the C2/m structure that had a reciprocal lattice vector in the [010] (with respect to the monoclinic cell) direction, and then running a 1x25x1 phonon q-point calculation, giving a very finely resolved phonon dispersion. In the top image, the lowest-lying optical branch can be seen to exhibit a negative phonon mode; this distortion pushes the unmodulated structure into the ICM phase. Above 96 GPa, we find that the mode becomes real, implying that the unmodulated structure is now stable. This coincides with the vanishing of the modulation amplitude in the geometryoptimised ICM approximant, as shown in figure 13. We note that Degtyareva et al, in their first-principles paper [14] (also using PBE), report the mode to be imaginary until 135 GPa, which is in better agreement with experiment. However, we have shown here that the modulation vanishes at 96 GPa using three separate approaches: Gibbs free energy considerations, observation of the geometry-optimised modulation amplitude, and the imaginary  real transition of the phonon mode-All of which give the same result.

Longitudinal-transverse phonon mode coupling
It was remarked earlier that, upon optimisation, the body-centred atoms in the ICM approximant were displaced slightly parallel to the b-direction. Figure 11 shows high-precision phonon dispersion curves in the direction of interest ([010]) at two pressures. These were obtained by constructing a primitive cell of the C2/m structure that had a reciprocal lattice vector in the [010] (with respect to the monoclinic cell) direction, and then running a 1x25x1 phonon q-point calculation, giving a very finely resolved phonon dispersion.
The phonon dispersion curves in figure 11 however, show only a softening of a single transverse optical branch, apparently suggesting that any longitudinal modes serve only to raise the energy.
In figure 14 below, we show that this longitudinal displacement of the body-centre atoms is coupled to the transverse optical mode; that is, the longitudinal displacement only lowers the energy of the structure when it is simultaneously displaced along the transverse optical mode-by itself, the longitudinal mode has a positive energy. Figure 14 also shows that at 70 GPa, where the ICM phase has just become the lowest energy structure, the longitudinal modulation (with amplitude ≈0.010 Å) lowers the energy by only an extra 0.35 meV compared to when no longitudinal modulation is applied, whereas the transverse modulation lowers the energy by ≈3.7 meV relative to the completely unmodulated structure.
The longitudinal displacement and corresponding energy reduction is larger at lower pressures (for example, 0.016 Åand 0.6 meV respectively at 50 GPa), however the ICM phase is not the lowest energy structure at these pressures. Therefore, if experimentally detectable, the longitudinal displacement is most likely to be observed at 70 GPa; the lowest pressure at which the ICM phase is the most favourable structure.

Electronic bandstructures
The transverse modulation of the C2/m structure opens up a gap [14,15] at the Fermi level along the [010] direction, shown below in figure 15. The size of this band gap decreases with increasing pressure, and becomes zero when the modulation vanishes at 96GPa. The Fermi energy is also lowered relative to the unmodulated case by ≈0.151 eV. direction. These curves are presented at 50 GPa, where the ICM phase is not the lowest energy structure (but is still a stable structure), since the band gap is much larger at this pressure and thus easier to see. The Fermi energy can also clearly be seen to reduce after application of the modulation.  The much smaller longitudinal modulation does not change the electronic bandstructure significantly and does not reorder any bands or open any gaps. Rather, its effect is to curve the conduction band minimum upward so that it intersects the Γ point at a higher energy (see figure 16). The longitudinal modulation further lowers the Fermi energy by ≈1.6 meV relative to the solely transverse modulated structure. Sulfur remains a metal in the ICM phase; although there is a gap in the [010] direction, the electronic density of states, shown in figure 17 is non-zero at the Fermi level due to contributions from states in other directions.

Simulated diffraction patterns
We present simulated diffraction patterns using the DEBYER [16] x-ray diffraction code for the ICM approximant at 75, 85 and 95GPa in figure 18, using the same wavelength as Degtyareva et al.
The most significant features of the ICM diffraction pattern, which vanish in the unmodulated pattern at 95 GPa, are (i) the presence of a small peak at 2θ≈8 (marked with a red ellipse in figure 18), and a small peak at 2θ≈10.9 (in-between the two large peaks) and Since the small peaks at 2θ≈8 and 2θ≈10.9, which in experiment (with slightly shifted 2θ values) serve as very clear identifiers of the ICM phase (as they are completely absent in the unmodulated pattern), are experimentally present up to ≈135 GPa, this would indeed suggest that the ICM phase survives up to a greater pressure than the 96 GPa we present here. Table 2 gives a summary of the transition pressures that we have derived.

Conclusions
We do not find a¯ R m 3 SC transition anywhere, instead finding that, after inclusion of ZPEs, theR m 3 phase is stable by ≈35 meV (although there is still no transition at the static lattice level). Given that this transition is reported by at least two different authors [6,7] with a pressure at the upper end of what is accessible with a diamond anvil cell setup, experimental work on sulfur in this pressure regime is required to resolve this discrepancy. The lack of such a transition would make sulfur analogous to the heavier Chalcogens, Selenium and Tellurium.
Our transition pressures for the I acd P 4 1 1 and¯ P1 ICM transitions are in good agreement with experiment, and our pressure for which theĪm m 3 phase becomes most favourable is comparable to previous works (although we haveR m 3 , notPm m 3 , as the preceeding phase). In disagreement with experiment, we find, through consideration of both Gibbs free energies and phonon dispersion curves of the unmodulated structure, that the ICM phase reverts to its unmodulated C2/m form at 96 GPa. The magnitude of the modulation can also be seen to reduce to zero very close to this pressure by analysing the geometry optimised approximant.
We are able to identify the strongly second order C m R m 2 3 transition through a small discontinuity in V(p), although the uncertainty in this value is relatively large. In disagreement with [6], we do not find the C2/m andR m 3 phases to be appreciably different. We have discovered previously unreported weak phonon coupling between transverse and longitudinal optical modes in the C2/m structure; displacement along the transverse mode gives by far the greatest contribution to the energy reduction relative to the unmodulated phase, but displacing along the longitudinal optical mode simultaneously provides a greater reduction still. This extra energy reduction (and corresponding longitudinal displacement) is greatest at lower pressures, and would be most detectable at 70 GPa (the lowest pressure at which the ICM phase is most favourable).