The Crystal Structure of Carbonic Acid

: Ubiquitous carbonic acid, H 2 CO 3 , a key molecule in biochemistry, geochemistry, and also extraterrestrial chemistry, is known from a plethora of physicochemical studies. Its crystal structure has now been determined from neutron-diffraction data on a deuterated sample in a specially built hybrid clamped cell. At 1.85 GPa, D 2 CO 3 crystallizes in the monoclinic space 15 group P 2 1 / c with a = 5.392(2), b = 6.661(4), c = 5.690(1) Å, β = 92.66(3)°, Z = 4, with one symmetry-inequivalent anti - anti shaped D 2 CO 3 molecule forming dimers, as previously predicted. Quantum chemistry evidences π bonding within the CO 3 molecular core, very strong hydrogen bonding between the molecules, and a massive influence of the crystal field on all bonds; phonon calculations emphasize the locality of the vibrations, being rather insensitive to 20 the extended structure. One-Sentence Summary: carbonic acid and its crystal structure firmly established by neutron-diffraction

Carbonic acid, H2CO3, is the simple hydrate of carbon dioxide and a ubiquitous molecule whose name has been mistakingly used for carbon dioxide itself (1) at the end of the 19 th century; the H2CO3 molecule must have been touched by the alchemists already. As given in any standard chemistry textbook (2,3), carbonic acid forms directly from water and carbon dioxide according to the equilibrium equation 5 H2O(l) + CO2(g) ⇄ H2CO3(aq) ⇄ H + (aq) + HCO3 -(aq) ⇄ 2 H + (aq) + CO3 2- with the two acid constants of KA1 = 1.3 · 10 -4 and KA2 = 4.84 · 10 -11 . Because H2CO3 not only behaves as a weak acid but seemingly decays quite easily into the two starting molecules and, hence, makes bubbles in everyday soda water, preparing pure carbonic acid has been falsely considered plainly impossible for many decades until a pioneering low-temperature direct 10 synthesis according to Na2CO3 + 2 HCl → H2CO3 + 2 NaCl in dimethyl ether was successfully carried out (4) on the gram scale in 1965 already, including an estimate of its exothermic formation enthalpy (-649±42 kJ mol -1 ). The astonishing kinetic stability of carbonic acid has also been demonstrated from both experiment and quantum 15 chemistry (5), including the crucial role of excess water molecules in catalyzing the H2CO3 decay, essentially an autocatalytic process.
It does not come as a surprise that characterizing carbonic acid has therefore been mainly performed by spectroscopic methods, in particular as regards infrared data (6,7), and even two different polymorphs (α and β) of carbonic acid were claimed based on matrix IR techniques 20 (8,9). Only recently, however, "α-H2CO3" was characterized as being a monoethyl ester of H2CO3 while β-H2CO3 appears to be the only form of carbonic acid known up to the present day (10).
Despite the fundamental character of H2CO3 and its (biochemical) paramount importance for the regulation of the blood's pH, the acidification of the world oceans, its predicted chemistry in the 25 outer planets, and a true plethora of other issues (see references above and further references within), its crystal structure is still unknown, seemingly mirroring tremendous experimental challenges. Hence, the shape of the molecule in the crystal, the intermolecular connectivity, the details of the hydrogen bonds, the entire solid-state energetics and other related chemical and also physical questions are still unsolved, at least experimentally. With regard to theory, 30 however, there has been a groundbreaking density-functional prediction for H2CO3 having orthorhombic Pnma symmetry at about 1 GPa pressure (and beyond) based on evolutionary algorithms (11), in harmony with high-pressure spectroscopic data (12). This has been followed by the proposal of a unit cell for high-pressure water-CO2 phases (13) being of triclinic, not orthorhombic symmetry, however. 35 To solve the puzzle, in particular the necessity of correctly localizing the hydrogen atoms, a neutron-diffraction study on deuterated carbonic acid, D2CO3, seemed like an optimum choice but this would require a comparatively large sample volume of about 0.4 ml for GPa conditions not being commercially available; earlier unpublished calculations of ours assuming approximate isotypism of H2CO3 with urea firmly indicated stability slightly below 2 GPa (14). Hence, we 40 started with almost equimolar solid mixtures of crystalline H2O and CO2 manipulated at liquidnitrogen temperature (15) for orienting tests using tiny amounts with a standard diamond anvil cell and laboratory X-ray diffraction. For the decisive neutron-diffraction experiment, a much larger D2O/CO2 mixture with an intentional CO2 excess was ground and inserted into a unique self-built pre-cooled hybrid pressure cell (16) which was closed and pressurized with a load of 4 tons. After returning to room temperature the cell was further pressurized to a final pressure of 1.85 GPa as measured in situ by ruby fluorescence (16). Neutron powder diffraction patterns were collected at the STRESS-SPEC instrument at the FRM II neutron source in Garching (17) using two wavelengths. Three phases (D2CO3, CO2, and the Ni-Cr-Al fcc alloy of the container) 5 were simultaneously taken into account for the neutron-data analysis.
Manual structure solution inspired by the previous structure proposal (11) and further Rietveld refinements were carried out (18) but despite qualitative similarity with the predicted Pnma structure, a reduction to the translationengleich subgroup P1121/a (= P21/c) was needed which eventually yielded unit cell parameters as given in Table S1. To save parameters, atomic 10 distances and angles were softly restrained, with all atoms also restrained to be in-plane. Likewise, due to the intended presence of excess carbon dioxide in the form of CO2-I, its pattern served as an internal standard for the isotropic atomic displacement parameters while the Uiso parameters of the deuterium atoms were based on the D2O-VIII neutron study (19). All refined atomic parameters are found in Table S2. An overview of the final Rietveld result is provided in 15 Fig. S1. The ruby-determined pressure of 1.85 GPa was also independently verified using the measured molar volume of CO2-I (20-22), yielding perfect numerical agreement.
As regards the crystal structure itself, the unit cell of carbonic acid contains just one symmetryindependent D2CO3 molecule adopting the anti-anti conformation and four molecular units in total, see Fig. 1 25 counts, both bond-length irregularities must be taken seriously and directly point towards other mesomeric structures in the spirit of valence-bond theory, see Fig. 2, involving delocalized π bonding in the central CO3 triangle through well-overlapping C 2pz and O 2pz atomic orbitals. Hence, it is not surprising that even with a dramatically reduced restraining weight factor, the atoms within D2CO3 are in an almost perfect planar arrangement, the π bonding being in control 30 of the molecular shape. We note that this phenomenon has been found before in related molecular units, for example, in guanidinate salts involving a central CN3 core such as in SrC(NH)3 (23). would be classified as strong in terms of distance and medium in terms of angle (25), the angle being less reliable in this case, however (see below). The dimeric structure seemingly resembles the most stable molecular arrangement of carbonic acid already theoretically predicted in 1997 (26) and further theoretically explored recently (27,28). At the same time, this very dimer shape makes alternative D split positions with half occupancy less likely, as also corroborated 20 experimentally from significantly worse neutron-data refinements.
Because of the restricted resolution, it is instructive and also advantageous to compare the experimental structural data with those theoretical ones derived from DFT optimization using plane waves, PAW pseudopotentials, and a DFT-D3(BJ) or HSE exchange-correlation functional, also using P21/c space-group symmetry (29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39). All geometrical data are found in Table S3. While this DFT-optimized experimental structure is mechanically stable in terms of phonons (40) showing no imaginary frequencies (Fig. S2), also in an NPT molecular-dynamics (MD) simulation (41), see Fig. S3, the C-O1/O2 bonds arrive at an even shorter 1.30 Å instead of the expected single-bond 1.43 Å whereas O-D is confirmed as 1.02 Å. The theoretical 5 O1-O3* distance and O1-D1···O3* angle (as well as O2-D2···O3* because the H-bonding mode gets slightly more symmetric) widen to 2.57 Å and 176°. Most significant, however, is the C=O3 bond which then comes out shorter (1.26 Å) but not as short as in the prototype double bond. As such, π bonding over the entire CO3 triangle is nicely corroborated from DFT, as analyzed in more detail below. In what follows, we will refer to these theoretical bond distances 10 because they look superior to the experimental ones while still capturing all the experimentally determined structural motifs. Before doing so, a look at Fig. 3 (left) shows the packing of the dimers which essentially form infinite stripes along the b direction.  (Table S3). On the right, we provide an alternative, energetically even lower polymorph found by density-functional theory in which the same stripes run along a but this time in a zig-zag stacking mode (Table S4).
The theoretical free energy of formation of D2CO3 at 2 GPa including zero-point energies arrives 20 at −684 kJ mol −1 at 300 K and slightly lowers to −687 kJ mol −1 as the temperature cools down to 0 K, assuming a direct reaction from solid D2 in P63/mmc (42), diamond-like C (43), and solid O2 in Cmcm (43). This theoretical energy exceeds the 1965 experimental estimate by only 6% and is within a single standard deviation of that value.
There are two more puzzling theoretical findings which relate to the existence of hypothetic 25 molecular polymorphs. First, the theoretically predicted crystal structure of H2CO3 in Pnma space-group symmetry (11), also found here by means of structural searches via particle swarm optimization (44,45), can be firmly ruled out given the experimental data but it nonetheless turns out as being electronically lower in energy, stabilized by about 5 kJ mol -1 as calculated by PBE-D3, surprisingly enough. While this small energy difference clearly challenges density-functional 30 theory by its minute size, it is at least thinkable that the Pnma polymorph may exist and can be made using other growing conditions, yet to be found, but only if sufficient activation barriers exist between the polymorphs. Likewise, there is yet another monoclinic P21/c polymorph including the same D2CO3 molecule with identical interatomic distances and O-D···O* angle; its data have been compiled in Table S4. This hypothetical polymorph, also stable in terms of phonons and NPT MD simulations (see Figs. S3 and S4), turns out as being even lower by 1.5 kJ mol -1 (PBE-D3) and by almost 1 kJ mol -1 in terms of Gibbs free energies; it exhibits identical dimers and practically the same connectivity leading to molecular stripes but this time allowing for another stacking, running along the a direction as shown in Fig. 3 (right). Further experimental and theoretical studies may evidence whether those two hypothetical polymorphs 5 go back to theoretical deficiencies (say, inaccurate dispersion correction and/or neglect of nuclear anharmonic vibrational contributions) or must be considered real; for an overview of relative energies including various functional approaches, corrections as well as (negligible) H/D mass impact by phonon contributions, see Fig. S2. If the here found experimental structure were only metastable, Ostwald's step rule would be at play and render it as the polymorph with the 10 smallest energy barrier for crystallization (46).
In terms of chemical bonding, an analytical unitary transformation from the delocalized planewave picture into a localized representation was carried out (47-51). The various bonds are easily evaluated by their integrated crystal-orbital bond indices (ICOBI) given in Table S5 (52). As anticipated already from the mesomeric structures in Fig. 2 and the bond lengths themselves, 15 the "single" bonds (ICOBI ≈ 1.21, not 1) come out stronger, and the "double" bond (ICOBI ≈ 1.39, not 2) is weaker. Indeed, the π bonding originating from the pz-pz interaction alone amounts to a bond order of about ⅓, in harmony with MO theory (23), see also below. The O-D bonds (ICOBI ≈ 0.66) appear as slightly weaker than expected, a clear consequence of the astonishingly strong hydrogen-bridging bonds (ICOBI ≈ 0.18) into the opposite direction. Upon 20 looking at the isolated monomer (also Table S5), the "single" bonds get closer to single (ICOBI ≈ 1.08) and the "double" bond gets closer to double bond character (ICOBI ≈ 1.72) while the isolated O-D hydroxyl bonds (devoid of H bonding) also become more typical (ICOBI ≈ 0.79). Clearly, the bonding of carbonic acid is heavily affected by the periodic structure since the crystal field equilibrates the C-O and O-D bond orders. 25 Energywise, similar conclusions can be drawn from Crystal Orbital Hamilton Population (COHP) plots depicted in Fig. 4 whose numerical data are found in Table S6. The strongest C-O3 "double" bond exclusively exhibits bonding levels over the entire energy region, and the highest π part, between -3 eV and the Fermi level, shows up antibondingly for the "single" C-O1 and C-O2 bonds. A small antibonding level of the hydroxyl O1-D1 and O2-D2 single 30 bonds in the proximity of the Fermi level reappears bondingly in the hydrogen bridging bonds O3···D1* and O3···D2* which come out exceedingly strong, with ICOHP values of about -1.84 eV (Table S5) corresponding to ca. 178 kJ mol -1 . A while ago we have shown how covalency contributes to hydrogen bonding, measured by ICOHP (53), and the current value is even larger than for the H bond in the DFT-optimized benzoylacetone molecule (53). While 35 ICOHP does not exactly equal the covalent bond energy but measures the covalent part of the band-structure energy, it nonetheless serves as a reliable numerical criterion, even when densitybased schemes (say, bond-critical points etc.) may fail. Clearly, this hydrogen bond is extraordinarily strong. The π bonding already alluded to in Fig. 4 is easily visualized by transforming those bonding interactions between -3 eV and the Fermi level from reciprocal into real space using a molecular-orbital based representation depicted in Fig. S5 (54,55). In the crystal, one finds three visually identical π orbitals of mostly bonding nature, as expected, localized on the three carbon-10 oxygen bonds, in harmony with the three mesomeric structures depicted in Fig. 2.
In terms of a general comparison of the covalent bonding energetics in carbonic acid in the crystal and in the isolated monomer (without any H bonding), the numerical entries (Table S6) provide an unexpected finding. When the monomer condenses into the crystal, all intramolecular bonds strengthen quite a lot energywise, in particular the C-O "single" bonds (ca. +30%) but 15 also the C=O "double" bond (ca. +2%) and also the O-D hydroxyl bonds (ca. +4%). As anticipated from the aforementioned bond-order data, the crystal field also affects the intramolecular bond energies, in this case by stabilizing everything collectively.
As regards ionic contributions, the wavefunction-based Löwdin charges as projected from plane waves are given in Table S7 and arrive at +0.64 for C, -0.42 for the hydroxyl O atoms, -0.50 for 20 the carbonyl O atom, and +0.35 for hydrogen; for the isolated molecule not profiting from a periodic Madelung field, the charges must be somewhat higher (and they are) to establish about the same Coulomb forces. In any case, either solid-state or molecular, ionicity is also at play, besides the dominating covalency. For the crystalline solid, ionicity straightforwardly leads to the aforementioned Madelung field and a (electrostatically defined) lattice energy for crystalline 25 carbonic acid. This electrostatic part of the lattice energy is easily calculated based on Löwdin charges (52) and arrives at about -980 kJ mol -1 , a truly substantial amount for a molecular crystal, on the order of rock salt, NaCl. Admittedly, such point-charge model for a molecular solid must be considered with some caution. Interestingly enough, the alternative lowest-energy polymorph proposed by DFT has a 1.9% larger lattice energy, so its stabilization does not only go back to better packing but reflects a collective ionicity-based electrostatic effect.
If we focus on an isolated dimer of carbonic acid in vacuo, the transition state for the concerted proton transfer was localized from ab initio molecular-orbital calculations (56, 33,57). The associated barrier arrives at +6.5 kJ mol -1 , which would result in rapid conversion at room 5 temperature. In the extended solid, however, there is a stabilization of the ground state or a destabilization of the transition state by the surrounding periodic units, that is, by the aforementioned Madelung field. Attempts localizing the product resulting from such concerted proton transfer on the periodic potential-energy surface failed, indicating that it is indeed destabilized to the point of not appearing as a stationary point on the latter, in perfect harmony 10 with the neutron data.
It is rather instructive to compare the calculated IR and Raman spectra of carbonic acid with those that have been published before. Fig. 5 shows such a comparison with the experimental results (58) at ambient pressure indicating moderate (IR) to good (Raman) agreement, in particular for the experimental DFT-optimized structure; residual deviations may be related to 15 the neglect of anharmonicity, temperature, in addition to other factors, as previously reported (11,59). On the other hand, the calculations clearly evidence that the positions of the IR and Raman signals are almost unaffected by the crystal symmetry, thereby unsurprisingly highlighting the local nature of the molecular vibrations. It is somewhat surprising, however, that the match between theory and experiment does not become better if one analyzes the high-20 pressure experimental data (12) depicted in Fig. S6, in particular because these lack some of the important Raman signals, probably due to experimental difficulties. Nonetheless, the influence of the pressure on the peak positions must be considered weak, simply by visually comparing Fig. 5 with Fig. S6.  30 The crystal structure of carbonic acid is fundamental in character and also unique, as other simple species for potential comparison are rather rare. One exception would be given by formic acid, H2CO2 = HCOOH, also exhibiting symmetric double hydrogen bonding in all known physical phases (60), namely gas, liquid, solid, matrix, and molecular beam, and its enthalpy of dimerization by two H bonds is 150.7±0.8 kJ mol -1 . A probably even more fitting species to compare to is given by trithio carbonic acid, H2CS3, thus including the higher homologue of oxygen. That molecule likewise crystallizes at lower temperature, but this time without external pressure, hence it was structurally clarified decades ago already (61,62). In that pioneering study 5 one finds planar SC(SH)2 molecules in the same P21/c space group with mean C=S = 1.65 and C-S = 1.73 Å distances, also linked through hydrogen-bridging bonds but leading to a helical arrangement of chains along the a direction and three-dimensionally connected by weaker contacts. We also note that this single-crystal X-ray structure analysis of H2CS3 made the authors recognize the C-S "single" bonds being shorter than the normal value (1.78 Å) and the C=S 10 "double" bonds longer than anticipated (1.61 Å) such that significant π bonding in the CS3 triangle originating from C 2pz and S 3pz orbitals is indicated, a simple consequence of the molecular shape.
Summarizing, the monoclinic crystal structure of carbonic acid has been elucidated from a highpressure neutron-diffraction study of deuterated D2CO3 and further analyzed from density-15 functional theory. One finds anti-anti shaped dimers with a symmetric double H bonding motif, the latter bonds being extraordinarily strong, and this structural motif further leads to molecular stripes along the monoclinic b axis; other polymorphs might be possible but there is no experimental evidence as yet. Careful inspection of the molecular shape reflects substantial π bonding in the central CO3 core, and a detailed comparison between molecule and periodic 20 structure evidences a massive influence of the internal crystal (Madelung) field affecting bond lengths, bond orders, bond energies, and also atomic charges. Hence, besides being a rather trivial molecular crystal, solid-state carbonic acid also resembles a typical inorganic salt. Prior experimental vibrational data are in moderate (IR) or good (Raman) accord with the theoretical simulations but theory evidences their local character and their insensitivity to the given 25 pressure. By means of this fundamental small-molecule structural study, we hope to catalyze further high-pressure experiments as regards the chemical reactivity of carbonic acid.

