The Raman jet spectrum of trans-formic acid and its deuterated isotopologs: Combining theory and experiment to extend the vibrational database

Since its first spectroscopic detection in 1938,¹ the internal dynamics of formic acid have continuously been motivating experimental^2–12 and theoretical^13–18 investigations. The number of experimentally determined (ro-)vibrational parameters of formic acid in its ground and vibrationally excited states is enormous (see Refs. 6–9, 19, and 20 and references therein), the experimental data further extending to its ¹³C,^21–28 18O,^3,29,30 and ²H isotopologs.^19,31–37 One of the reasons for this interest results from formic acid having two low-lying conformational isomers, depicted in Fig. 1. They are connected through the large-amplitude O–H torsion, which is responsible for many interesting and challenging dynamical features that we will touch upon in this work. It is well known that this motion is not appropriately treated using rectilinear normal coordinates, as the torsion corresponds to atomic displacements perpendicular to the figure plane in straight line paths. To correctly investigate this degree of freedom and its coupling to the remaining degrees of freedom, one has to treat this motion as a true internal torsional coordinate.^38,39 From the perspective of this paper, the large anharmonicity of the torsion causes this vibration to tune in and out of resonance interactions as it is excited to increasingly high vibrational levels. This tuning provides an opportunity to interrogate a wide range of anharmonic couplings by matching theory with experiment. As an example, several experimental¹² and computational^16,17 studies have added clarity to the long-debated question of the assignment of the O–H in-plane bend ν₅ and O–H torsional overtone 2ν₉, which are in resonance. The higher IR intensity of the overtone has obscured the assignment of this resonance pair in the past as the more intense band was believed to be the fundamental.^8,40 Recently published Raman jet spectra intuitively support the recent proposal,⁴⁰ i.e., to assign the fundamental to 1306 cm⁻¹ and the overtone to 1220 cm⁻¹, where the band integral of the fundamental is seven times higher compared to the overtone.

This reassignment was triggered by recent high level anharmonic vibrational calculations, namely, vibrational configuration interaction using an internal coordinate path Hamiltonian (ICPH)¹⁶ and multi-configuration time-dependent Hartree (MCTDH),¹⁷ further providing the community with two full-dimensional CCSD(T)-F12 (explicitly correlated coupled-cluster singles, doubles, and perturbative triples) quality potential energy surfaces of cis- and trans-formic acid. While the reported MCTDH and ICPH results for trans-HCOOH agree well with each other and available experimental reference data, deviations for cis are significant with root-mean-square deviations of 14 cm⁻¹ and 49 cm⁻¹ for one- and two-quantum states, respectively.^16,17 Richter and Carbonnière consulted VPT2 energies for the nine fundamentals using both surfaces to rule out significant contributions from both potentials to the above mentioned discrepancies.¹⁷ 

In this work, we use these two surfaces to make detailed comparisons with experiment for several isotopologs of both isomers. Our comparisons will detail the differences in the harmonic and anharmonic contributions, going beyond VPT2 and fundamentals. To solve for the vibrational eigenstates, we use high order canonical Van Vleck perturbation theory (CVPT) as implemented in curvilinear coordinates. We shall see that this approach provides accurate energies, allows us to describe couplings with relatively small matrices, and provides a test of perturbative methods.

To benchmark both potentials, we can build on a wide range of perturbation-free vibrational reference data for trans-HCOOH extending far beyond fundamentals, and particularly noteworthy is the enormous list of IR bands reported by Freytes et al.⁷ Comprehensive Raman spectra of formic acid and its three deuterated isotopologs were published four decades ago by Bertie et al.^4,41 As their significant work was concerned with the formic acid dimer, the assignment of non-fundamental monomeric bands from their spectra is impaired by overlapping dimer contributions.⁴⁰ Recently published Raman jet spectra, which build on earlier studies,⁴² are optimized for monomer signals and as an additional feature of the heatable setup that is used,⁴³ peaks due to cold monomer, clusters, or hot transitions can easily be distinguished.¹¹ These Raman spectra revealed a richness of hot and (weak) combination and overtone bands of trans-formic acid, many of which have not been reported in the previous literature, particularly for the deuterated isotopologs. These data that are, in part, published^40,44 could not be fully analyzed and interpreted so far using available methods, such as standard VPT2. Particularly promising is the interpretation of hot bands as they are direct probes of anharmonicities. Naturally, the gas phase spectroscopy of cis-formic acid is more complicated due to the low abundance (one per mill) at room temperature.^25,45 The cis database²⁵ was significantly extended in recent years using the heatable Raman jet setup.^11,40,46 Our goal is to unify and summarize recent developments regarding the theoretical and experimental vibrational spectrum of trans-formic acid and its deuterated isotopologs below 4000 cm⁻¹.

This paper is structured as follows: We begin with discussing newly reported CVPT results, examining convergence, and comparing fundamental transition energies of both conformers to computational and experimental reference data. The main body of discussion is then centered around the assignment of new trans-formic acid peaks from Raman jet spectra and reassignment of literature IR bands. Finally, the two formic acid potential energy surfaces are compared.

Using the same experimental protocol as before,⁴⁰ we have recorded new spectra to fill spectral gaps up to 3750 cm⁻¹. For experimental details, we therefore refer to Ref. 40. In short, different isotopologs of formic acid were seeded into helium and expanded at different temperatures through a vertical slit nozzle into an evacuated jet chamber. The expansion was probed with a continuous-wave 532 nm laser, which was operated at 20 W for the new measurements. The scattered light was collected perpendicular to the laser and nozzle flow and focused onto a monochromator, which disperses the photons onto a liquid nitrogen-cooled CCD-camera (1340 × 400 pixels) binning over 400 vertical pixels. To account for the final resolution from the combination of laser and monochromator, we generously assign band center errors of ±2 cm⁻¹. Spectra of a temperature series are intensity-scaled to a trans-formic acid fundamental in each spectral region. The advantage is that hot bands, i.e., transitions from thermally populated vibrational states, can be easily distinguished from the cold monomer and cluster bands as they increase in intensity with rising temperature, whereas the cluster bands decrease (for further reading, see Refs. 11 and 46). Overtone and combination bands of the trans-formic acid monomer can furthermore be distinguished as their intensity remains constant over a scaled temperature series, similar to fundamentals.

The two potential energy surfaces we use in this work are referred to as PES-2016¹⁶ and PES-2018.¹⁷ Both surfaces are full-dimensional and semi-global as they describe the cis and trans conformational space of formic acid. In the following, we briefly summarize technical aspects, and for details and references, see Refs. 16 and 17.

PES-2016 was presented by Tew and Mizukami in 2016. They fitted 17 076 single point energies at the CCSD(T)-F12c/VTZ-F12 level using a method similar to LASSO constrained optimization. Note that points along the O–H torsion as defined by the internal coordinate path were included in the fit. The potential is composed of a zero-order surface using Morse oscillators for the atom–atom distances to ensure proper asymptotic behavior. Perturbations are introduced by a more flexible correction surface that is a sum of distributed multivariate Gaussian functions of combinations of the atom–atom distances. The authors report an accuracy of the fit to 0.25% with a root-mean-square deviation of 9 cm⁻¹ to the ab initio data in the energy range 0 cm⁻¹–15 000 cm⁻¹.

In 2018, Richter and Carbonnière presented a second full-dimensional formic acid potential, which they fitted to 660 single point energies at the CCSD(T)-F12a/aVTZ level using the AGAPES program. The potential is defined in a sum-of-product form using internal valence coordinates that minimize the potential energy coupling, e.g., the torsional motion is mostly localized along one dihedral coordinate. The AGAPES algorithm iteratively converges the potential energy expansion to a predefined cut-off value, in the process adding more ab initio points. The dead branching variant of AGAPES was used and the authors report the fit to be accurate between 0 cm⁻¹ and 13 327 cm⁻¹.

In this work, we solve the vibrational eigenstates and their corresponding eigenvalues using a numerical implementation of high order canonical Van Vleck perturbation theory (CVPT).^47–49 Our implementation requires a Hamiltonian expanded in terms of the normal coordinates and their conjugate momenta written in terms of raising and lowering operators.

