Materials synthesis at terapascal static pressures

Theoretical modelling predicts very unusual structures and properties of materials at extreme pressure and temperature conditions1,2. Hitherto, their synthesis and investigation above 200 gigapascals have been hindered both by the technical complexity of ultrahigh-pressure experiments and by the absence of relevant in situ methods of materials analysis. Here we report on a methodology developed to enable experiments at static compression in the terapascal regime with laser heating. We apply this method to realize pressures of about 600 and 900 gigapascals in a laser-heated double-stage diamond anvil cell3, producing a rhenium–nitrogen alloy and achieving the synthesis of rhenium nitride Re7N3—which, as our theoretical analysis shows, is only stable under extreme compression. Full chemical and structural characterization of the materials, realized using synchrotron single-crystal X-ray diffraction on microcrystals in situ, demonstrates the capabilities of the methodology to extend high-pressure crystallography to the terapascal regime.

Due to the difference in electronegativity between Re and N, one can expect that the bonding in the system should be predominantly ionic. On the other hand, the overlap between dorbitals of Re and p-orbitals of N can be considered as an indicator of hybridization and the covalent component of the bonding. Indeed, the presence of unoccupied N p-states above the Fermi energy shows that the interaction cannot be purely ionic. To better understand the nature of atomic bonding in Re7N3, we calculated the electron localization function (ELF) ( Supplementary   Information Fig. S4), spatial distribution of electronic density maps and total charge-density contour (Fig. S5). The low localization of electrons in the region between adjacent N and Re atoms with almost spherically distributed ELF attractors around nitrogen atoms indicates the ionic component of the interatomic bond ( Supplementary Information Fig. S4). The same features can be seen from Supplementary