This PDF file includes:
Materials and Methods Figs. S1 to S6 Tables S1 to S7

Other Supplementary Materials for this manuscript include the following:
Data S1 and S2

Materials and Methods
Neutron powder diffraction patterns were collected at the STRESS-SPEC instrument at the FRM II neutron source in Garching (17) using two wavelengths, 1.351 and 2.118 Å, and two 2θ ranges of 22.441-75.543° and 22.387-75.599°, respectively. To avoid reflections caused by the pressure apparatus, a radial collimator (field view = 2 mm) was used. Despite the (111), (200), and (220) reflections of the inner-cylinder alloy being visible at 38.25°, 44.45°, and 64.67° in the 1.351 Å diffraction pattern, no other pressure-cell signals were recorded. Three phases (D2CO3, CO2, and the Ni-Cr-Al fcc alloy of the container) were simultaneously taken into account for the neutrondata analysis.
Manual structure solution inspired by the previous structure proposal (11) in space group Pnma (reproduced by our structure prediction) and further Rietveld refinement were carried out using the GSAS-II program system (18). Despite some qualitative similarity of the data with the predicted Pnma structure for the strongest reflections around d = 2.35-2.6 Å, orthorhombic symmetry was clearly missing, so a reduction to the translationengleich subgroup P1121/a (= P21/c) was needed which eventually yielded unit cell parameters of a = 5.392 (2)  reasonable values of 0 ′ , 0 and 0 , the molecular volume of CO2-I was extrapolated to ambient pressure using experimental data (22). As such, the molecular volume of 39.633 Å³ is equivalent to a calculated pressure of 1.85 GPa, in perfect numerical agreement with the pressure determined from ruby fluorescence.
Phonon calculations were conducted using the supercell approach in the PHONOPY code (40), with calculational parameters as in the static calculation. The phonon calculations also served to calculate the Born charges and dielectric constants for simulating the IR and Raman spectra (41). Molecular dynamics (MD) simulations were carried out using NPT ensembles with a time step of 1 fs as implemented in MedeA VASP code (31,32,36). A plane-wave energy cutoff of 500 eV and a single Γ point were used. The simulations were performed in 2×2×2 supercells containing 192 atoms for all structures at 2 GPa and 298 K.
Structural searches were conducted through a particle swarm optimization algorithm implemented in the CALYPSO code (44,45) resting on density-functional theory (see above). To ensure the diversity and reliability of the predicted structures, structure simulations were performed with models containing up to four formula units at 1, 2, 10, 20, 30, and 40 GPa, respectively.
The periodic single-point delocalized wavefunctions were unitarily transformed into a complete set of atomic-orbital like wavefunctions of contracted Slater type by means of LOBSTER 4.1.0 (47-51). The results were further analyzed by calculating Crystal Orbital Bond Indices (COBI), projected Crystal Orbital Hamilton Populations (COHP), wavefunction-based Löwdin charges, and automatic Madelung energies based on the latter. Molecular bonding orbitals were generated from the plane-wave data using an accurate top-down frozen-density embedding approach (54) and a subsequent localization following the Pipek-Mezey algorithm (55), implemented into a development version of LOBSTER.
Molecular optimizations of ground-state and transition-state structures were done using Turbomole 7.5.1 (56) employing a PBE0-D3 hybrid density-functional approach (33) based on the def2-QZVP basis set (57). Stationary points were verified by frequency analysis.         Table S7. Atomic charges of carbonic acid in the crystal and molecule. Wavefunction-based Löwdin atomic charges of carbonic acid extracted from plane waves for the experimentally determined crystal structure, the DFT-optimized experimental crystal structure (in bold face, see text), the lowest-energy theoretical crystal structure, and the isolated monomer (in bold face, see text).