Novel Stable Compounds in the C-H-O Ternary System at High Pressure

The chemistry of the elements is heavily altered by high pressure, with stabilization of many new and often unexpected compounds, the emergence of which can profoundly change models of planetary interiors, where high pressure reigns. The C-H-O system is one of the most important planet-forming systems, but its high-pressure chemistry is not well known. Here, using state-of-the-art variable-composition evolutionary searches combined with quantum-mechanical calculations, we explore the C-H-O system at pressures up to 400 GPa. Besides uncovering new stable polymorphs of high-pressure elements and known molecules, we predicted the formation of new compounds. A 2CH4:3H2 inclusion compound forms at low pressure and remains stable up to 215 GPa. Carbonic acid (H2CO3), highly unstable at ambient conditions, was predicted to form exothermically at mild pressure (about 1 GPa). As pressure rises, it polymerizes and, above 314 GPa, reacts with water to form orthocarbonic acid (H4CO4). This unexpected high-pressure chemistry is rationalized by analyzing charge density and electron localization function distributions, and implications for general chemistry and planetary science are also discussed.

known compounds made by a framework of hydrogen-bonded water molecules encapsulating guest molecules. Common guests within the C-H-O composition are oxygen, hydrogen, carbon dioxide and small hydrocarbons. All these inclusion compounds are known to decompose when pressure rises beyond certain limits, typically a few GPa. However, some exceptions were reported 16 , e.g. the H 2 O:2H 2 complex was predicted to be stable "at least up to 120 GPa 17 ". Within inclusion compounds, host and guest retain their molecular identities and interact through weak van der Waals forces. To the best of our knowledge, no thermodynamically stable molecules formed solely by C, H and O altogether are known to date.
In this work, the C-H-O phase diagram is explored in the pressure range 10-400 GPa by means of the powerful variable-composition evolutionary algorithm USPEX 18 coupled with periodic density-functional calculations. Thorough, unbiased searches were performed sampling all possible C-H-O compositions, and a total of more than 125000 structures, generated by the evolutionary method, were relaxed to the closest minimum-enthalpy configuration. Chemical bonding in the resulting compounds was investigated by analyzing the real-space distribution of their charge density, in the framework of the Quantum Theory of Atoms in Molecules (QTAIM 19 ), and of their electron localization function (ELF 20 ). According to QTAIM, chemical bonds are mirrored in the charge density distribution at special saddle points, called 'bond critical points' (bcps). The evaluation of certain scalar properties at bcps provides information about the bond type. ELF is a simple measure of electron localization. Maxima appear in regions associated with chemical bonds, lone pairs and atomic cores, and electronic populations can be assigned to these chemical entities by integrating charge density within the basins corresponding to those maxima.
A schematic representation of the C-H-O pressure-composition phase diagram is reported in Fig. 1a,b. Our calculations recovered all the stable compounds known from previous crystal structure prediction works [10][11][12]17,[21][22][23] . The crystal structures we obtained were either the same as those previously reported, or structurally and enthalpically (typically within 1 meV/atom) similar. We briefly summarize these findings ( (Fig. 2a) and the reported one have different H 2 orientation but their enthalpy is nearly indistinguishable (Δ E ≤ 0.5 meV/atom), suggesting that the guest molecules are rotationally disordered even at high pressure. For this compound, we determined for the first time its decomposition pressure: 153 GPa. We discovered a structure of ethane (Fig. 2e) that is more enthalpically favorable than the previously reported one. This reduces the previous estimate 12 of the pressure of the decomposition reaction 2CH 4 → C 2 H 6 + H 2 from 200 GPa down to 121 GPa. For water ice above 200 GPa, we obtained the Pbcm structure already For each compound, we report the space group(s) and the predicted phase transition pressures (in GPa), the latter in round brackets. Numbers in square (curly) brackets indicate the formation (decomposition) pressures. In (a), the formation pressure is not reported for those compounds already stable at 10 GPa. For elements, a detailed description of the known phases and the comparison with the ones obtained in our calculations is reported in Supplementary Section S2. A complete representation of the C-H-O phase diagram at various pressures, including Gibbs triangles, is reported in Supplementary Fig. S6. We also report the convex hulls at various temperatures for the H 2 O-CO 2 system at 1 GPa (c) and 2 GPa (d). 'no ZPE' represents the value obtained neglecting the zero-point (vibrational) energy correction. proposed in previous ab initio studies 24,25 (Fig. 2d). Concerning oxygen, our high-pressure calculations uncovered a previously unknown hexagonal structure reported in Fig. 2b (P6 3 /mmc space group), stable above 375 GPa. Similarly to its lower-pressure ξ -phase (space group C2/m) 21 , this P6 3 /mmc oxygen is a metallic molecular crystal ( Supplementary Fig. S2). In the lower-pressure ε -and ξ -phases, the intermolecular distances within the ab crystallographic plane have different values, ε -phase even having exotic (O 2 ) 4 clusters 21 . In the P6 3 /mmc structure, instead, they become equal thereby conferring to the crystal its hexagonal symmetry.
New compounds were predicted to form. A 2CH 4 :3H 2 (Fig. 2c) clathrate was found to be stable from < 10 GPa up to 215 GPa. It has the same crystal structure throughout its pressure range of stability. Curiously, the topology of the host framework ( Supplementary Fig. S13) is reminiscent of that of gas hydrates. Concerning the possible formation of molecules containing C, H and O, our results indicate that carbonic acid becomes thermodynamically stable above 0.95 GPa (effects of zero-point energy correction and temperature are discussed below). This was quite unexpected as this molecule is highly unstable at ambient conditions 26 . Its synthesis requires the use of high-energy radiation or strong acids and the resulting compound can only be isolated under high vacuum or in argon matrix and at very low temperatures 27,28 . Indeed, the decomposition of carbonic acid into water and carbon dioxide is highly exothermic, although in the absence of water it is hampered by a high kinetic barrier 29 . The low-pressure structure (space group Pnma, Fig. 3a) is composed of chains made by hydrogen-bonded molecules. The latter adopt an almost flat conformation in which hydrogens are in cis position with respect to the C=O bond. This type of HBs arrangement was shown to be the most energetically favorable at ambient conditions, 'cis-cis' being the most stable conformation of the isolated molecule 28 . The crystal structure of carbonic acid is significantly denser than water ice (2.147 vs 1.560 g/cm 3 at 1 GPa). This fact has important implications for planetary science (see below). As pressure rises, carbonic acid polymerizes (p > 44 GPa), forming the structure shown in Fig. 3b. The -CO-backbone of each polymer is parallel to the c crystallographic axis, while hydroxyl groups form HBs joining the adjacent polymeric chains along the b direction. This type of structure is stable up to the highest investigated pressure (400 GPa), although the HB network rearranges above 240 GPa, leaving the CO backbone and the crystallographic space group unchanged. In the higher-pressure conformation (Fig. 3d), each polymer is hydrogen-bonded to the nearest neighbors along the two ab diagonals. This implies that, differently from the lower-pressure Cmc2 1 structure, at high pressure the HB network joins all the polymers together. Finally, at 314 GPa, we detected an exothermic reaction between carbonic acid and water, leading to the formation of orthocarbonic acid (H 4 CO 4 , Fig. 3e). Besides being thermodynamically stable, the newly discovered compounds were ascertained to be dynamically stable by the absence of imaginary frequencies in their phonon dispersion curves (Supplementary Section S6).
The intricate chemistry described above can be rationalized by taking a closer look at the crystal structures and by analyzing the QTAIM properties and ELF distribution. First we note that, as molecular crystal Pnma-H 2 CO 3 is compressed, the C-O hydroxyl bonds (C=O carbonyl bonds) shorten (elongate) and their bcps ellipticity increases (decreases), as shown in Supplementary Table S2. This indicates a progressive delocalization of the π orbital over the two formally single C-O1 bonds, eventually making them shorter than the formally double C=O2 bond. This behavior can be explained with the chemical scheme of Fig. 3c. As pressure rises, the weight of the resonance forms II and III increases. ELF distribution supports this hypothesis: a fourth maximum appears around O2 above 40 GPa (Fig. 4c-f), indicating a progressive shift from sp 2 to sp 3 hybridization, and at  Table S3). This trend can be seen as a pressure-induced destabilization of the double bond, which is likely to be an important factor contributing to the polymerization of carbonic acid. The latter is indeed associated with the breaking of the π bond and the concurrent formation of a new C-O σ bond, in a way which is reminiscent of the high-pressure behavior of CO 2 10 . Accordingly, upon the Pnma→ Cmc2 1 phase transition, C-O bonds elongate and their ellipticity at bcp decreases (Table 1). HBs undergo important changes, too. Differently from the Pnma phase, in the Cmc2 1 structure O-H and H···O interactions display similar bond lengths, degree of covalency (as measured by the kinetic energy density per electron 30 ) and charge density at bcp (Table 1). Moreover, in passing from Pnma to Cmc2 1 phase, a significant decrease (increase) in the population of the ELF basins corresponding to acceptor lone pairs (O-H bonds) takes place. Such a charge transfer from the acceptor lone pair to the donor-H fragment generally occurs in passing from weak to strong HBs 31 . Thus, along the Pnma→ Cmc2 1 phase transition, HBs strengthen and become noticeably more symmetric. The symmetrization of HBs increases as pressure rises, and in the Cmc2 1 phase at 400 GPa O-H and H···O interactions are nearly indistinguishable (Table 1; in the following, these types of HBs will be referred to as 'O-H-O bonds' , to distinguish them from weaker HBs). The different nature of HBs in the two phases of carbonic acid is even more evident in the ELF distribution (Fig. 4a,b), which for Cmc2 1 is more akin to that of symmetric O-H-O bonds observed in the high-pressure phases of ice (ice X and Pbcm-ice, Supplementary Fig. S5). Noteworthy is the appearance of an ELF maximum on H atoms for all the investigated high-pressure HBs, a feature observed at ambient pressure only in 3c-4e bonds such as those of  (Fig. 5). Such HBs are not present in water, hence their number increases in passing from the reactant mixture to H 4 CO 4 . The formation of additional interactions is expected to be associated with a more efficient packing; in fact, by comparing the ELF basins of orthocarbonic acid with those of the reactants (Fig. 6), it emerges that the O-H-O bonds of water significantly shrink along the formation reaction, thereby overcompensating the volume increase observed for other basins. These results suggest that the occurrence of additional HBs upon the formation of orthocarbonic . We represent as dotted, blue lines those hydrogen bonds having OHO angle and H···O distance greater than 150° and shorter than 1.3 Å, respectively. For all the structures having more than one symmetryindependent atom of each type, the labels adopted in Figs 4,5 and 6 and Table 1 are shown. acid plays a fundamental role in the achievement of the volume reduction responsible for its high-pressure stability. We observed similar results (formation of O···H contacts and ELF trends) in a metastable polymorph of H 4 CO 4 (Supplementary Section S8), which is structurally different from the structure discussed above, but is also enthalpically favorable over the H 2 O + H 2 CO 3 mixture above 395 GPa.
We now move to discuss the implications of our results for general chemistry and planetary science. The long-standing view that inclusion compounds systematically decompose at low pressures of a few GPa has been     refuted in light of a number of counter-examples discovered during the last decades. However, only hydrates were known to persist above 50 GPa 16 . Our 2CH 4 :3H 2 compound not only sets a new upper limit for the stability of inclusion compounds, but also introduces a qualitative shift of views, for it broadens the classes of inclusion compounds stable at very high pressures. In fact, several nCH 4 :mH 2 crystals were experimentally found to be stable up to 8 GPa, and CH 4 :2H 2 was compressed up to 30 GPa without showing any sign of decomposition 14 . The difference in the stoichiometry between such compound and the one presented here might be traced back to a number of causes, for example the possible metastability of the experimentally detected phase above 10 GPa. Remarkably, the existence of 2CH 4 :3H 2 affects the high-pressure chemistry of methane: i) it lowers the decomposition pressure of pure methane crystals down to 93 GPa, much lower than previous estimates 12 (Fig. 1a), and ii) allows methane molecules to survive up to 215 GPa (Fig. 1b), i.e. at a higher pressure higher than previously reported. Concerning carbonic acid, its discovered stability already at moderate pressure opens new important possibilities for synthesizing and stabilizing this elusive molecule. However, it must be pointed out that direct comparison between the calculated and experimental formation pressures is complicated by two problems: the influence of lattice vibrations (temperature effects and zero-point energy correction) and the approximations made within the DFT approach. We have tackled these two issues by evaluating how the formation pressure of H 2 CO 3 is affected by phonons (within the harmonic approximation) and by the use of different computational settings (including changes in the basis set and exchange-correlation functional). The pressure required for stabilizing carbonic acid shifts to 1.45 GPa when zero-point vibrational energy is included, and then weakly increases with temperature (Fig. 1c,d). All the tested computational approaches reproduced the pressure-induced stabilization of carbonic acid, although the formation pressure showed some variations (Supplementary Table S4). Overall, a realistic estimate of the formation pressure of carbonic acid would be the range 0.6-1.6 (0.75-1.75) GPa at 100 (300) K. Similar pressures occur on the bed of water oceans of icy satellites 7 . There, both water ice and carbon dioxide are present, hence carbonic acid is likely to form. Moreover, its high density implies that, once formed, carbonic acid will sink to the bottom of ice layers, just above the rocky cores, thereby experiencing an even greater pressure. In such a scenario, carbonic acid insulates water ice from the core. This would modify the chemical composition models for these celestial bodies, which now include possible reactions between water and rocky compounds such as ferromagnesiansilicates 7 .
Very recently (in fact, while this paper was being finalized) Wang et al. 33 detected the formation of carbonic acid in compressed H 2 O:CO 2 mixtures, which confirms our prediction. The authors hypothesized carbonic acid to be a metastable product of a complex, non-equilibrium process, while our findings lay down a theoretical ground for this discovery and assert carbonic acid as a thermodynamically stable compound under pressure. Our calculated IR and Raman spectra agree well with the measured ones (Supplementary Section S9, which contains a detailed comparison between experimental findings from literature and our theoretical results). This indicates that the H 2 CO 3 crystal structure presented here is either the observed one, which is still not determined experimentally, or one with a very similar crystal packing. In fact, the energy landscape of H 2 CO 3 is very complex and features many low-lying minima with comparable energies 34 , hence a systematic exploration of all of them is needed in order to reliably determine the correct structure. Such extensive analysis will be the subject of our future research.
In conclusion, we have carried out a thorough DFT investigation of the C-H-O phase diagram up to 400 GPa. For each composition, the most stable structures at various pressures were obtained by the powerful evolutionary algorithm USPEX. Our calculations predicted the formation of several new compounds. An inclusion compound, 2CH 4 :3H 2 , was found to be stable up to the unprecedented pressure of 215 GPa. Carbonic acid was predicted to be stable above 1 GPa, remaining stable throughout the investigated pressure range and polymerizing above 44 GPa. At 314 GPa it reacts with water to form orthocarbonic acid, H 4 CO 4 . An extensive chemical bonding analysis was performed, which provided important insights such as the pressure-induced stabilization of resonance forms of molecular carbonic acid and the role played by HBs in the stabilization of orthocarbonic acid. This novel chemistry should have major implications for planetary science.

