The efficient calculation of electron impact ionization cross sections with effective core potentials

Electron impact ionization is a fundamental process that is used to initiate chemical processes in a variety of environments such as the spark plugs of the internal combustion engine or as an initial step in plasma processing. Electron impact ionization cross sections have long been measured, and there are a number of theoretical methods available for their estimation. In particular, Kim and Rudd¹ developed the semi-empirical binary encounter Bethe (BEB) method, which has been widely and successfully used to compute electron impact ionization cross sections as a function of electron energy.^2,3 The BEB method is relatively easy to implement, see below, and is found to give reliable answers in most cases; indeed, a recent compilation of electron impact cross sections for water chose the BEB ionization cross section in preference to various experimental determinations.⁴ There are a number of other similar simplified methods of calculating electron impact ionization cross sections including methods due to Deutsch and Märk (DM),^5,6 and Jain and Khare.^7–9 The whole area has been comprehensively reviewed by Tanaka et al.¹⁰ 

Recently, two of us collaborated on constructing Quantemol electron collisions (QEC),¹¹ an expert system that performs electron–molecule collision calculations using the latest version of the UK molecular R-matrix codes known as UKRmol+.¹² QEC is based on the use of an electronic structure code Molpro¹³ to provide accurate and robust information on the target molecule. Both QEC and its predecessor Quantemol-N¹⁴ provide an automated implementation of the BEB procedure; in the case of QEC, this is based on the use of Hartree–Fock (HF) orbitals and kinetic energy integrals generated by Molpro. In this work, we investigate two possible amendments to this procedure.

The first improvement concerns molecules containing heavy atoms. The wavefunctions of these atoms are not well represented by standard self-consistent field (SCF) orbital sets as they neglect relativistic effects, which strongly influence all the orbitals in the atom. A standard method of dealing with this issue is the use of effective core potentials (ECPs), which eliminate the need to explicitly represent the inner orbitals but instead provide an effective potential that allows for an improved representation of the valence orbitals. As demonstrated below, use of ECPs in BEB calculations leads to significantly improved cross sections compared to experiment. ECPs have been used previously in BEB calculations,^15–17 particularly for molecules of interest in nuclear fusion,^18–20 biologically significant molecules where electron impact ionization often plays a role in energy transfer mechanisms as a result of radiotherapy,^21,22 and molecules of interest for industrial plasma applications.^23,24 Our focus here is to make an implementation of ECPs for BEB as automated as possible to allow efficient generation of these cross sections for molecules containing heavy atoms. To this end, we investigate implementing a default choice of ECPs, balancing the availability across the periodic table with accuracy. We note that the ECP method cannot account for the ionization of inner shell electrons, but that for the majority of applications, ionization cross sections are needed for energies well below the threshold to inner shell ionization.

Additionally, we test the empirical observation that the peak of electron impact ionization cross section for a given molecule scales almost linearly with the polarizability of that molecule, α₀.^25,26 This has led to the suggestion that the prediction of the BEB model can be further improved by using a simple, universal scaling factor based on this linear relationship.²⁷ We explore this model, referred to as αBEB, further below. Our aim is to design a predictive method which can be used for molecules that have not been well-characterized experimentally; therefore, use of this scaling procedure also introduces a need for accurate polarizabilities, which, in the absence of empirical data, must be computed.

The following Sec. II details the changes made to the BEB procedure in QEC to allow ECPs to be used. We benchmark various ECPs and conclude on a default ECP based on accuracy and coverage of the periodic table. In Sec. III, we detail our treatment of polarizabilities used in the αBEB calculations and provide additional information on the sources of molecular geometries used. Section IV provides computational details. In Sec. V, we benchmark the ECP and αBEB implementations with detailed comparisons to experiment where available.

The BEB cross section can be calculated from the following formula:

