Vibrational transition moments of CH$_4$ from first principles

A new nine-dimensional (9D), \textit{ab initio} electric dipole moment surface (DMS) of methane in its ground electronic state is presented. The DMS is computed using an explicitly correlated coupled cluster CCSD(T)-F12 method in conjunction with an F12-optimized correlation consistent basis set of the TZ-family. A symmetrized molecular bond representation is used to parameterise the 9D DMS in terms of sixth-order polynomials. Vibrational transition moments as well as band intensities for a large number of IR-active vibrational bands of $^{12}$CH$_4$ are computed by vibrationally averaging the \textit{ab initio} dipole moment components. The vibrational wavefunctions required for these averages are computed variationally using the program TROVE and a new `spectroscopic' $^{12}$CH$_4$ potential energy surface. The new DMS will be used to produce a hot line list for $^{12}$CH$_4$.


Introduction
Methane plays an important role in atmospheric and astrophysical chemistry. Its rotationvibration spectrum is of key importance for models of the atmospheres of bodies ranging from Titan to cool stars. However the lack of precise data on methane spectra, particularly at higher temperatures, has severely limited models for atmospheres as diverse as Jupiter [1], exoplanets [2,3] and brown dwarfs [4]. Any temperature-dependent model of the methane spectrum requires a reliable treatment of the associated transition intensities. The construction of a reliable electric dipole moment surface (DMS) which can be used to generate such intensities concerns us in this paper.
Theoretically, a number of ab initio potential energy surfaces (PES) computed at different levels of theory have been reported [5][6][7][8][9][10]. The surface by Schwenke and Partridge [7] has been used in many benchmark studies as well as for accurate predictions of vibrational and ro-vibrational energies of CH 4 . However, to our knowledge no accurate and comprehensive ab initio DMSs of CH 4 have been published so far. In some cases ab initio DMSs were used, but not actually provided, for example by Oyanagi et al. [9] (at the CCSD(T)/cc-pVTZ and MP2/cc-pVTZ level of theory) and Warmbier et al. [11] (RCCSD(T)/aug-cc-pVTZ). An empirical DMS of limited accuracy was developed by Hollenstein et al. [12].
The goal of this paper is to bridge this gap and provide an accurate and detailed ab initio DMS of CH 4 . To this end we employed an explicitly correlated (F12) coupled cluster method in conjunction with an F12-optimized basis set to generate a nine-dimensional (9D) DMS of CH 4 . This surface is then used to compute vibrational transition moments and vibrational band intensities of 12 CH 4 for a large number of IR-active transitions in the range between 0 and 10000 cm −1 . This is done by vibrationally averaging the DMS over vibrational functions of 12 CH 4 computed through variationally solving the nuclear motion Schrödinger equation with the TROVE program [13]. This program has been used successfully for accurate and extensive ro-vibrational calculations on different tetratomic molecules [14][15][16][17][18][19][20][21]. The present work is the first application of TROVE to a larger system. Other variational approaches capable of treating this pentatomic molecule include the Lanczos-method developed by Carrington and Wang [22][23][24][25], MULTIMODE [26][27][28][29], and GENIUSH [30] as well as the methods of Halonen [31], Xie and Tennyson [32], Schwenke [33], and Nikitin [34]. The rigorous quantum dynamics algorithm employed by Yu [35,36] should be also mentioned as an alternative accurate approach employed for calculating vibrational energies of CH 4 .
On the experimental side, electric transition dipole moments for CH 4 are the subject of a large number of publications [37][38][39][40][41][42][43][44][45][46][47], where effective dipole moment parameters were derived 2 either through fits to intensity measurements or using perturbation theory. The most recent and comprehensive collection of the effective dipole moment parameters describing excitations in the region up to 4800 cm −1 was reported recently by Albert et al. [47]. Many of these studies also provide band intensities. We use some of these data to validate our vibrational transition moments as well as our DMS, where possible.
The paper is structured as follows. In Section 2 the ab initio method used for generating DMS of CH 4 is described and an analytical representation for this DMS is introduced. The PES is detailed in Section 3. The variational approach is described in Section 4 and theoretical vibrational transition moments and band intensities are reported in Section 5. Finally, Section 6 offers some discussion and future prospects.

Ab initio calculations
We used the recently proposed explicitly correlated F12 coupled cluster method with single and double excitations including a perturbational estimate of connected triple excitations, CCSD(T)-F12c [48] in conjunction with the corresponding F12-optimized correlation-consistent basis sets, namely aug-cc-pVTZ-F12 [49] (augmented correlation-consistent polarized valence triple-zeta basis) as implemented in Molpro 2010 [50].
Methane has nine vibrational degrees of freedom and thus requires a very extensive grid of geometries in the ab initio calculations. We have computed three ab initio DMSs (one for each Cartesian component) on a fine grid of 114,000 geometries covering the energy range up to 50,000 cm −1 with C-H distances ranging from 0.75 to 2.0Å, and interbond angles ranging from 50 to 140 • . The finite field method, with a field of 0.005 a.u., was used to compute the dipole moment component as first derivatives with respect to the external electric field in the dipole approximation.
In order to represent the ab initio dipole moment vector µ of CH 4 analytically we use the so-called symmetrized molecular bond (SMB) representation [8], where µ is projected onto three symmetrized combinations of vectors associated with the four molecular bonds C-H i (i = 1, 2, 3, 4). An analogous approach was used to represent the ab initio dipole moment of NH 3 [51] (see also Ref. [8] for a non-symmetrized MB representation of the methane DMS). The three components of the dipole moment of CH 4 span three components of the F 2 representation, i.e.
F 2x , F 2y , and F 2z in the irreducible representation of the T d (M) molecular symmetry group [52].
We define the following three symmetrically independent reference vectors n Γ spanning the same 3 representation F 2 : where r i is a molecular bond vector pointing from C to H i (i = 1, 2, 3, 4). The dipole moment vector µ can now be represented as where the projections µ F 2α (α = x, y, z) are the dipole moment functions also spanning the F 2 representation. Each component µ F 2α is then expanded as in terms of the following set of nine independent internal coordinates: where r i are the bond lengths and α ij are the interbond angles. Following Hollenstein et al. [12], we have introduced the factor exp −β(r i − r e ) 2 in order to keep the expansion in Eq. (5) from diverging at large r i . Here r e is a reference expansion center which is conveniently taken at a value close to the molecular equilibrium.
The properties of the components µ F 2α associated with F 2 impose symmetry relations between the expansion coefficients µ α i 1 ,...,i 9 . These relations can be reconstructed by applying the standard symmetry transformation rules [52,53] between the F 2x , F 2y , F 2z components of F 2 .
We employed the Maple program function 'solve' to find a set of independent coefficients µ α i 1 ,...,i 9 by solving an over-determined system of linear equations (see Ref. [54] for details on the computational procedure). For example, to first order the dipole moment functions µ α are given by 4 the expansions where the shorthand notations µ 7 = µ α 0,0,0,0,0,0,1,0,0 and µ 1 = µ α 1,0,0,0,0,0,0,0,0 are used. The complete expansion set is given in the supplementary material to this paper. A full sixth-order symmetrized expansion requires 680 symmetrically independent dipole moment parameters in total. Taking into account the symmetry relations, the expansion Eq. (5) can be written as: The actual values of expansion parameters µ α i 1 ,i 2 ,i 3 ,...,i 9 corresponding to our ab initio DMS were determined through least-squares fittings to the 114,000 ab initio dipole moment values.
Some parameters with large standard errors that exhibit strong correlation were constrained to zero. The ab initio data in these fits were weighted using the following expression proposed by Partridge and Schwenke [55]: where N is a normalization factor defined by i w 1 = 1 and V i is the corresponding ab initio energy at the ith geometry (in cm −1 ), measured relative to the equilibrium energy and computed using the same level of theory. The ab initio energy V i is weighted by the factor w i in the PES fitting; these weight factors favour the energies below 18 000 cm −1 . In the end 296 parameters were needed to obtain a weighted (w i ) root-mean-squares (rms) error of 0.00016 D.
The expansion parameters are given in the supplementary materials of the paper along with a fortran 95 routine incorporating the corresponding analytical form.

Potential energy surface
In this work we use a new 'spectroscopic' PES generated as follows. First, an ab initio PES was generated on the same grid of 114,000 geometries employing the same CCSD(T)-F12c approach [48] as for the DMS, but with the larger aug-cc-pVQZ-F12 basis set [49]. The frozen core approximation was applied in the coupled cluster calculations. We did not include corevalence corrections and higher-order coupled cluster corrections because they cancel to a large 5 extent (see the work by Yachmenev et al. [18]), and also because this would be very costly with the resources available to us. Relativistic corrections were computed using the Douglas-Kroll method as implemented in Molpro 2010. With this method we obtained the equilibrium bond length r e = 1.08734Å, which is close to the value r e = 1.08595Å derived in a combined experimental/ab initio analysis [56]. The ab initio potential energies have been represented using a symmetrized analytical expansion in terms of the four stretching coordinates 1 − exp(−a(r i − r e )) (i = 1, 2, 3, 4) and five bending coordinates ξ j (j = 5, 6,7,8,9) from Eqs. (6)(7)(8)(9)(10)(11). The symmetry relations between potential parameters in the T d (M) group were derived employing Maple with the procedure described above. In total 287 symmetrically independent parameters are needed for a full sixth-order expansion. These expansion parameters were obtained from a fitting procedure using the weighting scheme given in Eq. (16). Some parameters showing large standard errors were constrained to zero, which resulted in 268 potential parameters needed to get a weighted rms error of 0.3 cm −1 .
To improve the quality of the ab initio PES, we slightly refined it by fitting to the 'experimental' energies of 12 CH 4 extracted from the HITRAN 2008 [57] database for J = 0, 1, 2, 3, 4.
In these fittings the variational program TROVE [13] was used, closely following the approach described in detail by Yurchenko et al. [58]. This new 'spectroscopic' PES for 12 CH 4 will be the subject of further improvements and a future publication. Both the initial ab initio and refined potential parameters are given as supplementary material along with a fortran 95 program containing the corresponding potential energy function.

Variational calculations
The vibrational wavefunctions and energies of 12 CH 4 were calculated using the variational program TROVE [13]. In this approach the Hamiltonian operator is represented as an expansion around a reference geometry taken presently at the molecular equilibrium. The kinetic energy operator is expanded in terms of the nine linearized coordinates ξ ℓ i , which are linearized versions of the internal coordinates ξ i (i = 1, . . . , 9) defined above; for the potential energy function we use 1 − exp(−aξ ℓ i ) (i = 1, 2, 3, 4) for the four stretching and ξ ℓ j (j = 5, 6,7,8,9) for the five bending linearized coordinates. The latter expansions were applied to the refined PES introduced above. The kinetic and potential energy parts were truncated after the sixth and eighth order, respectively.
To construct the basis set, TROVE employs a multi-step contraction scheme based on the following polyad truncation. The polyad number in case of CH 4 is given by where v i is a vibrational quantum number associated with a one-dimensional primitive basis 9). At Step 1, nine sets of φ v i (ξ ℓ i ) are generated using the Numerov-Cooley [59,60] method by solving the 1D Schrödinger equations for each of the nine modes separately. The corresponding 1D Hamiltonian operators are obtained from the 9D Hamiltonian operator (J = 0) by freezing all but one vibrational coordinates at their equilibrium values. At Step 2 the 9D coordinate space is divided into the three reduced sub-spaces, . This division is dictated by symmetry: each of the sub-spaces is symmetrically independent and thus can be processed separately. For each of these sub-spaces, the Schrödinger equations are solved for the corresponding reduced Hamiltonian operators employing as basis the products of the corresponding primitive functions φ v i (ξ ℓ i ) generated at Step 1. The sub-space bases are truncated using the condition P ≤ P max , where P is given by Eq. (17). In this work we use P max = 12, i.e. our reduced sub-space basis sets are defined by The three sets of eigenfunctions Ψ (i) λ 3 , which are contracted with P ≤ 12 and symmetrized again using the standard technique for transformations to irreducible representations (see, for example, Ref. [52]). Thus our contracted vibrational basis set consists of ten symmetrically independent groups, one for each irreducible representation of the T d (M) group A 1 , A 2 , E a ,  Table 1 when used with the same approach and basis set employed in the refinements, i.e. 7 P max = 12. However the vibrational transition moments reported in the next section appear to be rather insensitive to the size of the basis set and the shape of the potential.

Transition moments
The transition moment connecting two states i ('initial') and f ('final') can be defined as the vibrational average of a dipole moment component µ α (α = x, y, z): where µ α (α = x, y, z) are the Cartesian components of the dipole moment in the molecule- . It should be noted that the transition moments defined by Eq. (21) depend on the choice of the molecular-fixed coordinate system [65] as well as on the particular choice of the transformational properties of the irreducible representations E, F 1 , and F 2 (we use irrep-matrices from Hougen's monograph [66]). The following quantity, averaged over the Cartesian components α = x, y, z and any degenerate components of the wavefunctions, is independent of these choices: and the vibrational transition momentμ if is defined as the length of the vector μ  Tables 1-3 for each transition satisfy the following set of relations: where indices i, f are omitted for simplicity. Therefore Tables 1-3 , such as the ones given in Tables 1-3.
The transition moments presented here can be used to derive the effective dipole moment parameters of 12 CH 4 , which are important for modelling intensities of this molecule (see Ref. [47] for the most recent collection of the effective dipole moment parameters). Examples of such an exercise can be found in Refs. [43,69]. For example, the parameter V 0(0,0A1) 0000A 1 −0001F 2 = 0.098677(88) D [47] is in good agreement with our value √ 3μ  Table 1). Here √ 3 is a conversion factor associated with the Dyad-GS and Pentad-GS zeroorder effective dipole constants in the tensorial representation (see, for example, Refs. [38,39]).
For higher vibrational order parameters the correspondence between our transition moments and effective dipole moment parameters is less straightforward owing to the complexity of the effective dipole moment operators which are usually given in a tensorial representation.
From the transition moments computed we also generated vibrational band intensities as given by Brown et al. [38] and [39]:  Tables 1-3, where they are compared to the experimental band intensities from Refs. [39,39,42,[45][46][47], where available. The agreement is generally very good. In some cases, particularly for combinational bands, there is some ambiguity over how 9 to distribute the individual lines between transitions with the same normal-mode vibrational quantum numbers (see also footnotes to Tables 1-3). This introduces an extra uncertainty into these comparisons.
As an illustration of a possible application of our transition moments, Fig. 1 shows the vibrational band intensities calculated at T = 298.15 K (empty blue squares) and T = 1500 K (filled red circles) using Eq. (25) (Q v = 3.038 at T = 1500 K).

Conclusion
In this work we present a new ab initio DMS for CH 4 which is used to generate the vibrational    [47], which includes both F 1 and F 2 sub-components agrees well with our S (calc) The total intensity of the 3ν 4 band, 0.83 cm −1 atm −2 Albert et al. [47], which includes A 1 and F 1 sub-components in addition to two F 2 ones can be compared to the sum of the two 'calc' band intensities S (calc) 0003F 2 . c The total intensity of the 2ν 2 + ν 4 band which includes three sub-bands (F 1 , F 2 and F 2 ) by Albert et al. [47]