The two potential surfaces we consider, as described above, are expanded in Taylor series in terms of the stretch, bend, and dihedral angles extensions following the coordinate choice of Richter and Carbonnière.¹⁷ One difference is that we use Simons–Parr–Finlan coordinates for the stretches to improve the convergence properties.⁵⁰ The potential includes up to four-body terms. These contributions include expansion terms up to 8th, 6th, 6th, and 4th for the one- to four-body terms, respectively. A comparison of the torsional potentials and their fits is given in Fig. S13 of the supplementary material. The exact kinetic energy operator⁵¹ is a function of the G-matrix elements. These elements are expanded through sixth order in the internal coordinates with the exception of four-body terms that are expanded through fourth order. The potential-like contribution to the kinetic energy⁵¹ is expanded through fourth order.

The potential and kinetic contributions to the Hamiltonian are re-expressed in normal coordinates and conjugate momenta. The normal coordinates are curvilinear, since they are written as a linear expansion of the internal coordinates R = LQ, where L is obtained by diagonalization of the F and G matrices following the work of Wilson, Decius, and Cross.⁵² As a final step, in order to carry out CVPT, the Hamiltonian is expressed in terms of the harmonic oscillator raising and lowering operators expressed in normal form.⁵³ 

The CVPT is carried out through a series of transformations that order by order remove off-resonance terms. The details of the calculations have been described previously, so they are not repeated here.^53,54 The central point is that the more off-diagonal terms that are removed, the poorer the agreement between subsequent orders of perturbation theory and the smaller the sizes of the matrices that need to be diagonalized. The two extremes are a final effective Hamiltonian that equals the initial Hamiltonian and a Hamiltonian with no coupling terms. There are many criteria that are used in deciding the form of the final Hamiltonian.^47,55,56 In this work, we follow an approach that was successful for describing the vibrations of HFCO,⁵⁴ a molecule with many possible resonances. In that approach, a rather large list of possible resonance terms is generated. Our criteria are based solely on energy mismatch; coupling strength is not considered. If a coupling term couples states with an energy difference less than W_i/cm⁻¹, where i is the order of the coupling terms, then the coupling term is retained. Here, we set {W₃, W₄, …, W₈} = {350, 200, 60, 50, 15, 15} and therefore retain many small coupling terms that could be transformed away, knowing that these terms are readily treated in the ensuing variational calculation of the transformed Hamiltonian.

The resulting transformed Hamiltonian is expressed as a matrix using a product basis of harmonic oscillators. The size of this basis is constrained by the number of quanta in each oscillator. We use the constraint ∑_ic_in_i ≤ N_t, where c = {1, 2, 2, 2, 3, 3, 6, 6} when the modes are ordered in terms of increasing frequency and N_t = 8. These numbers are roughly based on the frequencies and have not been optimized in a systematic way. These values lead to modest-sized matrices of 264 and 213 for trans-formic acid for the two symmetry blocks, respectively. Comparing results to N_t = 12 shows that all states are converged to better than 0.1 cm⁻¹ compared to the smaller bases for the sixth order CVPT results of trans-formic acid. Our approach allows for many possible resonance interactions through the large values of W that still yield representations for which the eigenfunctions of the effective Hamiltonians can be readily obtained.

In order to test the quality of the results, it is important to compare eigenvalues calculated at different orders of perturbation theory. We do this by comparing E(4) results to E(6) results where the order of the perturbative result is in parentheses. If one finds significant discrepancies, we consider the possible terms that are leading to the slow convergence and adjust the values of W accordingly. Atomic masses used in the CVPT calculations are reported in the supplementary material.

The convergence of newly reported CVPT vibrational eigenvalues is quantified with the three quantities {E_max, n_max, N_ord}, where N_ord is the order of perturbation theory and the other two parameters identify the states included in the comparison set. We set E_max = 4000 cm⁻¹ and n_max = 4, constraining states to have energies below the cut-off energy E_max and ∑_in_i ≤ n_max. The state identification is based on the leading coefficient in the expansion.

Convergence with respect to the level of perturbation theory is shown in Fig. 2 for trans- and cis-HCOOH where the energy difference between sixth and second as well as sixth and fourth order states is compared. In comparing two orders of perturbation theory, we compare states that are most similar to each other, and we use their overlap as a measure of this similarity. As an example, the overlap between the fourth order and sixth order states is defined by $S=U4TU6$⁠, where $UnTHnUn$ is the similarity transformation that diagonalizes the nth order Hamiltonian. At low energies, the diagonal elements S_ii are close to one. Whenever eigenvalues switch order between two levels of perturbation theory or if two states are in resonance at one level and not at another, then this information is encoded in the elements of S and can be used to match the states appropriately.

Figure 2 (left panels) shows that vibrational eigenstates of trans-HCOOH on PES-2016 are slightly better converged than on PES-2018. On the whole, the deviation between fourth and sixth order shows that eigenstates are converged within a few cm⁻¹. The convergence of individual states can be checked with data provided in Table S7 of the supplementary material. The mean absolute deviation (maximum absolute deviation) for n_max = 1, 2, 3, 4 states is below 1(2), 1(3), 1(6), and 2(7) on PES-2016 and below 1(5), 1(5), 2(13), and 3(21) on PES-2018, all in units of cm⁻¹. The comparison between sixth and second order is especially insightful as second order Van Vleck perturbation theory is the same as standard VPT2⁵⁷ with resonance treatment (VPT2+K).⁵⁸ On PES-2016, the performance for fundamentals and binary states is quite satisfactory with mean absolute deviations (maximum absolute deviations) below 1(3) cm⁻¹ and 2(6) cm⁻¹, respectively. The limitations of VPT2+K are revealed beyond two quanta where the maximum absolute deviations to sixth order CVPT amount to 19 cm⁻¹ and 16 cm⁻¹ for three- and four-quantum states, respectively (cf. Table S7).

The CVPT calculated fundamentals, as well as states with multiple quanta of excitation, converge more slowly for cis-formic acid on PES-2016 than PES-2018, which is reflected in Fig. 2 (right panels). Fundamental wavenumbers of cis-HCOOH are converged within 1 cm⁻¹ with respect to the order of perturbation on PES-2018, whereas the two high-frequency X–H stretching fundamentals on PES-2016 differ by 4 cm⁻¹ between fourth and sixth order. As expected, the O–H torsion of formic acid is very sensitive to the relative conformation, as reflected in the large down-shift from 640.73 cm⁻¹ to 493.42 cm⁻¹ for the fundamental. The frequency lowering reduces the involvement of the torsion in strong resonance polyads—groups of nearly degenerate states that are described as zero-order coupled states⁵⁹—which are characteristic of the trans species (vide infra).¹⁴ This decoupling for cis appears to improve the convergence of the remaining higher-frequency low-lying states so that overall CVPT converges faster for cis than trans, as can be observed by comparing the left- and right-hand side panels of Fig. 2.

The nine vibrations of trans-formic acid are schematically shown in Fig. 3. These localized representations are idealized, especially the C–H/D and O–H/D in-plane bending, and C–O stretching vibrations ν₄, ν₅, and ν₆ are mixed group vibrations. In this work, we apply the Herzberg nomenclature of HCOOH also to its three deuterated isotopologs. With a large number of perturbation-free experimental reference data available for fundamentals of both conformers, we can compare PES-2016 and PES-2018 using different anharmonic vibrational methods such as the internal coordinate path Hamiltonian (ICPH),¹⁶ multi-configuration time-dependent Hartree (MCTDH),^17,18 and high order canonical Van Vleck perturbation theory (CVPT).