where $S=4πa02N(R/B)2$⁠, with a₀ being the Bohr radius, R being the Rydberg energy, N being the orbital occupation, and B being the binding energy; t = T/B, where T is the energy of the incoming electron; and u = U/B, where U is the orbital kinetic energy. The total cross section is found by summing the BEB cross sections for each of the occupied orbitals in an atom or molecule. A key feature of this method is that it does not require fitting parameters and can be used when the dipole oscillator strengths are not known. This makes the method extremely flexible for estimating ionization cross sections as the only parameters required are the occupation numbers, binding energy, and kinetic energies of each of the orbitals, all of which can be obtained through a HF calculation.

Our focus is to make the calculation of BEB cross sections as highly automated as possible. To do this, we assume Koopman’s theorem for ionization thresholds both for single ionization and double ionization, where we assume similarly to Nishimura et al.²⁸ that where the incoming electron energy exceeds the double ionization threshold, a secondary electron will be emitted via an Auger process, resulting in doubly charged ions or two singly charged fragments. It is possible to change these thresholds for those known experimentally through the QEC interface.

For heavy atoms, performing all electron (AE) calculations in standard basis sets for BEB cross sections leads to poor agreement with experiment.^15,24 This is because relativistic effects on the orbital binding and particularly kinetic energies are not included in simple HF orbital calculations. Valence electrons in an all electron calculation can penetrate close to the core of the atom resulting in a high kinetic energy integral. As ionization is more likely to occur when the electron is further away from the nucleus, it is expected to be traveling slower on an average. An overestimation of the kinetic energy of the ionizing electron in a molecular orbital, particularly the valence orbitals, leads to a smaller cross section. Previous attempts to correct this involve dividing the kinetic energy integrals by a factor often chosen as the principle quantum number of the dominant ionizing orbital.²⁹ 

An alternative approach is to use an effective core potential on the heavy atoms in the molecule.¹⁵ ECPs replace the nucleus and core shell electrons of an atom with an effective potential, and valence orbitals are replaced with pseudo-valence orbitals that cannot penetrate the core. ECPs are routinely used in quantum chemical calculations for molecules with heavy atoms and have been extensively reviewed by Dolg and Cao.³⁰ ECPs have also been implemented to calculate elastic³¹ and inelastic electron impact cross sections, using a complex Kohn variational method,³² the Schwinger multichannel method,³³ or a combination of the modified additivity rule (MAR) and the spherical complex optical potential.³⁴ 

As our approach is to make ECP BEB as automated as possible, we focus on the ECPs available in the Molpro basis library, using GeH₄ as an example molecule, for which a variety of ECPs are defined, and an all electron calculation is possible in a standard double zeta basis. The pseudopotentials of the Stuttgart/Cologne group are of the form ECPnXY, where n is the number of core electrons that have been replaced with a pseudopotential and XY describes the level of theory used to fit the pseudo-valence orbitals numerically to reproduce valence energy spectra and, thus, termed energy-consistent. If X is S, a single ion has been used, whereas if X is M, a neutral atom has been used. Y can be HF, WB, or DF for Hartree–Fock, quasi-relativistic, or full relativistic, respectively. We evaluated the BEB cross section for ECP10MDF,³⁵ ECP28MDF,³⁶ and ECP28MWB³⁷ on the germanium atom. The complimentary basis for describing the valence orbitals for these ECPs is minimal; however, additional basis functions can be added, which can help describe polarization effects in the molecule. An example of such a valence basis is ECP28MWB_VTZ, which contains up to f-type basis functions compared to ECP28MWB, which is defined for only the valence s, p functions. LANL2DZ³⁸ is the effective core potential from the Los Alamos group of double zeta quality based on scalar relativistic all electron calculations. The SBKJC³⁹ effective core potential was also tested, which has a more compact set of basis functions for the valence orbitals compared to LANL2DZ. Both LANL2DZ and SBKJC effective potentials are shape-consistent, that is, the pseudo-valence orbitals produced preserve the shape of an all electron valence orbital after some critical radius that defines the start of the valence region of the atom.

