Carbonic Acid Revisited: Vibrational Spectra, Energetics and the Possibility of Detecting an Elusive Molecule

We calculate harmonic frequencies of the three most abundant carbonic acid conformers. For this, different model chemistries are investigated with respect to their benefits and shortcomings. Based on these results we use perturbation theory to calculate anharmonic corrections at the {\omega}B97XD/aug-cc-pVXZ, X=D,T,Q, level of theory and compare them with recent experimental data and theoretical predictions. A discrete variable representation method is used to predict the large anharmonic contributions to the frequencies of the stretching vibrations in the hydrogen bonds in the carbonic acid dimer. Moreover, we re-investigate the energetics of the formation of the carbonic acid dimer from its constituents water and carbon dioxide using a high-level extrapolation method. We find that the {\omega}B97XD functional performs well in estimating the fundamental frequencies of the carbonic acid conformers. Concerning the reaction energetics, the accuracy of {\omega}B97XD is even comparable to the high-level extrapolation method. We discuss possibilities to detect carbonic acid in various natural environments such as Earth's and Martian atmospheres.

spectra of gaseous carbonic acid species trapped in a solid matrix of either neon or argon. 9 Obviously, these spectra are thought to represent a benchmark for possible identification of naturally occurring carbonic acid vapor.
Accompanying the experiments, a significant amount of theoretical and computational work has been carried out in order to investigate the chemical properties of carbonic acid and its stability. A few of them, in particular those which are important for the present work, are discussed in the following. In 1995, Wight et al. 10 investigated the potential energy surface and the harmonic vibrational spectrum of carbonic acid monomers at the MP2/6-31++G** level of theory. It was shown that the global minimum structure of the monomer is cis-cis, but that the cis-trans isomer is less stable by only 1.6 kcal/mol including a ZPE correction. The two isomers differ by rotation of one OH bond, see FIG. 1a and b. Furthermore, it has been shown that these carbonic acid monomers are only metastable with respect to dissociation to water and carbon dioxide by 10.4 kcal/mol, but reactants and products are separated by a high potential barrier of about 51.9 kcal/mol. Hence, the carbonic acid cis-cis and cis-trans monomers are kinetically stable, but thermodynamically they would decompose. Two years later, in 1997, Liedl et al. 11 studied the energetics of the carbonic acid dimer consisting of two cis-cis monomers. Such a cyclic dimer is stabilized by two strong hydrogen bonds, see FIG. 1c. They investigated whether the dimer is thermodynamically stable with respect to water and carbon dioxide. It was concluded that with the best then feasible method, i.e. MP2/aug-cc-pVTZ calculations, that 'it is not possible to judge […] if the zero-point-corrected energy for the formation of the carbonic acid dimer from water and carbon dioxide is slightly positive or slightly negative; however, it is astonishingly close to zero'. 11 Hence, in gaseous carbonic acid one can most probably expect a mixture of the cis-cis and cis-trans monomers and the carbonic acid dimer shown in FIG. 1. A study of Tossell 12 in 2006 on the H 2 CO 3 cis-cis and cis-trans monomers (and on oligomers composed of cis-cis monomers, using the CBSB7 B3LYP model chemistry and approximations for anharmonic frequencies at the CCSD(T)/6-311G+(2d,p) level revealed the importance of going beyond the harmonic frequency approximation for a reliable calculation of infrared spectra of carbonic acid. Particularily the deviations from the harmonic values for stretching and bending of bonds that include hydrogens can be quite large and this is even more the case for hydrogen bonds. Tossell found that the oligomers are more stable the larger they are. 12 In the gas phase one can, however, expect that oligomers are mainly restricted to the dimer species, since larger ones require more collision events and are thus thought to be very limited due to probability considerations. 9 Bernard et al. calculated the harmonic frequencies at the MP2/aug-cc-pVDZ level of theory to compare with their experimental results and obtained a good agreement between theoretical harmonic frequencies and experimental fundamental ones. 9 In fact, this raises the question of the reliability of this model chemistry for predicting accurate frequencies, as including anharmonic corrections might cause a considerable redshift of the line positions. The assignment of spectral lines in the work of Bernard et al. 9 is based on isotope shifts and not on direct comparison of frequencies, i.e. the assignments of vibrational modes 9 are still thought to be well based and only the reliable prediction of infrared spectra with MP2/aug-cc-pVDZ calculations is questioned. Murillo et al. 13 investigated the energetics of different dimer geometries and especially the effect of the geometry on the strength of hydrogen bonds.
They found the strongest hydrogen bonds to be associated with the dimer composed of one cis-cis and one cis-trans monomer, not for the one with two cis-cis units. The latter dimer, however, remains the energetically most stable one at the CCSD(T)/aug-cc-pVDZ//MP2/6-311++G** level of theory. In 2011, Wu et al. 14 calculated a variety of molecular properties of carbonic acid isomers at the CCSD(T)//B3LYP/aug-cc-pVTZ level of theory. They reported harmonic frequencies and investigated the potential energy surface of H 2 CO 3 . In agreement with earlier work 10 they found the cis-cis and cis-trans monomers to be the most stable ones, but lying 6.72 kcal/mol (8.28 kcal/mol for the cis-trans monomer) above water and carbon dioxide and thus being only kinetically stabilized by a barrier of about 41.58 kcal/mol (43.14 kcal/mol for the cis-trans monomer).
FIG. 1. Optimized geometries at the ωB97XD/aug-cc-pVQZ level of theory for (a) the cis-cis monomer, (b) the cis-trans monomer and (c) the energetically lowest dimer of carbonic acid. Bond lengths are given in Ångstrom and bond angles in degrees.
In this work we aim to improve on the aforementioned results, especially concerning the infrared absorption spectra for gas-phase carbonic acid species and the energetics of carbonic acid dimer (referring for the rest of this paper to the one composed of two cis-cis monomers). The quantum chemical methods that are used are described in section 2.
In section 3.1, we shall compare several commonly used methods in quantum chemistry to calculate harmonic frequencies. In section 3.2, we derive anharmonic corrections and discuss the basis set dependencies of the results.
Concerning the energetics of the carbonic acid dimer we aim for giving a definite answer on the question whether it is lower in energy than its constituents water and carbon dioxide in section 3.3. To this aim, we restrict our calculations to the three lowest-energy carbonic acid species depicted in FIG. 1, namely the cis-cis and cis-trans monomer and the dimer composed of two cis-cis monomers.
We use these theoretical considerations for the aim to tackle a possible detection of carbonic acid in natural environments from an observational point of view as well in section 3.4. Particularly, we discuss Earth's and Martian's atmospheres as regions of naturally occurring carbonic acid and what problems and limitations one has to face aiming for a successful detection of this elusive molecule. In section 4, we summarize the work.

A. Frequencies
We calculate harmonic frequencies with three standard methods of quantum chemistry. Two of these stem from density functional theory, namely the B3LYP functional, which uses Becke's three parameter hybrid exchange functional 15 and a correlation functional of Lee, Yang and Parr with both local and non-local terms, 16 and the ωB97XD functional of Head-Gordon et al. 17 including long range corrections as well as empirical dispersion corrections, both of which are not included in B3LYP. We also perform MP2 18 calculations on our systems. All methods require basis sets to expand the density or wave function. We investigate the influence of the finite basis set size on the resulting IR spectra with Dunning's correlation consistent polarized basis sets augmented with diffuse functions, 19 aug-cc-pVXZ, with X=D, T, Q for the B3LYP and ωB97XD functionals and X=D only (for reasons of computational cost) for the MP2 method.
Anharmonic frequencies have been calculated at the ωB97XD/aug-cc-pVXZ level of theory with X=D, T, Q for the monomeric species and X=D, T only (for the same reasons) for the carbonic acid dimer using the MP2 approach 20 as implemented in the Gaussian 09 software. 21,22 For the very anharmonic O-H stretching modes of the hydrogen bonds stabilizing the carbonic acid dimer another approach is necessary for an accurate estimation of the frequency shifts in the asymmetric and symmetric mode of the two hydrogens in the eight-membered ring of the dimer. The full vibrational problem of two anharmonic O···H-O oscillators is solved, i.e., the time-independent Schrödinger equation in two dimensions was solved quantummechanically. As details of an analogous procedure can be found in an earlier work 23 , we will only give a short overview in the following. In a first step a two-dimensional potential energy surface 12 ( , ) V q q , where 1 q and 2 q denote the coordinates of the two hydrogen atoms with respect to their oxygen neighbors, i.e. the O-H bond lengths, has been calculated at the ωB97XD/aug-cc-pVXZ level of theory with X=D, T, Q. The two bond lengths were varied between 0.6 and 2.0 Å in steps of 0.1 Å. Altogether 15 15 225  single point calculations were performed to construct the potential energy surface for one basis set. Different basis sets have been used to determine the variation of the frequency shifts with basis set size. In a second step, a spline-type function consisting of smooth fourth-order polynomials connecting the ab initio generated grid points was calculated. In a third step, the vibrational Hamiltonian of the oscillator system was expressed as

B. Energetics
The energetics of the carbonic acid species have been investigated by calculating ground state energies of reactants and products at various levels of theory which are then combined in a way similar to the Gx extrapolation methods 26,27 to approximate the energies related to the dissociation of the carbonic species to water and carbon dioxide. In a first step the energies have been calculated at the G4(MP2) 28 level of theory in order to have reference values. Geometries and zero-point energies were already available from the frequency calculations. In a second step, we calculated energies at the CCSD(T)/aug-cc-pVXZ//B3LYP/aug-cc-pVQZ, CCSD(T)/aug-cc-pVXZ//ωB97XD/aug-cc-pVQZ and CCSD(T)/aug-cc-pVXZ//MP2/aug-cc-pVTZ levels of theory 29 with X=D, T to investigate the effect of different optimization levels. The third step consists of the extrapolation of these energies to the CCSD(T)/aug-cc-pVQZ//(optimization level) levels of theory using Truhlar's extrapolation scheme 30  All electronic structure calculations have been carried out using the Gaussian 09 software. 32

A. Harmonic Frequencies
We compare the harmonic frequencies obtained with B3LYP, ωB97XD and MP2 in conjunction with the aug-cc-pVXZ, with X=D, T, Q, basis sets in order to investigate dependencies on basis set size. Furthermore, we compare the results to the recently determined experimental values of Bernard et al. 9 In the latter work, experimental results are given only for a few groups of vibrational modes making no distinction between monomer and dimer species. For our purpose of a comparison of methods it is sufficient to adopt this strategy. The harmonic frequencies are given in  Assignments to specific molecules via the molecular symmetry of the vibrational modes reported in the paper of Bernard et al. 9 c i.p. and o.p. denotes in-plane and out-of-plane, respectively. d CO3 out-of-plane bending modes are associated with the vibration of the carbon out of the plane of the three oxygen atoms.
Except for the O-H stretching modes, the difference between double-and triple-zeta basis sets is rather large (up to about 20 cm -1 ), whereas the differences between triple-and quadruple-zeta basis sets remains already much smaller (below about 5 cm -1 ). Concerning the O-H stretching modes, however, increasing the basis set from triple-to quadruple-zeta can still lead to substantial changes in the resulting harmonic frequencies (up to 15 cm -1 ), especially with the MP2 method. In that case, increasing the basis set size from double-zeta to quadruple-zeta changes the O-H stretching vibrations by 30 cm -1 , 25 cm -1 in case of O=H stretching modes, 20 cm -1 in case of C-OH stretching modes, 5 cm -1 in case of C-O-H in-plane bending modes and 20 cm -1 for the CO 3 out-of-plane bending modes. In the same order, the B3LYP frequencies differ by 7 cm -1 , 6 cm -1 , 5 cm -1 , 11 cm -1 and 10 cm -1 and in case of the ωB97XD functional by 15 cm -1 , 5 cm -1 , 6 cm -1 , 14 cm -1 , and 9 cm -1 for the maximal shift from double-zeta to the quadruple-zeta basis set size. Thus, the deviations with respect to basis set size are smaller in case of DFT than in case of MP2. The generally weaker dependence on basis set size of DFT methods is rather well-known. Note that for most of the vibrational modes increasing the basis set size upshifts the harmonic frequencies.
We conclude that at all three levels of theory, one should be careful about the accuracy of the predicted harmonic frequencies. The aug-cc-pVDZ basis set is too small, especially with MP2 calculations. Except for the O-H stretching modes, the aug-cc-pVTZ basis set yields results already close to the ones obtained with the quadruple-zeta basis. However, for the the O-H stretching modes it appears to be necessary to use even larger basis sets. In addition to these basis set dependences one observes that, except for the O-H stretching modes, the harmonic frequencies obtained with B3LYP and MP2 appear to be rather too low in absolute values. In various cases they are smaller than the experimental fundamentals. Although the general trend of increasing the basis set size is a blueshift of the harmonic frequencies, they are still very close to the experimental results, which are thought to be slightly blueshifted due to interactions between carbonic acid molecules and the surrounding neon matrix. 9 Considering from earlier work on carbonic and formic acids 12,33 that the fundamental frequencies of such systems are mostly redshifted compared to harmonic ones, we conclude from this rough inspection that the ωB97XD functional is the most reliable method amongst the three ones tested. Of course, this does not come as a surprise since it includes long-range and empirical dispersion corrections, and has been reported to perform especially well in systems exhibiting weak bonds such as hydrogen bonds. 17 For this reason, we restricted the investigation of fundamental frequencies to this functional.

B. Anharmonic Corrections
Harmonic and fundamental frequencies, frequency shifts and fundamental-harmonic frequency ratios obtained at the ωB97XD/aug-cc-pVXZ, X=D, T, Q levels of theory are summarized in tables 2 to 4 for the carbonic acid cis-cis monomer, the cis-trans monomer and the carbonic acid dimer, respectively. For reasons of computational cost, only double-and triple zeta basis sets were used in case of the dimer. Those vibrational modes for which experimental data is available from the paper of Bernard et al. 9 are also given in the respective tables. Tossell, 12 already calculated fundamental frequencies for the cis-cis dimer at the CCSD(T)/6-311+G(2d,p) level and corrected these for anharmonicity by adding to them the difference between B3LYP/CBSB7 34 anharmonic and harmonic frequencies, see also the work of Carbonniere et al. 35 Since these fundamental frequencies are believed to be the best available theoretical ones, we compare our results to them in  Concerning the carbonic acid cis-cis monomer, we observe a good agreement within about 5 cm trend of increasing frequencies with increasing basis set size can also be seen from both harmonics and fundamentals. Therefore, due to error cancellation, frequencies obtained with the smaller basis sets fit better to experimental ones in those cases where the ωB97XD method seems to overestimate the fundamental frequencies.
Furthermore, we note that an ad-hoc inclusion of anharmonicity by empirical scaling factors might lead to significant errors, especially in case of the H-O-C out-of-plane bending modes at 300-400 cm -1 , where the ratio of fundamental to harmonic frequency is as small as 0.5-0.7 while the usual scaling factors are 0.95-0.97. 15,36 The ratios associated with the other vibrational modes are distributed in the interval 0.93-0.98. In case of the carbonic acid cis-trans monomer the agreement between the ωB97XD/aug-cc-pVQZ  stretching modes in case of the double-zeta basis, and by -21 and -42 cm -1 in case of the triple-zeta basis. As in the case of the cis-trans monomer we observe that especially for hydrogen stretching modes B3LYP appears to be well suited to estimate fundamental frequencies, whereas the ωB97XD functional seems to converge from below to slightly overestimated fundamental frequencies.  9 For comparison, the values for solid cryogenic ß-H2CO3, which has a one dimensional H-bonded chain structure, are given in parenthesis and have been adopted from the paper of Kohl et al. 37 In the latter work, the authors report without assignment spectral lines at 605, 307, 258 and 193 cm -1 .

C. Stability and energetics of the dimer
We discuss the energetics of the formation reaction of the carbonic acid dimer The contrary is true for the MP2 method. The energy obtained with the double-zeta basis set is quite significantly positive with 5.59 kcal/mol, but is significantly reduced by going to the triple-zeta basis to 1.9 kcal/mol. Note that for estimating the zero-point energy in case of the MP2/aug-cc-pVTZ level of theory, the value of ΔZPE at the MP2/aug-cc-pVDZ level of theory has been taken based on the observation that ΔZPE does only depend very weakly on basis set size. 11,13 In this case, the counterpoise correction yields a significantly different result of 3.93 kcal/mol.
The reason for this is probably again that for MP2 the triple-zeta basis is still too small to yield converged results. It has even been argued that one should not rely too heavily on counterpoise corrections when estimating the formation energy of the carbonic acid dimer. 11  The accuracy of our best estimate can be considered as a significant improvement over, for example, the G4(MP2) extrapolation scheme for the following reasons: The optimization and frequency calculations use the aug-cc-pVQZ basis set, which yields 824 basis functions in case of the dimer, in contrast to the 6-31G(2df,p) basis set, 38 which yields 244 basis functions in case of the dimer. In addition, the ωB97XD functional is used, which appears to be especially suited for molecules containing weak Hbonds due to the inclusion of long-range corrections and empirical dispersion. Moreover, fundamental frequencies are estimated to incorporate anharmonic corrections into the zero-point energy instead of using the harmonic approximation. wB97XD/aug-cc-pVTZ -6. It should be noted that the relevant quantity if discussing equilibria is, of course, the free energy of reaction, which is only partly given by the reaction energies discussed above. It should be kept in mind that the free energy of dimerization is more positive than the reaction energy due to the reduction of two molecules to one and is expected to be strongly temperature dependent.

D. Detection of Carbonic Acid
Water and carbon dioxide ices exist in various different environments in outer space. Carbonic acid is expected to be present in these environments as well. 5,6,8,9,14,39 It has been suggested that solid carbonic acid may be found on the Martian surface, on interstellar grains or on Jupiter's satellites Europa, Ganymede, and Callisto. Gaseous carbonic acid is expected to exist in our own atmosphere, 5 in the atmosphere of Mars and Venus and even on comets. 5,9 So far, however, no spectroscopic detections of carbonic acid have been reported in any of these environments.
In the following we discuss the reasons for this and the possibilities to detect carbonic acid via spectroscopic observations in the gas phase in the two most promising environments, namely Earth's atmosphere and Mars.
In Earth's atmosphere it is thought that carbonic acid may be formed in the upper troposphere, most likely in cirrus clouds. 5 Cirrus clouds typically occur in heights between 8 and 12 km, where one finds temperatures around 250 K and a pressure around 250 mbar. Following Hage et al., 5 carbonic acid can be sublimated and recondensed without decomposition into CO 2 and H 2 O. The temperature and pressure range relevant for the sublimation in the laboratory has been 180 to 220 K and 10 -7 to 10 -5 mbar. In our atmosphere we find regions with this temperature in a height of about 15 to 20 km, which is already above the highest cirrus clouds. Looking for a region with a pressure around 10 -5 mbar, we have already to go up to heights over 120 km. Although we expect carbonic acid to be found even at higher temperatures and pressure, like in the cirrus clouds, the collision rate is expected to be much higher leading to a much lower abundance of carbonic acid. Thus, it might be very challenging to detect carbonic acid in our atmosphere. Additionally, the detection is complicated by the fact that the spectroscopic lines of carbonic acid are not in a preferable part of the spectrum. We compared our model spectra with the atmospheric spectrum including all As the Atacama desert is known as one of the driest places we refined this equatorial profile by additionally selecting appropriate data from the Global Data Assimilation System (GDAS 44 provided by the National Oceanic and Atmospheric Administration, NOAA 45 ). These meteorological data are available on a 3 hour basis, as a global, 1 degree latitude/longitude data set (from December 2004 to present and ongoing), and contain height profiles for the temperature, pressure and relative humidity. In addition, the on-site meteo monitor of ESO was also incorporated to finally create an averaged atmospheric profile appropriate for Cerro Paranal. Unfortunately the spectral lines of carbonic acid are in a range where gluts of lines originating from other molecules are certain to be found. We were especially searching for the spectral feature of the O-H stretch of the carbonic acid dimer, which is expected to be the most prominent one. There are a few windows in the atmospheric spectrum, where the atmosphere is extensively transparent. There is no intermediate absorption and re-emission. In the range from about 4000 cm -1 to higher wavenumbers there is the astronomical K band (FIG. 2, right part). In the K band there is hardly any continuum emission, therefore the background is very small. The range from about 3500 cm -1 to lower wavenumbers is the L band, which has already more continuum emission compared to the K band, but at least there are a few emission-free km. The pressure in these heights is below 1 mbar at 20 km and even lowers to about 10 -3 mbar at 60 km. The temperature range is between 100 and 200 K depending on height and solar irradiation. Thus, the conditions for the detection of carbonic acid are promising.
In the following, we discuss the possibility of observations from ground based telescopes or satellites. Direct measurements of the atmosphere with the Martian surface in the background seem impossible because of the selfradiation of the planet. The spectroscopic measurements would be completely overlaid. For this reason, the possibility of grazing observations was taken into account. The thickness of the atmosphere that has to be observed is about 50 km. Taking the smallest possible distance Mars -Earth, which is about 56 10 6 km, leads to a required resolution of at least 0.1 arcsec. E.g., the Hubble Space Telescope (HST) has this resolution in the optical part, the resolution in the required infrared is worse. With all present instruments it is not possible to achieve the required resolution. The next generation space telescope James Webb Space Telescope (JWST) will be in the same range as HST, so it seems we have to wait for the E-ELT, the European Extremely Large Telescope. With this telescope we will reach a resolution below 1 milliarcsec. First light for the E-ELT is expected in 2022. 46 Another option is space exploration missions. The planned mission MAVEN (Mars Atmosphere and Volatile EvolutioN) will be launched in late 2013. 47 The aim of the mission is collecting data about the Martian atmosphere.
During this mission a measurement and detection of carbonic acid seems to be feasible.
In order to search for carbonic acid in gas phase in other astrophysical environments like the interstellar matter (ISM) or Jupiter's satellites, further studies should concentrate on providing high-resolution infrared spectra from laboratory measurements and investigations on temperature and pressure effects. These data are essential for extended future spectroscopic investigations in outer space.

IV. Conclusions
We find that, due to cancellation of errors, the B3LYP and MP2 methods lead to harmonic frequencies which are close to experimentally determined fundamental ones.. A more systematic prediction of fundamental frequencies based on the ωB97XD functional led to good agreement with experimental results for the carbonic acid cis-cis monomer, but only to a rough estimate for the cis-trans monomer and the carbonic acid dimer under consideration.
Nevertheless we observe a slightly better agreement with experiment than in earlier theoretical predictions where the B3LYP functional was employed for most vibrational modes of the carbonic acid dimer. B3LYP appears to perform very well for O-H stretching modes in the region around 3600 cm -1 and H-bond vibrations at 2600-2800 cm -1 . For this reason, both methods seem to provide a qualitatively good first estimate of fundamental frequencies. Basis sets of triple-zeta quality or larger are necessary in any case. Obviously, the estimation of accurate fundamental frequencies remains a difficult task and if aiming for reliable quantitative predictions, one has to go to significantly higher levels of theory. This has been shown in the comparison of theoretical and experimental data in case of the cis-trans monomer. Though the results of Tossell 12 were obtained at a high level of theory and reproduce the experimental data quite nicely, even they still exhibit differences of up to 40 cm -1 .
A re-investigation of the energetics of the formation of the carbonic acid dimer from water and carbon dioxide has affirmed that carbonic acid is not only experimentally a very elusive molecule. It turned out that a quite high level of theory is needed to draw a definite conclusion on the 15 years old question 11 whether the carbonic acid dimer is more stable than its constituents water and carbon dioxide. At least according to our best estimate this question can be answered with yes. Furthermore, we note that the clearly negative reaction enthalpy of about -5.14 kcal/mol at zero temperature can be obtained quite accurately at the ωB97XD/aug-cc-pVQZ level of theory (leading to -5.48 kcal/mol) at much reduced computational cost compared to CCSD(T) and related methods. This confirms the suitability of that functional, especially for molecules containing weak bonds such as hydrogen bonds.
Despite its presence is very likely, the detection of carbonic acid in astrophysical environments seems quite challenging. With present spectroscopic methods and instruments it will not be possible to detect carbonic acid in Earth's atmosphere. An upcoming space exploration mission makes it imaginable to detect it in the Martian atmosphere.