The anharmonic fundamentals of HCOOH are listed together with experimental reference data in Table I; deviations between experiment and anharmonic calculations employing both PES are visualized in Fig. 4. These figure results include all deuterated isotopologs. Table I shows that the CVPT6 results for trans-HCOOH compare favorably with ICPH and MCTDH reference data on the same surface as they agree within 2 cm⁻¹. Comparison for cis-HCOOH indicates that the internal coordinate path Hamiltonian (ICPH) eigenvalues are not fully converged. In addition, the CVPT6/PES-2016 calculation overestimates the O–H stretch ν₁(cis) transition energy by 16 cm⁻¹. Given that CVPT4 and CVPT6 agree within 4 cm⁻¹, this discrepancy is probably due to the surface. Note that ν₁(cis) was not reported in Ref. 16 and is therefore missing in Fig. 4. MCTDH and CVPT6 results for cis compare favorably with deviations below 3 cm⁻¹ except for ν₁, which differs more (3631 cm⁻¹ MCTDH, 3636 cm⁻¹ CVPT6). Upon C–H deuteration,¹⁸ ν₁ agreement between MCTDH (3625 cm⁻¹) and CVPT6 (3637 cm⁻¹) worsens where the MCTDH predictions seem somewhat more sensitive to isotopic substitution (−6 cm⁻¹, +1 cm⁻¹ CVPT6) than the experimental data indicate (−2 cm⁻¹). On the whole, Fig. 4 demonstrates that vibrational energies for both conformers tend to be higher than experiment on PES-2016, whereas PES-2018 tends to underestimate the experiment.

The careful reader will have noticed that in Fig. 4, data for 35 instead of 4 × 9 = 36 trans fundamentals are shown. To this day, the C–H out-of-plane vibration ν₈ of HCOOD remains the only fundamental not reported in the gas phase. The anharmonic calculations on PES-2016 and PES-2018 consistently predict ν₈ to be slightly shifted down by 1 cm⁻¹–2 cm⁻¹ upon O-D isotopic substitution. Experimental data in a weakly perturbing neon matrix indicate a similar shift in the opposite direction, however.^3,10 Both overtones could be observed in this work at 2060 cm⁻¹ (HCOOH) and 2055 cm⁻¹ (HCOOD). The overall weak experimental band integrals and “pure” CI coefficients (>0.99) indicate 2ν₈ not to be perturbed (significantly) by resonances, as is the case for DCOOH and DCOOD (cf. Table IV and Fig. S9). In the context of second order perturbation theory,⁵⁷ 

one can extract an x_i,i value if one knows the transition energy of both the fundamental and first overtone of the ith mode. Using the fundamental (1033.47 cm⁻¹ from Ref. 7) and first overtone (2060 cm⁻¹) of trans-HCOOH, one obtains a diagonal anharmonicity constant of 2x_8,8 = −7 cm⁻¹. Assuming x_8,8 to remain unchanged upon O-D isotopic substitution, the gas phase band center of trans-HCOOD can be estimated to be shifted down by 2 cm⁻¹ relative to trans-HCOOH, in line with anharmonic predictions (Table IV).

While trans-HCOOH experimental vibrational data are available up to 14 000 cm⁻¹ (Ref. 7) and select OH overtone data beyond,⁶⁴ we restrict our discussion in this work to vibrational levels up to 4000 cm⁻¹. On the basis of the available calculations and symmetry information from depolarized Raman spectra, we assign peaks in the Raman jet spectra reported in Refs. 40 and 44 and this work and review previously assigned IR bands. As some of the reassignments make use of depolarization ratios, we give a brief explanation of depolarization ratios and their relation to the symmetry of a vibration as note in the references.⁶⁵ For further reading, see Refs. 66 and 67.

In order to motivate these assignments, we have examined a wide range of data. We have compiled available benchmark quality, i.e., perturbation-free, experimental reference data from the literature and show them together with Raman jet band centers in Table II for HCOOH; data for the deuterated isotopologs are shown in  Appendix B (Table IV). Newly assigned bands are marked with an asterisk and assignments that deviate from Refs. 16 and 17 (or the experimental reference if not discussed in either) are indicated by a dagger. (Re)Assignments are based on the combined predictions from ICPH,¹⁶ MCTDH,^17,18 and CVPT6 calculations and taking into account the expected error progression for each state. We use the deviation between the experiment and calculation for a fundamental as a measure of the error associated with this vibrational degree of freedom per quantum of excitation. In the following, we discuss (re)assignments concerning trans-HCOOH and touch upon deuterated isotopologs whenever necessary for the interpretation of HCOOH bands. Convergence data and CI mixing coefficients are shown for states relevant to the discussion in  Appendix A (Table III). The Raman jet spectra in the full spectral range between 870 cm⁻¹ and 3750 cm⁻¹ are shown in Figs. S1–S12 of the supplementary material. The discussion of (re)assignments is sorted by energy, starting with the low-frequency end of the spectrum. The O–H stretching spectrum between 3510 cm⁻¹ and 3630 cm⁻¹ is discussed separately at the end.

The first reassignment in Table II, as indicated by the daggers (†), concerns the bands at 1727 cm⁻¹ and 1738 cm⁻¹. These bands are assigned as the combination bands ν₆ + ν₇ and ν₆ + ν₉, respectively. One might expect the latter to be higher in energy since ν₇ lies below ν₉. However, Perrin et al. analyzed several interacting states in the ν₃ region by high-resolution IR spectroscopy and fitted ro-vibrational lines with an effective Hamiltonian. They reported (among others) band centers at 1726.40 cm⁻¹ and 1737.96 cm⁻¹ for ν₆ + ν₉ and ν₆ + ν₇, respectively.⁹ This “inverse” order was attributed to resonance interactions with ν₃ shifting ν₆ + ν₇ above ν₆ + ν₉; the latter can only interact via Coriolis couplings, while the former can additionally interact via Fermi couplings. In contrast, all four calculations in Table II indicate that the naive assignment based on the order of fundamentals is correct as was recognized by Richter and Carbonnière. They had proposed a reassignment based on improved root-mean-square deviations among both ICPH and MCTDH calculations.¹⁷ Moreover, we find that the above postulated Fermi coupling is small, being on the order of 1 cm⁻¹, this providing further evidence that the lower band should be assigned to ν₆ + ν₇.

Coincidentally, in the Raman jet spectrum of HCOOH, peaks at 1737 cm⁻¹ and (much weaker) at 1726 cm⁻¹ are observed, and both correspond to a′ transitions. The agreement with the band assigned to ν₆ + ν₇ by Perrin et al., however, is coincidental, as the Raman signal at 1737 cm⁻¹ can be ascribed to ν₃ of H¹³COOH due to the naturally occurring ¹³C isotope. The isotope shift was measured in the gas phase²³ and weakly perturbing neon matrix,³ amounting to −40 cm⁻¹, in perfect agreement with the observed −39 cm⁻¹. The corresponding ¹³C peak can also be identified in the HCOOD spectrum at 1733 cm⁻¹, shifted by −39 cm⁻¹. Using CVPT6 predictions, we can furthermore identify one band of the ν₃/2ν₈ doublet for C-deuterated and even ν₆ for O-deuterated isotopologs of formic acid-¹³C. Band centers are tabulated in Table S3 of the supplementary material

For further confirmation that the ro-vibrational lines observed by Perrin et al., that correspond to a band center at 1726 cm⁻¹, should be assigned to ν₆ + ν₇, we turn to hot bands involving the two energy levels in question. Although we cannot observe ν₆ + ν₉ directly in the Raman spectrum of HCOOH, we can calculate hot bands of ν₆ originating from the two (by far) lowest energy levels ν₇ and ν₉ using experimental data and compare those results with the observed hot bands in Fig. 5. With the assignment by Perrin et al., we expect one hot band to be blue-shifted (+6.94 cm⁻¹) and one red-shifted (−19.18 cm⁻¹) relative to the ν₆ fundamental. In the Raman spectrum, three distinctly separated hot bands can be found, shifted from the fundamental by −3 cm⁻¹, −7 cm⁻¹, and −11 cm⁻¹—all three red-shifted from ν₆—where the latter was assigned to ν₆ of cis-HCOOH by Meyer and Suhm.^11,40 An assignment of ν₆ + ν₇ and ν₆ + ν₉ inverse to Perrin et al., however, would yield shifts of −4.62 cm⁻¹ and −7.62 cm⁻¹, which are compatible with the observed hot band structure as was recognized by Meyer⁴⁴ and which is consistent with all anharmonic predictions. We therefore confirm the reassignment proposed by Richter and Carbonnière.