Re-based solution phase
We use VASP to calculate the mixing enthalpies for 2×1×2 supercells with an underlying fcc crystal structure with 16 Re atoms and various amounts of either N or C (Re16(C,N)x with x=2, 3,4) occupying octapores in the supercells to simulate the Re-N and Re-C cubic phases with NaCl (B1) type structure. To calculate the formation enthalpy ΔH for Re-N solid solutions, we used hcprhenium and nitrogen in the cg-N structure and for Re-C carbon in the diamond structure.
Analysis of the 2×1×2 fcc configurations containing nitrogen shows that it is beneficial for it to be either in neighboring octopores (touching along the edge), or in the octopores that are on the third coordination sphere (that is, those which don't touch at all, but are not far away). Such configurations lead to negative formation enthalpies ( Supplementary Information Fig. S8, Table S5). The situation when the nitrogen octopores touch only in a single vertex (the case, where each unit cell stoichiometry is Re4N) is the least energetically favorable.
We note that for each individual supercell structural relaxation leads to slight orthorhombic distortions (< 4%). However, the presence of several nearly energetically degenerate supercells (Supplementary Information Table S5) suggests that the distortions could be out in a macroscopic sample. Therefore, we perform additional structural relaxations in 2×2×2 supercells using Quantum Espresso 14 . Similar to the VASP calculations , we use the projector augmented wave (PAW) method and consider exchange and correlation effects in terms of the generalized gradient approximation (GGA) with PBE parametrization 4 . We use a kinetic energy cutoff of 80 Ry for wavefunctions and 800 Ry for the charge density and 2×2×2 k-points.
We generate all possible structures based on a 2×2×2 supercell of fcc-Re with stoichiometry Re32N8, which showed the most negative formation enthalpies of the probed 2×1×2 supercells ( Supplementary Information Fig. S8, Supplementary Information Table S5). When occupying 8 out of 32 possible octopores in the cell, it is possible to build 10518300 structures.
We reduce the number by considering only structures in which the center of gravity of the nitrogen coordinates coincides with the center of the supercell, which reduces the data set to ~50000 structures and ensures a relatively uniform distribution of N atoms ( Supplementary Information   Fig. S9). From this set, we randomly choose 100 structures and perform a structural relaxation with no constraints on the lattice or cell parameters until forces are < 0.5 meV/Å 2 per atom. Within the data set, we obtain a small enthalpy variation, depending on the relative position of the octopores occupied by N, in agreement with the 2×1×2unit cell calculations. All calculations maintain the fcc-Re host lattice with varying orthorhombic distortion between 0.4% and 4% and < 0.3% variation in cell volume. Averaging over all coordinates leads to a nearly cubic structure (average orthorhombic distortion is below 0.02%).
Thus, the formation of ReN0.20 solution phase with NaCl (B1) type structure is quite plausible, even though the alloy should be metastable at the synthesis pressure as its mixing enthalpy is above the convex hull. For all ReC solid solutions, the formation enthalpies are positive, making the formation of the Re-C alloys in the experiment highly unlikely ( Supplementary Information Fig. 8, Supplementary Information Table S6).

Thermodynamic stability of Re 7 N 3
Investigation of the influence of pressure on the thermodynamic stability of Re7N3 is a highly non-trivial task. In the maximum pressure range 100 GPa-200 GPa that has been achieved so far, nitrides of rhenium can be found in a wide range of compositions as a consequence of the multiple oxidation states of rhenium, while the phase diagram in the TPa pressure range is unknown and an identification of all the competing phases, as well as a treatment of the offstoichiometric phases is beyond the subject matter of the present study. Therefore, the enthalpies of formation were calculated for Re7N3 and compared to experimentally known and theoretically predicted stoichiometric high-pressure phases in Re-N system from the literature. In more details, we used hcp rhenium in the P63/mmc structure, and nitrogen in the cg-N structure. Zhao et al. 15 previously constructed a convex hull at a pressure of 100 GPa 16 . Using the evolutionary structure search method in a region with a high rhenium content, two thermodynamically stable structures, Re3N ( 6 � 2) and Re2N (P63/mmc), were identified by Zhao et al. theoretically. In fact, these two compounds were synthesized earlier by Friedrich et al. 10 at lower pressures (13-31 GPa). Two nitrides Re3N2 ( 6 � 2) and ReN-NiAs were found in ref. 15 to be quite close to the convex hull, both slightly above it. Thus, these four structures together with Re7N3 were considered in the present study.
In the region enriched with nitrogen we considered the following phases: ReN2 (C2/m), ReN3 (Imm2), ReN4 (Cmmm), ReN2 (P21/c), ReN2 (P4/mbm), ReN10 (Immm). The three former phases, ReN2 (C2/m), ReN3 (Imm2) and ReN4 (Cmmm), were predicted theoretically as thermodynamically stable phases at 100 GPa in ref. 15  Since the enthalpy of formation for Re7N3 is very close to the ground state line, a change in pressure can affect the thermodynamic stability of this nitride. To study this further, we constructed a ground state line at higher pressure of 900 GPa (Fig. 3). In the concentration interval enriched with nitrogen the results changed very little. However, for nitrides with a high rhenium content changes were significant. It can be seen that at this pressure the Re2N and ReN phases are still thermodynamically stable, but Re3N turned out to lie slightly higher than the ground state line.
Most importantly, for Re7N3 the change in pressure placed its enthalpy of formation at the convex hull line. Note that at P=100 GPa Re7N3 has formation enthalpy wich is well above the convex hull line (Fig. 3).

Lattice dynamics of Re 7 N 3
It is necessary that the Re7N3 phase is at least metastable to be synthesized in experiment.
This can be evaluated theoretically in a study of the dynamical stability of a material: in a metastable phase all the frequencies of its lattice vibrations are real. To check the dynamical stability of the Re7N3, we calculated its phonon dispersion relations in the harmonic approximation.
We observed that at 100 GPa the dynamical stability condition is fulfilled (Extended Data Fig. 9), despite the fact that the formation enthalpy of Re7N3 is well above the convex hull (Fig. 3).
Surprisingly, calculations carried out at ~730 GPa, that corresponds to volume per atom at the experimental synthesis pressures, showed the presence of imaginary frequencies along Г-M-K-Г direction of the Brillouin zone (Extended Data Fig. 9). This means that in the harmonic approximation Re7N3 would be classified as dynamically unstable, impossible to synthesize compound. However, calculations of the vibrational spectral function carried out at T=300 K that take into account anharmonic effects of lattice vibrations do not show any sign of the dynamical instability and confirm that Re7N3 is at least metastable at this pressure (Extended Data Fig. 9).
This effect is interesting, giving the fact that the temperature 300 K corresponds to ~25 meV, while the PV term at 1 TPa is ~750 eV. Moreover, Re7N3 is not a strongly anharmonic solid: phonon lifetimes are quite long, as seen from a relatively small broadening of the phonon lines in (Extended Data Fig. 9). The contribution of anharmonic effects of lattice vibrations is expected to be tiny in such materials. To understand this observation, we consider in more details the electronic structure of the Re7N3 and its interplay with lattice dynamics of this system. We recall that Re electron configuration is Xe 4f 14 5d 5 6s 2 . Nine electrons are transferred from seven Re atoms to three N atoms. The remaining sp-and d-electrons are weakly hybridized, leading to a nearly halffilled Re d-band, similar to Re metal. Indeed, comparing the Re d-band of Re7N3 with the DOS of pure hcp-Re calculated at nearly the same pressure ( Supplementary Information Fig. S3). The latter, being half-filled by 5 Re d-electrons that occupy all the bonding states leaving all the antibonding states unoccupied explains very high formation energy of Re-metal in the framework of the Friedel rectangular band model 19 . Even more well-developed separation of the bonding states below the Fermi energy from antibonding states above the Fermi energy by a pseudogap located in a vicinity of EF is seen in the DOS of Re7N3.
In fact, comparing its electronic DOS at the synthesis pressure with that calculated at P~100 GPa one observes ( Supplementary Information Fig. 2), besides a typical pressure induced broadening of the bands, a clear shift of the occupied peaks down in energy (relative to EF, see insets in Supplementary Information Fig. 2). Such a position of the pseudogap, in general, is a characteristic of a stable system, while the pressure-induced shift of the occupied peaks reduces the one-electron contribution to the total energy contributing to the pressure-induced stabilization of Re7N3.
However, as we deal with a chemically complex compound with many atoms per unit cell, the details of the electronic structure are non-trivial. We identify the presence of a Van-Hove singularity seen as a peak of the electronic DOS at the Fermi level (inset in Supplementary   Information Fig. S2), associated with rather flat bands along Г-M line ( Supplementary Information   Fig. S1). Though the singularity is small, its position at the Fermi energy is quite unfavorable from an energetic point of view, contributing to the observed dynamical instability of Re7N3 in the harmonic approximation at P~730 GPa (Extended Data Fig. 9) which effectively uses the static ideal crystal lattice (with small displacements of selected ions) for calculations of the force constants. Note that at P~100 GPa the Van-Hove singularity is above EF, and Re7N3 is predicted to be dynamically stable in the harmonic approximation. On the other hand, it has been established that lattice vibrations at finite temperature smear out the singularities of the electronic structure, leading to a disappearance of anomalies at the phonon dispersion relations 20 . As the Van-Hove singularity at the Fermi level observed at P~730 GPa is quite small in Re7N3, the TDEP calculations at temperature 300 K predict that the materials is dynamically stable at the synthesis pressure (Extended Data Fig. 9).

Brief overview of the double-stage DAC (dsDAC) technique
The ds-DAC original purpose was to generate ultra-high static pressures, those beyond the limit of a conventional DAC of about 300 to maximum 400 GPa. At first, pressures of about 600 GPa were achieved in a ds-DAC 21 , as determined on the lattice parameters of gold using powder X-ray diffraction (XRD). Then static equations of states (EOSes) of several metals (platinum, osmium, tungsten, and tantalum) were cross-calibrated up to 500 GPa 22 based on the EOS of gold by Yokoo et al. 23 and the compressional behavior of Os was studied up to 774 GPa 22 . Hitherto, the highest static pressure achieved using ds-DACs is 1065 GPa 24 , as determined on the gold pressure scale 23 . Modified ds-DAC designs tried by various research groups [25][26][27] have not led to generating pressures comparable to those reported by Bayreuth team 24 . Sakai et al. 28 conducted experiments in a ds-DAC with the secondary anvils manufactured from single crystal diamond using FIB and reported the maxim pressure of 460 GPa. Although these authors achieved the same degree of compression of rhenium as in Dubrovinsky et al. 21 , they reported different pressures, as they used the Re pressure scale of Anzellini et al. 29 . There is a significant inconstancy between the data of Anzellini et al. 29 and Dubrovinsky et al. 21 and this reflects some general difficulties in the pressure characterization at multimegabar pressures. A unique estimate of pressure becomes difficult, if there are discrepancies in the literature. This is the case for Re, as the results by Sakai et al. 28 up to ∼300 GPa pressure are close to those of Anzellini et al. 29 , while Jenei's et al. 30 measurements in the same pressure range agree with the EOS of Dubrovinsky et al. 21 .
Recently diamond anvils of a toroidal shape (t-DACs) were described by Jenei et al. 30 and Dewaele et al. 31 The both groups 30,31 were able to generate static pressures above 600 GPa, thus confirming conclusions of Bayreuth team that pressures above 0.5 TPa may be achieved using the DAC technique. As pointed out by Dewaele et al. 31 , the toroidal shape of the diamond anvil of the t-DAC has similarities with the second-stage diamond anvil of the Bayreuth's ds-DAC 21,22,24 .