Figure 1 shows the BEB cross section for each of the ECPs described above on germane (GeH₄) with cc-pVDZ on the hydrogen atoms, as well as an all electron calculation using cc-pVDZ for all atoms. Unfortunately, for this molecule, there is only a single experimental measurement for the ionization cross section available at 100 eV from Perrin and Aarts.⁴⁰ Therefore, we compare to other theoretical methods for calculating the electron impact ionization cross section using a spherical-complex-optical-potential (SCOP) approach,⁴¹ modified additivity rule (MAR), and DM,⁴² as well as the values from NIST, which are BEB values that use the method of dividing the kinetic energy integrals by the principle quantum number of the dominant orbital involved in ionization.²⁹ 

The number of electrons included in the pseudopotential describing the core has the largest effect on the calculated BEB cross sections. The larger the core, the smaller the kinetic energy integrals for the pseudo-valence orbitals are, leading to a larger cross section at low energies compared to the all electron calculation. However, as there are no core electrons that can be removed, the high energy tail of the cross section is underestimated. There are no discernible differences between ECPs with a similarly sized core, that is, the cross sections obtained with ECP28MWB, ECP28MWB_VTZ, ECP28MDF, SBKJC, and LANLZDZ are all similar. Adding additional functions (ECP28MWB compared to ECP28MWB_VTZ) to help describe the valence orbitals appears to make little difference to the cross section. An important consideration for a highly automated implementation is that the same family of ECPs are defined for the vast majority of atoms in the periodic table, as such it was chosen to use the Stuttgart group ECPs with the largest possible core as the default ECP in QEC for BEB calculations. Where possible, this uses the pseudopotentials defined using the quasi-relativistic MWB, unless this is unavailable for a particular atom (often atoms in groups 1 and 2), where MDF or SDF is used.

It should be noted that a disadvantage of using an ECP is that all the inner shell orbitals are removed. These orbitals may contribute little to the BEB cross section at low energies; however, the high energy tail of an ECP calculated BEB cross section can be expected to be underestimated as a result of the missing core ionization. It is possible to perform an all electron calculation and use the contribution to BEB from core orbitals to supplement an ECP BEB calculation.¹⁶ However, this approach becomes more challenging for a black box implementation, particularly for heavy atoms due to the availability of all electron bases in the Molpro basis set library.

A number of authors have suggested that scaling the BEB cross section by some (fixed) ratio of the peak of the cross section to the target polarizability improves agreement with experiment.^25,27,43 We refer to this approach as αBEB below and assume that the relevant polarizability is the spherically averaged one, α₀, for the given system. There are broadly two issues with this approach. First, that not all authors use the same slope for their αBEB relationship, and second, the source of polarizabilities used. The CRC (chemical rubber company) provides a compilation of polarizabilities,⁴⁴ but these do not agree well with the supposedly similar compilation provided by the NIST (National Institute of Science and Technology) in the experimental section of their Computational Chemistry Comparison and Benchmark Database (CCCBDB).⁴⁵ In general, we find broad agreement between CCCBDB values and theoretical calculations performed by us and others; we, therefore, adopt CCCBDB values, where available, below. In the absence of measured polarizabilities, it is of course possible to compute them, and this would be necessary for an αBEB procedure to form part of a black box implementation; we note that such an implementation is particularly useful for studying radicals for which measured polarizabilities are generally not available.

In Table I, we tabulate computational polarizabilities calculated using the Molpro 2019 package⁴⁶ and compare them to the experimental values found in the NIST.⁴⁵ It was found that the d-aug-cc-pVDZ or, if not possible, the aug-cc-pVDZ basis set using the second order Møller–Plesset perturbation theory, MP2, method resulted in average analytical static dipole polarizabilities, which were closest to experimental values at a reasonable computational cost. It is evident that our calculations systematically underestimate the polarizabilities, giving an average of 96.3% of the observed value with a range of 11.2%. As the αBEB typically results in an increase in the BEB cross section of roughly this magnitude, these differences suggest that these theoretical polarizabilities may need to be improved if one want to use them in a scaling procedure.

The relationship between the dipole polarizabilities and the total ionization cross section maximum is shown by a plot of experimental polarizabilities^45,47–49 against the experimental total ionization cross section maximum^50–58 for the nine molecules, which are shown in Fig. 2. The gradient of a linear fit through the origin was used to predict an absolute value for the maximum of the total ionization cross section.^25,27,43 The difference between this predicted maximum and the maximum of the BEB curve can be used as a scaling factor, and thus, the BEB curves were corrected to pass through this maximum,