Contrary to HCOOH, the hot band structure around the C–O stretching fundamental ν₆ is much more complicated for DCOOH, featuring much less intensity but an overall higher density of peaks and a blue-shifted hot band. In the absence of anharmonic Raman intensities, we simulate the ν₆ spectrum using a simple model assigning a linear transition term to the C–O stretch. The experimental and simulated spectra using CVPT6/PES-2016 anharmonic energies are compared for all four H/D isotopologs in Fig. 5. We scale the trans fundamental with 3/20 to reproduce the relative experimental intensities and shift the wavenumber axis to match experiment. For the cis fundamental, the scaled transition oscillator strength of the trans fundamental is multiplied with the appropriate Boltzmann factor and the relative cis–trans Raman intensity factor from Ref. 40 (ESI, Table S3).

Considering the crudeness of the model, the simulated spectra predict the relative intensities and band contours between different isotopologs surprisingly well. Note that the choice of width for the Gaussian functions and the underestimated redshift for ν₆ + ν₇ − ν₇ (simulation −2 cm⁻¹, experiment −4 cm⁻¹) causes the hot band for HCOOH to overlap with the fundamental peak in Fig. 5. Hot bands originating from ν₇ and ν₉ overlap for O-deuterated formic acid and an accidental resonance between ν₆ + ν₉ and 3ν₉ causes ν₆ + ν₉ − ν₉ to be shifted above the fundamental for DCOOH. With a refined model including higher (anharmonic) terms in the polarizability tensor, it seems possible to assign ν₆(cis) for all three deuterated isotopologs in the Raman jet spectra shown in Fig. 5.

From the above mentioned ro-vibrational analysis of ν₃, Perrin et al. obtained band centers for several interacting states by simultaneously fitting the interacting vibrational states ν₃, ν₆ + ν₇, ν₆ + ν₉, 3ν₉, and ν₇ + 2ν₉ (the latter two originally labeled ν₅ + ν₉ and ν₅ + ν₇, respectively).⁹ We find that our Raman band centers for the ν₅ + ν₇/ν₇ + 2ν₉ Fermi doublet at 1931 cm⁻¹ and 1847 cm⁻¹ agree well with band centers reported by Freytes et al. at 1931.1 cm⁻¹ and 1847.8 cm⁻¹ (originally assigned to 3ν₇ and ν₅ + ν₇, respectively), which were also determined from the respective Q branches.⁷ The band center that Perrin et al. reported at 1843.48 cm⁻¹,⁹ however, significantly deviates from the high-resolution IR value reported by Freytes et al. This discrepancy may be resolved by updating the state assignments and therefore the form of the off-diagonal coupling elements in the fit Hamiltonian.

So far, the discussion has focused on the assignment of peaks that are observed in the Raman jet spectrum. Looking at Table II, it is surprising that while the two-quantum states ν₅ + ν₇ (1931 cm⁻¹) and 2ν₈ (2060 cm⁻¹) are observed, the state ν₄ + ν₇, which is predicted to lie between the two, is not. Following the method by Dübal and Quack,^59,68 we analyze nearby a′ symmetric states for the CI contribution of the ν₃ fundamental to the respective state, assuming that “dark” states gain intensity via resonance mixing with ν₃. The CVPT6/PES-2016 calculation predicts contributions of ν₃ to seven other states with squared CI coefficients of 0.87% (ν₇ + 2ν₉), 0.31% (2ν₈), 0.30% (ν₅ + ν₇), 0.16% (ν₈ + ν₉), 0.09% (ν₆ + ν₇), 0.02% (3ν₇), and 0.01% (ν₄ + ν₇). This simple model satisfactorily explains the absence of ν₄ + ν₇, 3ν₇, and the very weak intensity of ν₆ + ν₇. In analogy, it further illustrates the absence of 2ν₇, which has a total contribution below 0.1% from the much weaker (relative to ν₃) fundamentals ν₄, ν₅, and ν₆, and the general lack of a″ symmetric combination/overtones (cf. Tables II and IV). The (un)observed peak structure of the Raman jet spectrum of trans-HCOOH up to 2150 cm⁻¹ can entirely be understood within this simple model where “bright” fundamentals light up nearby dark states via (weak) resonance mixing.

Segments of the HCOOH and DCOOH Raman jet spectra between 2210 cm⁻¹ and 2450 cm⁻¹ are shown together with CVPT6 predictions on both surfaces in Fig. 6 to aid the following discussion. In this spectral window, four HCOOH bands are reported by Freytes et al. at 2400.2 cm⁻¹, 2376 cm⁻¹, 2338 cm⁻¹, and 2298.6 cm⁻¹. We propose reassignments in all four cases.

The initial assignment of the band at 2400.2 cm⁻¹ to 2ν₅ by Freytes et al. was rejected by Tew and Mizukami as well as Richter and Carbonnière on the grounds that their calculations predict 2ν₅ at much higher energies.^16,17 Richter and Carbonnière pointed out that taking into account the expected error progression of the two possible candidates ν₃ + ν₇ and ν₅ + ν₆, the latter is a better match since the former is expected at slightly lower energies around 2395 cm⁻¹–2396 cm⁻¹. We can confirm this reassignment as we observe two bands in our Raman spectra at 2400 cm⁻¹ and 2395 cm⁻¹, which we assign to ν₅ + ν₆ and ν₃ + ν₇, respectively.

By re-visiting early IR studies on formic acid with high quality anharmonic calculations, we are also able to assign ν₅ + ν₆ for trans-HCOOD. In one of the first spectroscopic investigations of deuterated formic acid monomer, Williams reported two bands of trans-HCOOD at 2178.8 cm⁻¹ and 2142.4 cm⁻¹, which he could not assign.⁶⁹ Nearly 70 years later, these bands remain unassigned, as far as we know. Using the anharmonic predictions on both PESs (Tables III and IV), the lower energy band can either be assigned to ν₅ + ν₆ or 2ν₇ + ν₈ and the higher energy band to ν₆ + 2ν₉ or 3ν₇ + ν₉, where ν₅ + ν₆ and ν₆ + 2ν₉ are predicted to be strongly coupled on the basis of their respective CI coefficients. The best energy matches for the bands at 2142.4 cm⁻¹ and 2178.8 cm⁻¹, and further confirmed by correction for expected error progression, are ν₅ + ν₆ and ν₆ + 2ν₉, respectively. ν₅ + ν₆ of trans-HCOOD was reported by Marushkevich et al. at 2139.8 cm⁻¹ and 2176.8 cm⁻¹ in neon and 2144.1 cm⁻¹ and 2181.0 cm⁻¹ in argon matrices—the similarity to the gas phase band centers introduces the possibility of interpreting the matrix data as the Fermi pair ν₅ + ν₆ and ν₆ + 2ν₉, as opposed to a matrix site effect.⁷⁰ 

The band at 2376 cm⁻¹ was initially assigned by Freytes et al.⁷ to ν₃ + ν₇ and tentatively reassigned by Tew and Mizukami to ν₆ + ν₇ + ν₉.¹⁶ Since ν₃ + ν₇ is assigned to 2395 cm⁻¹, this leaves ν₆ + ν₇ + ν₉ as the only plausible assignment considering that the energetically next higher state is predicted to be greater than 2410 cm⁻¹. Even taking into account expected error progressions, the predictions are still 9 cm⁻¹–14 cm⁻¹ below experiment. Based on the remarkable agreement with experiment for states up to three and four quanta, and as we shall see below, even for states involving resonances, this assignment seems unlikely, and we believe that this band does not correspond to a ground state trans-HCOOH transition.

A critical reassignment concerns the two IR bands at 2338 cm⁻¹ and 2298.6 cm⁻¹.⁷ Freytes et al. assigned the latter tentatively to ν₅ + ν₆, which we have already reassigned to 2400.2 cm⁻¹. Tew and Mizukami reassigned both bands to ν₅ + ν₈ (a″) and ν₇ + ν₈ + ν₉ (a′),¹⁶ respectively; Richter and Carbonnière noted that in light of their MCTDH results, alternative reassignments to 4ν₉ (a′) and ν₆ + 2ν₉ (a′) seem plausible.¹⁷ In the Raman spectrum of trans-HCOOH (Fig. 6, top), we observe three peaks at 2336 cm⁻¹, 2301 cm⁻¹, and 2298 cm⁻¹, all of which correspond to totally symmetric a′ transitions. This yields three or four bands in this spectral region, depending on whether the band Freytes et al. observed at 2338 cm⁻¹ corresponds with the signal we observe at 2336 cm⁻¹—the deviation of 2 cm⁻¹ is just within our experimental resolution of ±2 cm⁻¹. The ICPH calculation cannot explain the a′ band at 2336 cm⁻¹, since the “best match” ν₅ + ν₈ can be ruled out by symmetry considerations. The energetically next possible candidates are 4ν₉ at 2358 cm⁻¹ and ν₆ + 2ν₉ at 2312 cm⁻¹, which (even after expected error correction) are energetically too far away to plausibly be assigned to the Raman peak at 2336 cm⁻¹. Looking at Table II, it can be seen that the MCTDH predictions using PES-2018 are compatible with the CVPT results on PES-2016; CVPT and MCTDH predict the two aforementioned states 4ν₉ and ν₆ + 2ν₉ to be strongly coupled where the higher energy band is predicted between 2331 cm⁻¹ and 2338 cm⁻¹ and the lower energy band between 2294 cm⁻¹ and 2305 cm⁻¹.