Methods
The variable-composition evolutionary approach implemented in USPEX, scans the structural and chemical spaces and seeks the thermodynamically stable compounds. The first generation of structures was mostly produced randomly (stable structures obtained from previous runs were also added). The subsequent generations were produced by both symmetric random generator and by applying variation operators (described in ref. 35) to the most promising structures (65% of the current generation are allowed to produce the next generation). Details of all USPEX runs can be found in Supplementary Table S1. Exhaustive description of the method 35,36 and studies where the approach reliability was confirmed by comparison with experiments 3,12,21,37 are reported elsewhere. The VASP code 38 was used for structure relaxations and enthalpy calculations. We adopted the Perdew-Burke-Ernzerhof (PBE) functional 39 in the framework of the all-electron projector augmented wave (PAW) method 40 , with 'hard' PAW potentials, plane wave kinetic energy cutoff of 850 eV and a uniform, Γ -centered grid with 2π *0.056 Å −1 spacing for reciprocal space sampling. In the 0-10 GPa range, where dispersion interactions play an important role, we employed for VASP calculations a van der Waals functional (optB88-vdw 41 ). At these pressures, a slightly looser reciprocal space sampling of 2π *0.064 Å −1 was adopted. Note that at each pressure the stability of each compound was ascertained by considering the most stable form of the reactants, including structures from literature. To perform the chemical bonding analysis, we carried out single-point calculations (at the geometry obtained from VASP) with the CRYSTAL14 code 42 . Within the latter, crystal orbitals are described in terms of atom-centered functions. The employed basis set was of 'triple-ξ plus polarization' quality, whose functions were optimized for solid-state calculations 43 . In order to better describe intermolecular interactions, we augmented the basis set with d orbitals on H atoms taken from ref. 44. To make the basis set apt for high-pressure calculations, the exponents of the radial part of the outermost functions of each shell were contracted by a factor 1.44 and 1.69 for calculations below and above 200 GPa, respectively. The bond critical point analysis and the determination of atomic basins within the QTAIM framework was done using the TOPOND code, now implemented within CRYSTAL14. The latter was also exploited to evaluate grid files of scalar properties (such as charge density and ELF), subsequently converted in the standard 'cube' format using the NCImilano code 45 . The integration of quantities within ELF basins, not implemented in TOPOND, was performed with the critic2 code 46 , and in particular the grid-based Yu-Trinkle algorithm 47 was exploited. The input grids were 400 × 400 × 400, spanning the whole unit cell. Phonon dispersion curves, and phonon contribution to the free energy of formation were calculated by means of the finite displacement method implemented in the PHONOPY code 48 . Further details concerning phonon calculations are reported in Supplementary Section S6. Images of structures and 2D maps/isosurfaces were produced with Diamond 49 and VESTA 50 software, respectively. Infrared and Raman spectra (reported in Supplementary Section S9) were calculated by means of the Coupled Perturbed Kohn-Sham approach implemented in CRYSTAL14 [51][52][53] .