where M is the gradient from Fig. 2.

Calculations were conducted for the all electron (AE) BEB curves and the ECP BEB curves using QEC 1.2.⁵⁹ Series of molecules were chosen, which were designed to show the increasing importance of using ECPs. The molecular geometries were taken from the NIST Computational Chemistry Comparison and Benchmark database⁴⁵ (CCCBDB) or, when unavailable, were estimated using similar known structures. The geometries were then optimized at the HF/cc-pVDZ level using Molpro2019 through the QEC optimization facility. For molecules containing heavy elements meaning only ECP calculations could be conducted, geometries were optimized using the same ECP basis set that was used to obtain molecular orbital parameters for the BEB calculation. HF orbitals were used in all calculations. For the AE BEB curves, all calculations were completed to a cc-pVDZ level with the exception of SnH₄, which was studied at the STO-3G level due to the availability of basis sets. For the ECP BEB curves, cc-pVDZ was used for hydrogen atoms and the appropriate ECP for all other atoms with Z > 3. In the case of the atoms used in this work, QEC applied the appropriate ECPnMWB basis set, where “n” represents the number of core electrons, to the selected atoms. See caption headings for more details. Where possible, these curves are compared to experimental absolute values. Results from other computational methods have also been included. These figures are shown in Figs. 3–15.

A file giving details of the BEB calculation in the form of orbital energies and kinetic energy integrals is given in the supplementary material. In general, we consistently see ∼2% difference in the orbital energies obtained using an ECP compared to an AE calculation, with the ECP giving higher binding energies. There is more variation as expected in the kinetic energy integrals. The most notable effect comes when the outer most valence orbitals in an AE calculation penetrate the core, the case when the HOMO is σ for instance, or when n the principle quantum number is high. In such cases, there is a significant reduction in the kinetic energy integrals when an ECP is used, which can be as much as 60% lower. When the HOMO orbital is π or in the case of low n such as for species containing F, the HOMO often has a node at the atomic center and so does not significantly penetrate the core region of the atom. This leads to little difference, or slight increase is seen in the kinetic energy integrals of the outermost valence electrons. As the outermost valence orbitals dominate the electron impact ionization cross sections, this can lead to very similar or slightly reduced BEB cross sections for ECP BEB than AE BEB.

The use of an ECP is expected to improve the BEB curve against the AE BEB with increasing effect for heavy atoms (Z > 10). However, a crucial limitation of the BEB method is that multiple ionization events are not included. Therefore, despite an improvement being expected to be observed, the ECP BEB cross sections will be lower than the experimental results when multiple ionization events contribute significantly to the total ionization cross section. Furthermore, at high energies, core ionization can occur, which will be considered by the AE BEB curves but not by the ECP BEB curves.

The result from the calculated polarizabilities compared to experimental values can be seen in Table I. As this work was to investigate whether the αBEB method results in a constant improvement to the total ionization cross section, the values used for the polarizabilities were experimental. The investigation into the best method for polarizability calculations gives insight into the level of accuracy that could be expected if computational methods were used.

The relationship between the dipole polarizabilities and the total ionization cross section maximum has been shown by a plot of experimental polarizabilities^45,47–49 against experimental total ionization cross section maximum^50–58 for the nine molecules, which are shown in Fig. 2. The gradient of a linear fit through the origin was used to predict an absolute value for the maximum of the total ionization cross section.^25,27,43 The difference between this predicted maximum and the maximum of the BEB curve can been used as a scaling factor, and thus, the BEB curves were corrected to pass through this maximum. The application of Eq. (2) to the ECP BEB curves to produce the αBEB curves can be seen in Figs. 6–14. Additionally, it is evident from Fig. 2 that there is a linear relationship between ECP BEB curves and the experimental polarizabilities for the 12 molecules in this study.^45,47–49