We can gain further insights by turning to the Raman spectrum of DCOOH, as the CVPT calculations indicate that this resonance persists upon C–H deuteration. In the same spectral region, we observe five a′ bands at 2428 cm⁻¹, 2342 cm⁻¹, 2290 cm⁻¹, 2277 cm⁻¹, and 2271 cm⁻¹. The spectrum is shown in Fig. 6 alongside CVPT6 predictions in this spectral region. The former peak can be identified as ν₅ + ν₆, which is shifted toward higher energies upon C–H deuteration, as expected from the sum of fundamentals. The second last and most intense band is identified as 2ν₆. According to the CVPT results, this state has the largest CI overlap with ν₂ (cf. Table III) and therefore is expected to gain the most intensity via Fermi resonance mixing. Moreover, from the data in Table IV, we observe that if one uses the difference between the CVPT6 and experimental fundamental values for ν₆ to adjust the CVPT6 energy of the overtone, PES-2016 transitions yield 2282.2–2 × (2.4) = 2277.4 cm⁻¹. The similarly calculated value for the PES-2018 is 2272.3 + 2 × (3.1) = 2278.5 cm⁻¹, both being in excellent agreement with experiment. The two signals in the vicinity of 2ν₆ can then straightforwardly be assigned to ν₄ + ν₅ (2271 cm⁻¹) and 4ν₉ (2290 cm⁻¹) as only a′ states within ±50 cm⁻¹ (Fig. 6). The calculations allow two possible assignments for the remaining signal at 2342 cm⁻¹; ν₃ + ν₇ or ν₆ + 2ν₉. Due to the strong resonance between ν₃ and 2ν₈ in C-deuterated formic acid, combination bands in the Raman spectrum associated with this polyad are always observed as doublets (cf. Table IV). The absence of a second signal therefore rules out ν₃ + ν₇ as possible assignment. In this context, we note that the Raman band center for 2ν₆ (2277 cm⁻¹) significantly deviates from the high-resolution IR value reported by Tan et al. (2254.24 cm⁻¹).⁷¹ A likely explanation is the missing Fermi resonance coupling between ν₂ and 2ν₆, as Tan et al. included only Coriolis-type coupling in the effective fit Hamiltonian (cf. Ref. 71).

With ν₆ + 2ν₉ assigned to 2342 cm⁻¹ and 4ν₉ to 2290 cm⁻¹ for DCOOH, we are now able to assign the remaining HCOOH Raman signals at 2336 cm⁻¹, 2301 cm⁻¹, and 2298 cm⁻¹. We assign the higher energy band to ν₆ + 2ν₉ and the two remaining bands to 4ν₉ and ν₇ + ν₈ + ν₉—the order depending on whether the band observed by Freytes et al. at 2338 cm⁻¹ corresponds to the Raman band at 2336 cm⁻¹ or is indeed ν₅ + ν₈. Since ν₆ has a very strong IR oscillator strength,⁷² it seems plausible that combination bands involving ν₆ are not completely “dark” states, as 2ν₆, ν₆ + ν₈, and ν₅ + ν₆ are all observed by Freytes et al. (cf. experimental intensities reported in Table II of Ref. 7). It is therefore more convincing to assign the two IR bands at 2338 cm⁻¹ and 2298.6 cm⁻¹ to ν₆ + 2ν₉ and 4ν₉, respectively. The resulting deviation of −2 cm⁻¹ (2338 cm⁻¹ IR to 2336 cm⁻¹ Raman) in turn corroborates the assignment of the Raman signal at 2298 cm⁻¹ to 4ν₉.

A second although experimentally less certain indication for discrepancies between ICPH and CVPT/MCTDH predictions for higher excited torsional trans states comes from a blue-shifted hot band of the ν₃ fundamental of trans-HCOOH. We observe four distinct peaks at 1783 cm⁻¹, 1770 cm⁻¹, 1763 cm⁻¹, and 1757 cm⁻¹, which are shifted by +7 cm⁻¹, −6 cm⁻¹, −13 cm⁻¹, and −19 cm⁻¹ relative to the ν₃ fundamental at 1776 cm⁻¹ (see Fig. S9 of the supplementary material). Increased intensity between 1770 cm⁻¹ and 1776 cm⁻¹ indicates other hot contributions where no Q branch can be identified due to the spectral congestion. Focusing first to the down-shifted hot bands, Meyer and Suhm assigned the hot band at 1770 cm⁻¹ to ν₃ + ν₇ − ν₇, which is consistent with our assignment of ν₃ + ν₇ to 2395 cm⁻¹.¹¹ We assign the weakest band shifted by −19 cm⁻¹ to 2ν₃ − ν₃, in agreement with 2ν₃ band centers from FTIR jet⁶³ and helium nanodroplet^73,74 measurements and anharmonic calculations that predict a large redshift of 17 cm⁻¹–21 cm⁻¹. Meyer and Suhm were not able to assign the blue-shifted hot band since all DFT VPT2 anharmonicity constants x_3,j were negative, predicting exclusively redshifts. Upon O-D isotopic substitution (Fig. S9), we no longer observe a blue-shifted hot band of ν₃. This absence could be an indication of a Fermi resonance between a dark state (involving O–H motion) and a combination state involving ν₃, based on our previous discussion of a blue-shifted hot band of ν₆ in the Raman spectrum of DCOOH. The ICPH and CVPT calculations on PES-2016 indeed predict strong mixing between ν₃ + ν₇ and ν₅ + ν₆. Since an experimental band center is available for ν₅ + ν₆, we can directly compute the hot transition ν₅ + ν₆ − ν₇, which is shifted by −3 cm⁻¹ and can therefore not explain the observed blue-shifted band. The MCTDH calculation predicts strong mixing between ν₇ + 3ν₉ and ν₃ + ν₉, as does CVPT6 on both surfaces. CVPT6/PES-2016 predicts ν₄ + ν₈ to be involved in this resonance most likely due to the higher energy associated with ν₃. The CI coefficients do not allow us to unambiguously assign ν₃ + ν₉ and ν₇ + 3ν₉, and we therefore refer to them as “higher” and “lower” energy state. We obtain predicted shifts of +(3–5) cm⁻¹ and −(4–7) cm⁻¹ for hot transitions originating from ν₉ into the higher and lower energy state, respectively. The blue-shifted prediction is therefore compatible with the experimental data where the lower energy hot band could be assigned to the increased intensity between 1776 cm⁻¹ and 1770 cm⁻¹. This mixing is not predicted by the ICPH calculation, which predicts shifts of +27 cm⁻¹ and 0 cm⁻¹.

The previously discussed assignments concerning 4ν₉ and possibly ν₇ + 3ν₉ illustrate the difficulty of high level variational calculations that explicitly take into account wavefunction delocalization effects associated with the necessity to converge more and higher energy eigenvalues compared to single-reference methods.^16,75 In the unified treatment of cis- and trans-formic acid in the ICPH model, cis states correspond to higher (torsional) excitations of the global minimum trans isomer. In that treatment, 4ν₉ of trans corresponds to nν₉ with n = 6, as n = 4, 5 correspond to the ground and torsionally first excited state of cis-HCOOH.¹⁶ In light of the not fully converged cis fundamentals discussed in Sec. III B, it appears plausible that the ICPH eigenstates for 4ν₉ and ν₇ + 3ν₉ are not fully converged either—mind that the eigenvalues are above, not below the experimental value—and extension of the vibrational basis set is expected to converge these higher excited torsional ICPH eigenvalues toward experiment.

