R-matrix calculations of low-energy electron collisions with methane

R-matrix calculations are performed for electron collision with CH4 at energies between 0.02 and 15 eV using a series of different ab initio models for both the target and the full scattering system. A target model similar to the standard multi-reference configuration interaction used in electronic structure calculations is found to give the best results. Results are presented for elastic scattering, with particular emphasis on the Ramsauer–Townsend miminum, and for rotational excitation, momentum transfer and electron impact dissociation. Extensive comparisons are made with previous studies.


Introduction
Cross sections for electron collisions with methane are important for a number of different applications including combustion (Goodings et al 1979, Prager et al 2007 and plasma-enhanced combustion (Wisman et al 2007), the atmosphere of Titan (Cravens et al 2010) and chemical vapour deposition (Baek et al 2013). As a result a number of compilations of recommended values for these cross sections have been performed (Morgan 1992, Shirai et al 2002, Kato et al 2009, Reiter and Janev 2010, Fuss et al 2010, Song et al 2014. Individual experimental studies are discussed later in the paper. Methane has become a standard system for testing theoretical methods (Lengsfield III et al 1991, Winstead et al 1993, Bettega et al 1993, Nestmann et al 1994. However it is well-established that close-coupling based methods have difficulty converging the polarization potential (Gil et al 1994, Varambhia et al 2008. Theoretical treatments have considered elastic scattering (Boesten and Tanaka 1991, Jain and Baluja 1992, Machado et al 2002, as well as electron impact rotational (Jain and Thompson 1983, Brescansin et al 1989, Machado et al 2002, vibrational (Althorpe et al 1995(Althorpe et al , Čurík et al 2008 and electronic excitation (Gil et al 1994, Bettega et al 1998, Winstead et al 1993, Kato et al 2009. Recently Ziółlkowski et al (2012) used a closecoupling R-matrix calculation to obtain electronically inelastic collisional excitation cross sections; they then used a high-level electronic structure calculation to determine excited state energies and derivative couplings, and trajectory surface hopping to determine branching in the dissociation of the methane. None of these theoretical treatments provide a comprehensive solution to the low-energy electron scattering problem. Such a solution would, for example, provide a good representation of the well-known Ramsauer-Townsend minimum, which is very sensitive to the treatment of polarization, and at the same time consider electron impact electronic excitation which, in the case of methane, means representing diffuse electronically excited states which have a significant Rydberg character.
Rydberg states are not usually well represented in closecoupling expansions based upon standard treatments of the target molecule electronic structure. Special procedures have been proposed by Gil et al (1994) and Rozum et al (2003) for treating Rydberg states in close-coupling studies. More generally the molecular R-matrix with pseudostates (RMPS) method Tennyson 2004, 2005) has been demonstrated to give an excellent representation of polarization effects Tennyson 2008, Zhang et al 2011). In this work we develop a comprehensive, ab initio model for the low energy scattering of electrons from methane which considers elastic scattering and rotational excitation at energies spanning the Ramsauer-Townsend minimum, as well as electron impact electronic excitation and hence impact dissociation.
Our aim was to create a single model which replicated the specific features of the low-energy electron-methane collision system simultaneously. To do this we concentrated on three main aspects of the problem: the location (and presence) of the Ramsauer-Townsend minimum, the polarizability of the methane target, and the thresholds for the electronic excitations of the target. The following section gives an overview of the theory used, and then describes our attempts to develop a unified model. Section 3 presents and discusses our results with the final section giving conclusions.

Theory
Calculations were all performed using the R-matrix method as implemented in the Quantemol-N expert system (Tennyson et al 2007) and the recently updated UKRMol programs (Carr et al 2012). This methodology has been thoroughly reviewed (Tennyson 2010, Burke 2011) and we will therefore only consider aspects of the problem associated with this paper.
The basic idea of the R-matrix method is the division of space into an inner region, which encloses the entire charge cloud of the N-electron target molecule, and an outer region. Within the inner region the interaction of the scattering electron with the target is complicated: both correlation and exchange effects need to be considered in detail. The + N 1 electron wavefunction in this region is generally written: where Φ i N is the wave function of ith target state. For a many electron system such as methane, the target wave function itself is represented by a sum over configurations: (1), u ij are the extra orbitals introduced to represent the scattering electron, here represented by Gaussian type orbitals (GTOs) (Faure et al 2002). The electrons in the scattering wave function must obey the Pauli principle and are therefore anti-symmetrized by operator . The UKRMol codes use a particularly efficient procedure for treating wavefunctions in this form (Tennyson 1996a).
The second summation in equation (1) involves configurations which have no amplitude on the R-matrix boundary and where all electrons are placed in orbitals associated with the target. Since they are confined to a finite volume of space they are often referred to as L 2 configurations. Such configurations allow for relaxation of the orthogonalization between the continuum orbitals and those belonging to the target and are also used to model the effects of target polarization. Different models are discussed extensively below; L 2 configurations are generated for each these models by placing an extra electron (the scattering electron) in any of the orbitals specified subject to the constraints of occupancy and overall symmetry. Doing this requires care with the phase of the overall wavefunction (Tennyson 1997).
In standard close-coupling treatments, the first summation over i in equation (1) runs over the target electronic states included in the model. Given that even without considering the target continuum, there are an infinite number of target states to consider, this sum is always truncated. The RMPS method (Bartschat et al 1996) uses the properties of the Rmatrix method to try and create effectively complete closecoupling expansions. In terms of equation (1), this is done by extending the sum over i to run over a set of pseudo-states. These states are designed to give a complete representation of all the target electronic states, including the continuum, up to some total energy but only within the confines of the R-matrix sphere. In practice this is done by adding an extra basis set at the centre of the system; for molecules this involves an additional set of even-temperered GTOs (Gorfinkiel and Tennyson 2004).
Methane in its equilibrium geometry has T d symmetry. However the polyatomic implementation of the UKRMol codes only treats Abelian groups which in practice means D h 2 and its subgroups (Morgan et al 1997). Here methane was treated using C v 2 symmetry. D 2 symmetry can also be used and allows only a single H atom to be defined. However, tests found that D 2 calculations yield the same results, with no noticable computational advantages. Care was taken in all models to preserve the degeneracies present in a fully symmetrized calculation. Where possible our results are presented below using T d symmetry. A C-H bondlength of 1.093 95 Å was used for all calculations. All calculations used Hartree-Fock orbitals in the target representation.

Target wavefunctions
A scattering model comprises a target representation and a treatment of the ' + N 1 electron' scattering system. Clearly the target polarizability and electronic excitation thresholds are properties of the target model alone, whereas the treatment of the Ramsauer-Townsend (RT) minimum depends on both steps. Below we discuss a number of different target models we tested, in choosing between them it is also necessary to consider results obtained using them in scattering calculations.
target wavefunctions determined using configuration interaction (CI) based procedure. In practice each term is defined using configuration state functions (CSFs) which define the distribution of the electrons between orbitals and the associated spin-couplings. Use of a complete active space (CAS) CI representation of the target has certain advantages when balancing the N and + N 1 electron calculations (Tennyson 1996b). A CASCI target wavefunction includes all possible CSFs generated using a given orbital set and can be described as: where n electrons are frozen in core orbitals and the remaining − N n ( ) target electrons are distributed freely across a set of suitable selected valence orbitals. All our calculations froze at least the electrons in the a 1 1 (carbon 1 s) orbital. A small CASCI calculation for methane might be given by: which comprises the ground configuration plus the lowest unoccupied molecular orbital (LUMO). In practice this model gives a target polarizability which is far too small, no RT minimum in the associated scattering calculation, and electronic excitations which are too high. Given that the RT minimum is a feature of low-energy electron scattering caused by a cancellation between the repulsive, static exchange potential and the attractive polarization potential, it is to be expected that underestimating the polarizability will lead to a poor representation of the minimum. A more reasonable CAS distributes the eight valence electrons amongst the 8 active or valence orbitals, giving the following model: a a t t 1 2 , 3 , 1 , 2 . In a CASCI calculation increasing the number of active orbitals improves all aspects of the calculation. In practice, this approach offers rapidly diminishing returns as using an enlarged CAS gives only a modest improvement to the scattering calculation for a large increase in computational cost (Tennyson 1996b). Increasing the size of the active space has a larger effect on the excitation energies than the polarizability. Including extra orbitals in the active space does give an RT minimum eventually, where the number of orbitals required for it to present depends on the basis set. Using a larger basis set with more diffuse functions tended to give the RT minimum sooner. 6-311G* gave an RT minimum for all models examined, as did cc-pVDZ (though surprisingly cc-pVTZ did not). Including the first 4 LUMOs meant that all basis sets tested produced an RT minimum (with the exception of TZ), although they were all located too low in energy. The smaller basis sets also produced larger total cross sections above about 2 eV (larger by about 1 Å 2 at the peak). The same shifting of thresholds and the peak is seen when including more orbitals.
Excitation energies and polarizabilties as a function of basis set used are shown in table 1. The value of the vertical excitation threshold for the first excited state, which is a repulsive triplet state, is not particularly well determined experimentally but would appear to be in the region of 8.8 eV (Brongersma and Oosterhoff 1969) and 9.0 eV (Kanik et al 1993). As can be seen from the table, theory tends to overestimate these values with Ziółlkowski et al (2012) obtaining 10.01 eV. Our final scattering calculations presented below use a model for which the threshold to excitation is 10.58 eV, although some of our other models gave values lower than this.
The polarizabilities presented in table 1 are consistently too low, ranging from 7.5-11.5 a 0 3 , where the experimental value is 16.52 a 0 3 (Olney et al 1997). The polarizabilities can be improved beyond those in the table by including more states in the target region, but the gains are small and the computational cost high. Even if all the states are included, this approach still leads to polarizabilities which converge to less than the true value (Jones and Tennyson 2010). The dependence of the polarizability on basis and model is slightly complicated as approximate wavefunctions are often more polarizable than their more exact counterparts. This means that the calculated polarizabilty may decrease as the calculation is improved (Jones and Tennyson 2010 2.2.2. Rydberg model. Methaneʼs low-lying excited electron states are diffuse and have a strong Rydberg character. An attempt at modelling this was made including several additional diffuse functions using standard basis sets on the carbon without changing the CAS. These treatments lowered the excitation energies, but not significantly. We tested two different basis sets to represent Rydberg-like orbitals, which were added to the carbon basis set since the carbon is located at the centre-of-mass. The first basis was a set of quasi-Slater type orbitals recreated from Rozum et al (2003) as a sum of six GTOs each. The second set were GTOs taken from Nestmann et al (1994). These attempts gave large orbital spaces, the largest being a t t e [8 , 5 , 7 , 3 ] 1 1 2 ; these procedures gave significantly lower ground state energies but proved computationally expensive, where the target calculation alone took a whole day on a high end workstation; although these calculations were performed prior to the diagonalization routines being improved (Zhang et al 2011). The results of the two different basis sets were very similar. Increasing the number of orbitals in the extended space decreases the ground energy monotonically, with larger calculations giving lower results-the calculations are variational, so this is exactly as expected. This model improved upon all aspects of the CASCI target model. The results of this model were promising, however to obtain a reasonable polarizability, the number of states required in the inner region was very large-with some tests reaching 1000 states. The RMPS method is proposed as a solution to this.
2.2.3. RMPS. The previous calculations which attempted to represent Rydberg states predate the development of the molecular RMPS method Tennyson 2004, 2005). In fact the two procedures have some similarities since the RMPS method involves using the same CSFs as in the Rydberg model, but uses a distinct basis set for the RMPS orbitals-a set of even-tempered GTOs placed on the molecular centre-of-mass. These orbitals are designed to fill up all the space inside the R-matrix sphere, representing the Rydberg states leading up to ionization and continuum of states, found above the ionization limit, within this sphere.
We therefore undertook a series of RMPS calculations. For these calculations, the RMPS orbitals were represented using even-tempered GTOs with 14 functions for s, p, and d functions, each with α = 0.05 0 and β = 1.4. The orthogonalization thresholds used in the N, and + N 1 electron regions were 2 × 10 −4 and 1 × 10 −6 , respectively.
As with the Rydberg Model, the polarizability is dependent on the number of states included, and increases aymptotically with increasing numbers of states. Including all states up to 30 eV gave a value within 15% of experiment. The RMPS approach again improved all aspects of the target model over the Rydberg Model, ground energies, first excitation thresholds, and polarizabilities are shown in table 2.
The RMPS method gives a good target description but we encountered problems when it was used as part of a scattering calculation: the model predicted an unphysical bound anionic state. First, it was thought this might be a linear dependance problem (Little and Tennyson 2014), particularly because this was the first time the molecular RMPS method had been used with an atom placed on the centre-of-mass. The integral codes had to be adapted to cope with this, and the deletion thresholds to deal with linear dependence closely monitored, and adjusted. Increasing the deletion thresholds did not fix the model. Second, we tested the effects of the RMPS orbitals being used. After several other tests the spacefilling, even-tempered GTOs were replaced by the virtual orbitals generated by a cc-pVTZ calculation. However this also had little effect: the cross sections obtained after this change were very similar to the standard RMPS ones. Third, a closer inspection of the configurations used in the target model compared to those used in the scattering calculation suggested that the calculation was over correlated as the L 2 terms from equation (1) used in the RMPS calculation can contribute as double excitations of the target.
The unbalanced nature of the RMPS model is clearly shown by the behaviour of the low-energy eigenphase sums as the number of orbitals included in the RMPS procedure is increased. As shown in figure 1, the behaviour of the eigenphases changes abruptly with larger models. The eigenphases of the larger models are characteristic of the presence of a bound anionic state; these models also feature an R-matrix pole which lies below the energy of the target ground state. A standard CASCI calculation is balanced, by construction. An RMPS calculation may not be. How balanced a calculation must be to give physical results is not a simple problem, and the level of correlation between the N and + N 1 electron problems is a subtle problem (Rescigno 1994).

Multi-reference configuration interaction (MRCI).
To try and balance the RMPS procedure we tested target models which had multiple excitations out of the CAS. Given that the cc-pVTZ basis set alone gave similar results to the cc-pVTZ + RMPS GTOs sets for the single excitation RMPS model, and given the extra complexity of using multiple basis sets and multiple orthogonalizations for the target, just the cc-pVTZ basis set alone was used to produce all orbitals. This approach may not be general as the central carbon in methane locates GTOs at the centre-of-mass.
The  figure 2 with those of the RMPS model given in figure 1. The corresponding elastic cross sections show an even more dramatic change at low-energy, see figure 3. We therefore decided to employ a model based on a CASCI and double excitations from it. This form of the wavefunction is similar to that produced by quantum chemical MRCI calculations. We therefore dub this the MRCI model.
Of course the MRCI model can regarded as step towards doing a complete CASCI run with a CAS of − − a t e [1 5 , 1 4 , 1 ] 1 2 10 . The next step would be to include triple excitations. However that appears unneccesary for a balanced calculation; the improvement of double over single excitations is very significantly more than the changes found when including triples over doubles. Furthermore, the computational demands of including triple excitations is too great to be contemplated in a full calculation. The model Table 2. Ground energies, first excitation thresholds, and spherical polarizabilities (obtained by summing states up to 30 eV) for various RMPS and MRCI models. MRCI models are defined by the number of frozen orbitals in the core, and the number of possible excitations out of the CAS.

RMPS orbitals included
Ground including only single excitations gives incorrect behaviour at low energy unlike that with doubles, see the total cross sections shown in figure 3.

Outer region.
For the smaller calculations the outer region is computationally straightforward. However when the number of target states involved becomes large, this is no longer so and to obtain a good polarizability, a large number of states are required in the sum over states expression (Jones and Tennyson 2010). For low energy, this is not an issue, as high-energy states can simply be dropped in the outer region calculation with no detriment to lower energy cross sections (Rabadán and Tennyson 1997). However this is not possible for when considering electronic excitation, as these states represent the electron impact dissociation channel. This means that all the states needed to be included, which is not feasible using the standard UKRMol codes. This problem was solved for energies up to 15 eV using the PFARM , Sunderland et al 2010 codes, which perform the outer region work in a massively parallel environment. PFARM, originally designed to enable full CI calculations of electron collisions with open d-shell ions (and atoms) relevant to astrophysics, with many open channels, has recently been adapted with an interface to the UKRMol suite (Carr et al 2012). The current set of calculations was thoroughly tested for convergence of results with respect to variations of the long-range propagation distance at which the R-matrix is matched to asymptotic coupled radial solutions and with respect to the accuracy and stability of the series expansions used to form these solutions . As a result of using PFARM, with distributed pipelined task-groups of cores and parallelized linear algebra (Sunderland et al 2010), the calculation time for a set of energies with a given symmetry using (for example) 1152 cores was reduced from months to one or two hours, making the calculations possible.
Finally, the program polyDCS (Sanna and Gianturco 1998) was run to produce rotationally resolved differential cross sections (DCS).

Results and discussions
A number of target basis sets and R-matrix radii were tested, however the calculations were not particularly sensitive to these. The final model used a cc-pVTZ GTO basis set, an Rmatrix sphere radius of 13 a 0 and the MRCI model with the core, active and virtual spaces defined by a a [1 , 2 ] 1 , respectively. In the outer region R-matrices were propagated to 100.1 a 0 , plus or minus 10-20 a 0 for convergence tests. All subsequent results are from this model only, with 223 states included in the outer region (all states found up to 30 eV). Unless otherwise stated, the results are summed over the four symmetries given in a C v 2 calculation: 2 A 1 , 2 B 1 , 2 B 2 and 2 A 2 .

Eigenphases
Eigenphases are a useful diagnostic tool for comparing various models and theories. Sohn et al (1986) provide some empircal eigenphases extracted from their measured DCS using phaseshift analysis due to Nesbet (1979). Their analysis assumed that only partial waves up to d are important which is a dubious assumption as they considered energies as high as 5 eV. Ferch et al (1985) performed a modified effective range analysis of low-energy cross section measurements to give low-energy eigenphases. Very recently Fedus and Karwasz (2014) also performed an effective range analysis yielding a set of low-energy s, p and d eigenphases. Figure 4 compares our results with those of Ferch et al (1985), Sohn et al (1986) and Fedus and Karwasz (2014) for energies below 1 eV. Our 'A 1 ' calculations are actually for A 1 in C v 2 symmetry so contain contributions from other symmetries in T d . Considering only ⩽ ℓ 2, these eigenphases correspond to s+p+2d. It can be seen that even at 0.1 eV, partial waves with > ℓ 0 make a significant contribution. A comparison of our eigenphases with the 'A 1 ' eigenphases of Fedus and Karwasz (2014) gives good agreement at higher energies, at low energy the peak in our eigenphase sum appears slightly too low. This is consistent with the fact that  our model still does not recover the full polarizability, as discussed further below.

DCS
As seen in figure 5, the DCS do not agree well with experiment below 1 eV for a scattering angle below°30 . Above 1 eV, but below 5 eV, the agreement is good, but only above°3 0 . Above 5 eV, see figure 6, the agreement with experiment is excellent, and all the theories agree with just small differences in the forward scattering.
The behaviour for small angles at low energy may be due to the still incomplete representation of polarization effects. Previous calculations have shown that it is indeed at low angles and low energies that are particularly sensitive to the inclusion of polarization effects (Fujimoto et al 2014). However our calculation uses a truncated partial wave expansion with ⩽ ℓ 4. The long-range polarization potential will induce minor contributions from higher partial waves. This problem was studied by OʼMalley (1963) for collisions with rare gas atoms for which the polarizability is also the 3.5eV 4.2eV 5.0eV 10.0eV 7.5eV 6.0eV σ DCS ( A 2 sr -1)°π π π π π π π − 2 π − 2 π − 2 π − 2 π − 2 π − 2 Figure 6. Differential cross sections for energies between 3.5 and 10.0 eV. Red lines: present work; orange lines: Machado et al (2002); purple circles: Müller et al (1985); light blue circles: Shyn and Cravens (1990); yellow circles: Curry et al (1985); dark blue circles: Tanaka et al (1982).
leading long-range interaction. OʼMalley derives formulae for the high-ℓ contribution at low energy. Isaacs and Morrison (1996) considered this problem for non-polar linear molecules but we are not aware of a treatment for spherical tops like methane. A new version of the UKRMol codes is under development which will allow treatment of partial waves of arbitrarily high ℓ. This problem will be used as test case for these codes. In this context it is interesting to note that Machado et al (2002), who find significantly higher cross sections at low angles than us and hence better agreement with experiment, use ℓ up to 16 in their calculations.

Rotational excitation
Rotationally resolved DCS are also calculated, and agree well with experimental results of Müller et al (1985), as shown in figure 7. The elastic → 0 0 cross section is much larger than the excitation cross sections and the calculations show differences at low angles similar to those observed for the rotationally unresolved elastic DCS. The small excitation cross sections all have a relatively large experimental uncertainty and our calculations agree reasonably well with these measurements. Figure 8 shows that our calculated rotational excitation cross sections are in good agreement with the measured ones due to Abusalbi et al (1983) and Brescansin et al (1989).

Elastic cross section
As seen in figure 9, the calculated elastic cross section agrees well with the various experiments, the Ramsauer-Townsend minimum is present at the correct energy-albeit with a slightly lower minimum, which could be due to either to experimental cross sections being blurred by the energy resolution or our neglect of nuclear motion. Otherwise we obtain good agreement with the measured cross sections and in particular they agree with the recommendations of Song et al (2014), which were compiled from consideration of multiple experiment results, at all energies.    (2001), and Nestmann et al (1994), as well as the Schwinger variational iterative method calculations by Machado et al (2002). The recommended values of Song et al (2014) are also given. We note that our cross section is smooth and monotonic in the energy range 1 to 10 eV, while the previous R-matrix calculations gave undulant curves. Such undulations are not physical and almost certainly are the result of an incomplete representation of the continuum due to too aggressive orthogonalization. This is illustrated in figure 11 which shows the effect of changing the deletion threshold. The results with the higher threshold show the undulations characteristic of an incomplete basis. Figure 12 compares our momentum transfer cross section with the experimental results of Tanaka et al (1982), and the theoretical results of Machado et al (2002). Both the present work, and Machado et al (2002) show reasonable agreement     with the experimental values recommended by Song et al (2014), although in both cases the peak is shifted to higher energy.

Electron impact dissociation
All electronic excitations in methane lie above the dissociation limit and the low-lying excited states are all known to be dissociative. Therefore the dissociation cross section calculated here is simply a sum of all of the excitation cross sections, as they are all assumed to be dissociative. Figure 13 compares results for electron impact dissociation which shows that both experiment and theory conflict. As seen earlier in figure 10, the total cross section results differ only a small amount from one another; however the dissociation cross sections are much less consistent.
There are two issues with the near-threshold behaviour of the calculated electron impact dissociation cross sections. The first is that our theoretical model overestimates the vertical excitation energy for the lowest excited state. The second is that inclusion of nuclear motion effects is likely to significantly lower this threshold, see Stibbe and Tennyson (1998). Ziółlkowski et al (2012) addressed these issues by shifting their cross sections to lower energy by 3.5 eV which they justified because of their overestimate of the vertical excitation threshold and because the peak in their elastic scattering calculation was too high. As our vertical excitation energy is 0.3 eV lower than theirs, figure 13 also gives our cross sections shifted to lower energy by 3.2 eV. Our shifted results give excellent agreement with the near-threshold measurements of Nakano et al (1991) and Makochekanwa et al (2006); this must however be regarded as somewhat fortuitous.

Conclusion
It is difficult to develop a completely ab initio theoretical model for low-energy electron collisions with methane which recovers all the key properties of the system such as the Ramsauer-Townsend minimum and the electronic excitation energy of low-lying, diffuse excited states. In this work we test a number of possible collision models and find that the best results are obtained using a target representation which involves use of a complete active space for the valence electrons and excitations of up to two electrons from this valence space into an extended set of virtual orbitals. By analogy with quantum chemical electronic structure calculations, we call this a MRCI model. This model gives a good representation of the Ramsauer-Townsend minimum. It also appears to give a good treatment fixed-nuclei of the electronimpact dissociation problem, although a full study of this would require considering both a better treatment of the excited Rydberg states and of nuclear motion.
The treatment of polarization could perhaps be further improved by including more states in the calculation, or by freezing fewer electrons. The approaches are not exclusive, but both increase the computational cost.
The MRCI model represents a move beyond either the use of a complete active space target of or the standard RMPS representation, although it has some similarities with the latter approach. Whether it works equally well for low-energy electron collisions with other molecules will be a matter for further study. However it is clear that use of the MRCI model is computationally demanding and its application to electron collisions with targets with many active electrons will require further algorithmic developments to make the method generally tractable.