The correlation coefficient was 0.99 and 0.97 for the ECP BEB data points and the experimental data points, respectively, showing a strong linear relationship for both datasets. The gradient of the experimental data was determined to be 1.563 Å⁻¹. This is similar to the value of 1.478 Å⁻¹ determined by Bull et al.,⁶⁰ who used 63 medium sized organic and halocarbon species and obtained a correlation coefficient of 0.98. However, we note that some molecules are significantly far from this linear relationship, particularly BCl₃, CF₄, PH₃, and SiH₄. This, combined with issues in computing polarizabilities with the accuracy required to improve the BEB predictions, means that we do not pursue the αBEB method as a black box procedure.

Electron impact ionization cross sections for the series, CF₄, CCl₄, and CBr₄, can be seen in Figs. 3–5, respectively. Additionally, cross sections for the final two molecules in this series, CI₄ and CAt₄, are shown in Fig. 15 as only ECP BEB calculation could be completed due to the computational expense of the AE BEB calculations.

It can be seen throughout the series that the ECP BEB curve results in a cross section that is larger than the AE BEB. An exception to this is the slight decrease observed for CF₄, for which the kinetic energy integrals for the outermost valence orbitals are slightly larger than the comparative AE calculation. However, as both carbon and fluorine are light elements, a large change would not be expected, and in this case, the outermost valence orbitals have nodes at the atomic centers. In cases where experimental values are available, both the AE BEB and ECP BEB can be seen to underestimate the cross section in comparison to the experiment. For CBr₄, the ECP BEB curve is in much closer agreement to the computational results of Naghma et al.,⁶¹ who used a SCOP method, compared to the AE calculation.

Experimental polarizability values are available for the lightest two molecules (CF₄ and CCl₄) in the series, and hence, αBEB curves were calculated. For CF₄, this results in a cross section that is further from experimental values, whereas in the case of CCl₄, the cross section matches the experimental values of Lindsay et al.⁶² within the experimental uncertainty at almost all energies.

Results for the series of CH₄, SiH₄, GeH₄, and SnH₄ are shown in Figs. 6–8. The final molecule in this series, PbH₄, can be seen in Fig. 15 as only ECP calculations could be completed.

Once again, the expected trends can be observed where the ECP BEB curves are an improvement over the AE BEB, especially for molecules that contain heavy atoms. In the case of the smaller two molecules of the series, CH₄ and SiH₄, the use of an ECP shows a small change compared to the AE curve and brings the cross sections closer to experimental results. For GeH₄, seen in Fig. 1, only one experimental value could be found.⁴⁰ For comparison, data from the NIST database²⁹ have also been included. It can be seen, as expected, that the ECP BEB values are larger than the AE ones bringing them closer to experiment. However, for SnH₄, the AE BEB curve has a cross section that is unexpectedly larger. This irregularity is due to the STO-3G basis set used for the AE calculation, which is a small basis set that overestimates BEB cross sections, largely due to the orbital energies being significantly less bound in the AE calculation. This small basis set has also lead to differences at the threshold that is not observed for the other molecules, as well as changes to the shape of the cross section, as the peak of the ionization cross section is shifted to higher energy in comparison to the ECP curve.

Two of the molecules, CH₄ and SiH₄, also had αBEB curves calculated. In the case of CH₄, this method brings the cross section to be within the experimental uncertainty of the recommended values by Song et al.²⁷ However, for the latter molecule, the use of αBEB results in a large overestimate of the cross section and, hence, does not provide an improvement over either the AE BEB or ECP BEB results.

The final series considered is PH₃, PF₃, and PCl₃, shown in Figs. 9–11. Unlike the previous two series, all molecules in this series contain at least one heavy atom with the final molecules, PCl₃, containing entirely heavy elements. This could lead to the assumption that the improvement seen due to an ECP for PH₃ and PF₃ would be similar, much like the small change that can be observed in both CH₄ and CF₄. However, for PF₃, the AE BEB curve and ECP BEB curve are almost identical at all energies, and hence, a larger difference can be observed within PH₃. The minute difference observed for PF₃ can be accounted for by the similarity in the kinetic energy integrals for this molecule for AE and ECP basis sets being very similar. This is again due to the nodal structure of the outermost valence orbitals at the atomic centers.