The assignment of the band at 2504 cm⁻¹ is particularly interesting as it was previously reassigned to the nearly degenerate Fermi doublet 2ν₅/ν₅ + 2ν₉, which is predicted around 2500 cm⁻¹ and 2600 cm⁻¹. In this spectral window, Freytes et al. reported two bands at 2504 cm⁻¹ and 2600 cm⁻¹, respectively.⁷ In the Raman jet spectra, we observe a very weak signal at 2504 cm⁻¹, which nicely matches the value reported by Freytes et al., nothing around 2600 cm⁻¹, and then again a band at 2678 cm⁻¹ (Figs. S5 and S6 in the supplementary material). With the aid of the anharmonic calculations, the HCOOH peak at 2678 cm⁻¹ can unambiguously be assigned to ν₄ + ν₅ (Table III). Its resonance partner, ν₄ + 2ν₉, is also predicted around 2600 cm⁻¹, complicating the assignment of the IR band at 2600 cm⁻¹.

In order to get a better understanding of these couplings, we examine the CVPT6 results in more detail. Inspection of the CVPT6 Hamiltonian in this energetic region reveals an interesting situation of different ν₅/2ν₉ polyads intersecting that we now explore. The presence of significant mixing between states sharing the polyad quantum number N_p = n₅ + n₉/2, which has previously been analyzed for N_p = 1,^12,16,17 makes it challenging to identify additional resonant interactions between zero-order states. One way to circumvent this issue is to examine the couplings in the polyad eigenstate representation. To obtain this representation, we subsume into the zero-order CVPT6 effective Hamiltonian all the couplings between states sharing common {n₁, …, n₄, n₆, …, n₈, N_p} quantum numbers. The Hamiltonian matrix of select states in this representation, which is shown in the supplementary material, is visualized in Fig. 7. The states are labeled with {n₄, n₆, N_p} quantum numbers; the blue points correspond to N_p = 2, the red and green points correspond to N_p = 1, and the remaining states have N_p = 0. The y-axis represents the energy difference between the polyad eigenvalues and either those of the 10 × 10 matrix or the fully coupled CVPT6 eigenvalues. The smaller this distance, the better the polyad representation. States within each set of common quantum numbers, i.e., {1, 0, 1}, {0, 1, 1}, and {0, 0, 2}, are uncoupled due to the choice of the zero-order Hamiltonian.

There are three important results of this representation. First, the diagonal elements of the matrix match reasonably well with the eigenvalues of the 10 × 10 matrix, with only two exceptions. The exceptions, indicated by the double headed arrow, result from an isolated resonance between the lowest energy states of the blue and red polyads. The cubic terms responsible for this interaction are the same as those that couple 2v₉ to v₆ in the original normal mode representation that was discussed earlier. We conclude that the current representation is an excellent basis for describing the coupling. Second, although not shown, we note that the interpolyad couplings are substantially different from those of the original harmonic oscillator representation. Finally, the similarities (±1 cm⁻¹) between CVPT6 eigenvalues of the full and effective 10 × 10 Hamiltonian indicate that there is only minor coupling between the states shown here and all other states. There is one exception: the C–H bending overtone 2ν₄ ({2, 0, 0}) is missing contributions from the C–H stretch (ν₂). Upon inclusion of ν₂ into the effective Hamiltonian, the {2, 0, 0} eigenvalue agrees to within 1 cm⁻¹ with the CVPT6 value.

Beside 2ν₄, there are other nearby states that seem to interact with ν₂, as reflected in the CI coefficients of CVPT6 states. We now turn our attention to those peaks. In the vicinity of ±300 cm⁻¹ around the ν₂ band center, we observe overall five peaks corresponding to cold formic acid monomer (Fig. 8, top). Based on symmetry information (all a′) and the calculated anharmonic energies, the assignments are straightforward (cf. Table II). It is instructive to return to the bright state picture, which offers a more intuitive explanation for the peak structure in the C–H stretching spectrum. The computed CVPT6/PES-2016 eigenvalues plotted alongside the spectrum in Fig. 8 show that all observed trans peaks correspond to eigenstates with significant contributions from ν₂. The bright state model further predicts three more Raman signals to be observable, which are anticipated at 2934 cm⁻¹ (2ν₇ + ν₈ + ν₉), 2934 cm⁻¹ (ν₆ + ν₇ + 2ν₉), and 2964 cm⁻¹ (ν₇ + 4ν₉). These states interact strongly among each other and one can easily recognize their connection to the interacting states ν₇ + ν₈ + ν₉, ν₆ + 2ν₉, and 4ν₉ of a similar polyad, which was discussed in Sec. IV B. Meyer was able to show that two peaks at 2960 cm⁻¹ and 2928 cm⁻¹ are hidden underneath the rotational contour of ν₂, which she could completely remove by subtracting two Raman spectra with different incident laser polarizations from each other (Fig. 8, bottom).⁴⁴ Since ν₂ of H¹³COOH is shifted down by 11.3 cm⁻¹ in the gas phase,²³ the latter peak can most likely be assigned to the ¹³C isotopolog, whereas 2ν₇ + ν₈ + ν₉ and ν₆ + ν₇ + 2ν₉ are probably hidden under nearby hot bands. As the only a′ symmetric state between 2940 cm⁻¹ and 2980 cm⁻¹ (Table III), we assign the Raman signal at 2960 cm⁻¹ to ν₇ + 4ν₉.

There is evidence in the Raman jet spectrum of DCOOH for a similar combination band involving 4ν₉. At 3352 cm⁻¹, 3404 cm⁻¹, and 3421 cm⁻¹, three peaks are observed, corresponding to cold monomer transitions (Fig. S2). Depolarized measurements reveal the former to be a totally symmetric transition. The overall weak intensity of the latter two does not allow any conclusions; however, it seems plausible to assume the same as we have no evidence of any observable a″ symmetric combination or overtone peak in the Raman jet spectrum of any isotopolog. The CVPT6 calculations allow two different assignments of the signal at 3352 cm⁻¹; ν₂ + ν₆ or ν₃ + ν₄ + ν₇. The former assignment seems more plausible, as ν₂ + ν₆ is predicted to significantly mix with ν₆ + 4ν₉ and 3ν₆, which are predicted at 3352 cm⁻¹, 3408 cm⁻¹, and 3424 cm⁻¹, respectively, accounting for the other two weak peaks observed at 3404 cm⁻¹ and 3421 cm⁻¹.

In light of the newly assigned bands, such as 4ν₉, ν₆ + 2ν₉, and ν₅ + 2ν₉/2ν₅ paired with mixing coefficients from CVPT calculations, a consistent picture emerges, which indicates that, while the interaction between ν₅ and 2ν₉ remains strong, additional coupling to ν₆ becomes pronounced above 2200 cm⁻¹ and beyond. The available experimental data for torsional states of trans-formic acid are graphically summarized in Fig. 9. The clumps of nearly degenerate states that are observed in this figure, and that are most evident for trans-DCOOD, are conveniently described as nearly degenerate zero-order coupled states known as polyads.

Before separately discussing the O–H stretching spectrum of formic acid, we finally review five IR assignments of trans-HCOOH concerning bands at 3057 cm⁻¹, 3106.5 cm⁻¹, 3275 cm⁻¹, 3826 cm⁻¹, and 3963.6 cm⁻¹ previously reported by Freytes et al.⁷ 

The former two bands are only tentatively assigned by Tew and Mizukami to ν₄ + ν₈ + ν₉ (a″) and ν₄ + ν₆ + ν₇ (a′), respectively.¹⁶ In light of all four anharmonic calculations, these assignments are plausible. However, in each case, there is a state of the opposite symmetry that seems equally plausible in terms of energetic matching. Neither band is observed in the Raman jet spectrum of HCOOH, which would otherwise help resolve both assignments with additional depolarization information.

The assignment by Freytes et al. regarding the band at 3275 cm⁻¹ (3ν₆) was tentatively adopted by Tew and Mizukami. From the perspective of the anharmonic calculations, this assignment is plausible but not unambiguous. There are several other nearby states with 4–6 quanta of excitation. However, in light of the high IR activity of the C–O stretch ν₆ and the apparent variety of observed IR combination/overtone bands that correspond to vibrations associated with C–O stretching motion (cf. Table II), this assignment seems more plausible than any alternative, so we retain it as tentative.

Freytes et al. reported an IR band at 3826 cm⁻¹, which they tentatively assigned to 2ν₄ + ν₆.⁷ Tew and Mizukami adopted this tentative assignment, and it was not further discussed by Richter and Carbonnière.^16,17 On the basis of the ICPH and CVPT6 calculations, two assignments seem plausible: 2ν₄ + ν₆ and ν₃ + 2ν₈. These states are even isoenergetic in one of the calculations (Table II). Taking into account the expected error progression, however, 2ν₄ + ν₆ is expected at slightly higher energies 3831 cm⁻¹–3835 cm⁻¹, whereas ν₃ + 2ν₈ is uniformly expected at 3826 cm⁻¹–3827 cm⁻¹—in excellent agreement with the observed 3826 cm⁻¹.

The last reassignment concerns the band observed at 3963.6 cm⁻¹.⁷ Freytes et al. originally assigned it to ν₂ + ν₈, which Tew and Mizukami tentatively adopted.¹⁶ As before, this assignment can be checked by comparing it to hot bands that involve this state. Meyer and Suhm assigned a hot band in the Raman jet spectrum of HCOOH at 2925 cm⁻¹ to ν₂ + ν₈ − ν₈ (cf. Fig. 8) based on intensity arguments from the expected Boltzmann population of ν₈ and off-diagonal anharmonicity constants from DFT VPT2 calculations.¹¹ Their assignment is a much better match with the available anharmonic calculations. Using the literature value for ν₈ (1033.47 cm⁻¹ from Ref. 7), we obtain a band center for ν₂ + ν₈ at 3958 cm⁻¹, which nicely matches the error-corrected predictions between 3955 cm⁻¹ and 3958 cm⁻¹. Table III shows that for the band at 3963.6 cm⁻¹, two assignments are possible; ν₃ + 2ν₆ and a member of a set of strongly interacting states, including ν₄, ν₅, ν₇, and ν₉. In analogy to the previous discussion regarding 3ν₆, we tentatively assign the band to ν₃ + 2ν₆, as both CO stretches are the two most IR active fundamentals of formic acid.⁷ 

Complex vibrational spectra of carboxylic acids in the O–H stretching range are usually associated with their cyclic dimers that form two nearly unstrained O–H⋯O hydrogen bonds^76,77 and not their monomers. However, even for monomeric HCOOH, perturbations by skeletal modes seem to complicate the interpretation of the O–H stretching fundamental (ν₁). This complexity lead Hurtmans et al. to postpone the ro-vibrational analysis of the high-resolution gas phase spectrum of ν₁, which “is extremely dense.”⁶ Shortly after, in 2002, Madeja et al. used helium nanodroplets to obtain rotationally resolved spectra of the ν₁ band. For trans-HCOOH, they observed perturbations of ν₁ and reported three vibrational band centers at 3566.35 cm⁻¹, 3568.63 cm⁻¹, and 3570.66 cm⁻¹. They attributed the presence of the two additional bands as being due to Fermi and Coriolis interactions with ν₂ + ν₇ and ν₂ + ν₉.⁷³ Freytes et al. suggested resonance interactions of ν₁ as a likely explanation for a higher density of ro-vibrational lines but overall less than the expected intensity in their room temperature IR spectra, noting similarities between the gas phase and helium nanodroplet spectrum.⁷ Two decades later, the ν₁ gas phase high-resolution spectrum of HCOOH remains unanalyzed, to the best of our knowledge.

The Raman scattering spectrum of HCOOH, in constrast to the low-resolution FTIR jet spectrum shown alongside the Raman jet spectrum in Fig. 10, is spectrally less congested and allows the discrimination between different HCOOH bands via their respective Q branches at 3570 cm⁻¹, 3567 cm⁻¹, and 3559 cm⁻¹ with an intensity ratio of 20:13:2 and a fourth much weaker band at 3609 cm⁻¹.⁴⁰ ν₁ can be assigned to 3570 cm⁻¹, in agreement with the literature value.⁷ The excerpt of the CVPT/PES-2016 eigenvalues in Table III shows that the density of states in this spectral region is considerably large, making assignments purely based on the anharmonic predictions very difficult for the other bands. Depolarized spectra [Fig. 10, traces 1(b) and 1(c)] reveal that all three satellites are totally symmetric and C −D substitution [Fig. 10, trace 2(a)] further shows that of these three bands, two persist with only modest shifts of −2 cm⁻¹ (3609 cm⁻¹ HCOOH to 3607 cm⁻¹ DCOOH) and +4 cm⁻¹ (3559 cm⁻¹ HCOOH to 3563 cm⁻¹ DCOOH). The band at 3567 cm⁻¹, which is sensitive to C-D isotopic substitution, can then straightforwardly be assigned to ν₂ + ν₇. The high intensity of two-thirds of ν₁ suggests possible intensity stealing via Fermi resonance interaction, which in any case can only be of modest strength due to the small shift of −3 cm⁻¹ relative to ν₁; in the simplistic picture of a two level interaction, the magnitude of the off-diagonal coupling parameter represents the lower limit of the observed shift. Significant intrinsic oscillator strength can be ruled out; otherwise, ν₂ + ν₇ would be observed in the HCOOD spectrum, which is not the case (cf. Figs. S1 and S2 in the supplementary material). While the CI coefficients in all four calculations do not suggest significant resonance mixing between ν₂ + ν₇ and ν₁, only slight modifications to the O–H stretching harmonic force constant to reproduce the experimental value show that this resonance is very sensitive to the energy of the O–H stretch and we obtain mixing ratios of ∼50%:50% and a shift of −2 cm⁻¹, using CVPT6 on PES-2016.

The currently available calculations do not allow us to unambiguously assign the other two satellites at 3609 cm⁻¹ and 3559 cm⁻¹. Experimentally, we can narrow down the list of possible resonance partners: Isotopic C-D substitution excludes the C–H stretch, in-plane bending, and out-of-plane bending vibrations ν₂, ν₄, and ν₈. The very strong Fermi resonance between the C=O stretch ν₃ and 2ν₈ in DCOOH also makes ν₃ an unlikely candidate. Depolarized spectra shown in Fig. 10 further exclude odd quantum numbers for the torsion nν₉, as all satellites are totally symmetric. We therefore conclude that these satellites correspond to a member of the interacting states ν₅, ν₆, 2ν₇, and 2ν₉, which seem to be strongly coupled in this energetic regime, as can readily be seen by inspection of the CI coefficients reported in Table III and the supplementary material to Ref. 16.

As expected, O–H deuteration detunes these resonance perturbations of ν₁ and we only observe one peak for trans-HCOOD and trans-DCOOD, which agree within 1 cm⁻¹ with high-resolution values reported in the literature.^37,78 For DCOOD, Goh et al. observed small perturbations of ro-vibrational lines with a high K value that were too weak to be analyzed.⁷⁸ For HCOOD, A’dawiah et al. observed slight perturbations of ν₁, which they ascribed to Coriolis interactions with 3ν₇ + ν₈.³⁷ By using an effective fit Hamiltonian, they obtained a band center at 2601.13 cm⁻¹ for the perturbing level. Our CVPT calculations disagree with this assignment, as 3ν₇ + ν₈ is predicted to be significantly higher in energy around 2700 cm⁻¹. The best energy match by far is ν₇ + ν₈ + 2ν₉, which is predicted between 2589 cm⁻¹ and 2604 cm⁻¹ (Tables III and IV).

Using experimental data from the literature and newly assigned in this work, Fig. 4 can be extended beyond fundamentals and their strong resonance partners. Using the same vibrational framework (CVPT), the performance of PES-2016 and PES-2018 is compared for 155 energy levels in Fig. 11, including the three deuterated isotopologs. Similar to the conclusions drawn in Secs. III and IV, it can be seen that PES-2016 almost exclusively overestimates and PES-2018 often underestimates the energy. The experiment is in many cases (halfway) between both, which is also exemplified in Fig. 6. These discrepancies between both PESs seem surprising at first sight, as the underlying ab initio methods are of similar quality.