The only experimental values found for this series were for PH₃. It can be seen that the AE BEB curve matches the experimental values of Märk and Egger.⁵⁴ Furthermore, computational results by Kumar,⁶⁸ who used the Jain and Khare method, and the SCOP results by Vindokumar et al.⁶⁹ are also included. The AE curve matches experiment better than the ECP BEB curve does. One set of theoretical data by Kumar,⁷⁰ who used the modified Jain and Khare method, for PF₃ is included. It should be noted that Katyal et al.⁷¹ also used the modified Jain and Khare method and, thus, obtained identical results. For PCl₃, the results are as expected, the ECP BEB has a larger cross section than the AE.

The αBEB method was tested for all three of the molecules in this series. The unpredictability of this method along with the lack of experimental results makes the accuracy difficult to evaluate. However, the same trend can be observed as previously, where the αBEB curves are larger for all three molecules than the AE BEB and the ECP BEB curves.

The final four molecules studied in this work are SiCl₄, BCl₃, CH₃Br, and CF₃I. As the latter molecule contains an iodine atom, only the ECP BEB curve could be calculated, and hence, this molecule is included in Fig. 15 along with CAt₄, CI₄, and PbH₄. The other three final molecules are shown in Figs. 12 and 13.

All constituent atoms in SiCl₄ are considered as heavy atoms, much like for PCl₃, but this time absolute experimental cross sections were measured by Basner et al.,⁵⁵ and hence, a comparison could be conducted. Figure 12 shows that all BEB methods tried here underestimate the electron impact ionization cross section in comparison to the experimental data. The use of an ECP results in cross sections closer to experimental values, which is further improved upon by using the αBEB method. For CH₃Br in Fig. 13, the ECP BEB can be seen to be an improvement over the AE BEB bringing the cross section to almost within an experimental uncertainty. However, the αBEB method results in an overestimation in the cross section for almost all energies. Despite absolute experimental values being found for BCl₃, for energies below 60 eV, these values do not cover the cross section maximum, and thus, the accuracy, magnitude, and shape of the curve are difficult to analyze. Christophorou and Olthoff⁷² suggest that the values by Jiao et al.⁷³ should be taken as a lower limit and that significant discrepancies exist among available data that need further investigation. However, the expected trends are still observed, where the ECP BEB curve produces a cross section larger than the AE BEB, and the αBEB curve is larger than the ECP BEB curve.

For the four molecules shown in Fig. 15, no experimental values could be found, and only ECP BEB calculations could be completed. As no experimental polarizabilities were available either, the αBEB curves were also not calculated. Despite this, it can be seen that the molecule CAt₄ has the largest cross section and is by far the heaviest molecule; this is followed by CI₄, which is the second largest molecule.

We test various methods of computing the electron impact ionization cross section within the BEB framework with a view to implementing an automated and predictive procedure within the QEC expert system. A consistent improvement can be seen through the use of an effective core potential (ECP) against the all-electron (AE) BEB curves for molecules containing at least one heavy atom. For molecules containing light atoms only, there is either little difference in the resulting cross sections using an ECP compared to an AE calculation or the ECP BEB slightly smaller. Therefore, an AE BEB calculation should be recommended for molecules containing atoms from the first two rows of the periodic table.

Scaling the maximum of the BEB curve according to the polarizability (αBEB) proved to be inconsistent meaning that this method in its current form is not suitable for the implementation as part of a black box procedure. We have, therefore, implemented only the ECP option in the new release of BEB in the QEC software.

Another important consideration following electron impact ionization is the fragmentation pattern of the resulting molecular ion. Hamilton et al.⁷⁵ provided a method of estimating these patterns based on data from mass spectrometry. However, this procedure cannot be applied for species, such as most radicals, for which such data are not available. We are currently working on a general procedure to automatically estimate fragmentation patterns for use in plasma modeling.

See the supplementary material for the details of the BEB calculation in the form of orbital energies and kinetic energy integrals.

This work was supported by STFC (Grant No. ST/R005133/1).

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