To highlight the role of anharmonic effects, we return to the second order result of Eq. (1) and focus on the x_i,i/x_i,j values. One can extract these values experimentally using select high-resolution or jet-cooled combination or overtone transitions and compare them to values obtained using the same transition energies calculated with CVPT6. It is more instructive, however, to compare to x_i,i/x_i,j values obtained with second order perturbation theory. These values are familiar to workers in the field and allow for a global comparison of anharmonic effects with 45 parameters. In contrast, the CVPT6 Hamiltonian has over 4000 terms. Before we make these comparisons, we note that the differences between the two theoretical results arise through resonance effects and higher order perturbative contributions. Even at the second order, with the equivalent treatment of resonances, one will find differences between curvilinear and rectilinear normal mode representations of the x_i,i/x_i,j values.⁴⁸ 

We report the full x matrices for both PESs computed with CVPT2 in the supplementary material for cis- and trans-HCOOH (Table S5). The values of these terms depend on the resonance terms included in the CVPT transformations, but the same resonant terms are present for both potential surfaces. In comparing these values, we first focus on the C=O stretching degree of freedom ν₃, since there are notable differences in the fundamental values (cf. Table I). Comparing the PES-2016 and PES-2018 CVPT6 results, we find 1821 cm⁻¹ and 1810 cm⁻¹, respectively, for the cis isomer and 1783 cm⁻¹ and 1773 cm⁻¹, respectively, for the trans isomer. Inspection of individual anharmonicity constants x_3,i, obtained directly with CVPT2, reveals maximum absolute deviations of 2 cm⁻¹ for trans but 7 cm⁻¹ for cis-HCOOH. More generally, the anharmonic terms for trans-formic acid are similar for the two surfaces, and the resulting deviations in anharmonic energy levels can primarily be ascribed to deviations in the underlying harmonic potential. The argument holds when comparisons are made to experimentally extracted⁶³ x_1,j values for ν₁ of trans-HCOOH (see Table S6). For cis-formic acid, the differences between the surfaces in both x_i,j values and harmonic values are more pronounced. In several instances, these differences tend to cancel to yield similar fundamental transitions (Table I).

Finally, we turn to the large-amplitude O–H torsion that connects the two conformational isomers. In the energetic regime so far experimentally explored, diagonal anharmonic contributions and resonance couplings are equally well captured by both surfaces. This is elucidated in Fig. 12 where the energetic differences between subsequent energy levels in the trans well are compared. Deviations between higher and lower orders of perturbation theory (CVPT6 vs CVPT2) and anharmonic contributions from both surfaces (CVPT6 on PES-2016 vs PES-2018) appear to become substantial starting with 5ν₉. Delocalization effects between the cis and trans well are expected to start with 6ν₉.¹⁶ The precise energy of 6ν₉ should be a sensitive marker for high level calculations due to additional wavefunction mixing with the cis well and anharmonic potential contributions close to the isomerization threshold.

Utilizing the previously published CCSD(T)-F12 quality potential energy surfaces of formic acid in conjunction with high order canonical Van Vleck perturbation theory, we have revisited the IR and Raman spectra of monomeric trans-formic acid and its three deuterated isotopologs. Overall, we were able to add 11 new vibrational band centers of trans-HCOOH to the database and reassign seven previous IR assignments. We also assigned 53 new vibrational band centers of deuterated trans isotopologs that significantly extend the database, which until now was mostly composed of fundamentals.

Three strategies were pursued in making these assignments. In the absence of anharmonic Raman intensities, the simple model of dark states gaining intensity via resonance mixing with nearby fundamentals proved successful in providing a qualitative explanation for a plethora of (un)observed combination and overtone peaks in the Raman jet spectra of formic acid monomer. Using temperature as a population control, we have identified many hot bands and, in selected cases, used them to critically review other assignments. These bands allow direct insights into the anharmonic contributions to the Hamiltonian. Still many more observed hot bands remain unassigned, and a future direction is to incorporate the calculation of Raman activities into the CVPT6 framework. We have also made use of multiple isotopologs to probe the quality of the potential surface to aid in peak assignments. The approach included identifying bands in the Raman spectra that correspond to the two CO stretching vibrations (ν₃ and ν₆) of naturally occurring formic acid-¹³C. The resulting extensive vibrational reference data, particularly highly excited torsional and many resonance coupled states, ultimately will help establish formic acid as a model system to test and evaluate anharmonic vibrational methods for polyatomic molecules.

Comparison between second through sixth order Van Vleck perturbation theory validates the popular approach of second order vibrational perturbation theory with additional resonance treatment for fundamentals and binary combinations/overtones but also shows its limitations for quantitatively predicting highly excited vibrational states beyond two quanta. Errors for three- and four-quantum states were shown to be up to 20 cm⁻¹, even more with the increase in energy and excitation of vibrational quanta. The close experiment-theory interplay revealed increasing deviations of internal coordinate path Hamiltonian eigenvalues to the experiment for highly excited torsional states, further demonstrating the importance of close collaboration between experiment and theory.⁷⁹ While the question to this discrepancy remains open, several indications point to unconverged torsional eigenvalues.

As we have discussed and shown in Fig. 9, the large-amplitude O–H torsion leads to interesting and challenging dynamical features of trans-formic acid. The data shown in Fig. 9 illustrate how this torsion coupling extends to multiple isotopologs. The strong coupling among states, which share the polyad quantum number N_p = n₅ + n₉/2, is present for many of the isotopologs. Moreover, our analysis of these states (see Fig. 7), which shows a breakdown of the polyad quantum number N_p = n₅ + n₉/2 for N_p ≥ 2 due to mixing with the C–O stretch ν₆, extends to trans-DCOOH starting with N_p = 1.5. The resonance interactions involving the torsion extend to the cis well despite its lower frequency, and they are expected for 6ν₉,¹⁶ which can be investigated with appropriate methods, such as the internal coordinate path Hamiltonian. This work is intended to help calibrate future endeavors near the isomerization threshold.

Rotationally cold Raman jet spectra of HCOOH and DCOOH allowed us to identify perturbing states via their respective Q branches that interact with the O–H stretching fundamental ν₁, leading to a complicated ro-vibrational spectrum.⁶ We have presented experimental evidence for the involvement of the aforementioned polyad (possibly including the OCO bend) in this resonance with ν₁. Next steps in understanding the O–H stretching and torsion dynamics include refinement of the existing potential energy surfaces and analysis of the high-resolution spectrum of ν₁, which we strongly encourage. Better understanding of the interaction between O–H stretching and skeletal modes of isolated formic acid is an important stepping stone in elucidating the complex vibrational spectrum of the cyclic formic acid dimer and gaining new physical insight into the dynamics of strongly hydrogen bonded systems. This paves the way for theoretical work revisiting the vibrational spectrum of monomeric formic acid beyond 4000 cm⁻¹, as many of the band centers reported by Freytes et al. remain tentative.⁷ Further investigation along the lines presented in this work will deepen the understanding of coupling across chemical bonds as resonance mixing between the CH and OH stretches, and short time intramolecular vibration redistribution (IVR) dynamics become significant in the OH overtone spectra of formic acid.^13,64

See the supplementary material for a full set of computed CVPT eigenvalues for H/D isotopologs of trans-formic acid, fundamentals of cis-formic acid and trans-formic acid-¹³C, and Raman jet spectra in the spectral range between 870 cm⁻¹ and 3750 cm⁻¹.

A.N. gratefully acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG; Project No. 389479699/GRK2455), thanks the Fonds der Chemischen Industrie (FCI) for a generous scholarship, and thanks Katharina A. E. Meyer, Martin A. Suhm, and David P. Tew for valuable discussions. E.L.S. gratefully acknowledges support from the NSF via Grant No. CHE-1900095. The authors thank Falk Richter for kindly supplying ab initio harmonic frequencies of trans-HCOOH at the CCSD(T)-F12a/aVTZ level.

The data that support the findings of this study are available within the article (and its supplementary material).

To aid the discussion of trans-formic acid band assignments in Sec. IV, we report selected vibrational eigenstates together with squared coefficients of the two leading contributing basis states in Table III.