Electron impact ionisation cross sections of iron hydrogen clusters

We computed electron impact ionisation cross sections (EICSs) of iron hydrogen clusters, FeHn with n = 1,2,...,10, from the ionisation threshold to 10 keV using the Deutsch-Märk (DM) and the binary-encounter-Bethe (BEB) formalisms. The maxima of the cross sections for the iron hydrogen clusters range from 6.13 × 10-16 cm2 at 60 eV to 8.76 × 10-16 cm2 at 76 eV for BEB-AE (BEB method based on quantum-chemical data from all-electron basis sets) calculations, from 4.15 × 10-16 cm2 at 77 eV to 7.61 × 10-16 cm2 at 80 eV for BEB-ECP (BEB method based on quantum-chemical data from effective-core potentials for inner-core electrons) calculations and from 2.49 × 10-16 cm2 at 43.5 eV to 7.04 × 10-16 cm2 at 51 eV for the DM method. Cross sections calculated via the BEB method are substantially higher than the ones obtained via the DM method, up to a factor of about two for FeH and FeH2. The formation of Fe-H bonds depopulates the iron 4s orbital, causing significantly lower cross sections for the small iron hydrides compared to atomic iron. Both the DM and BEB cross sections can be fitted perfectly against a simple expression used in modelling and simulation codes in the framework of nuclear fusion research. The energetics of the iron hydrogen clusters change substantially when exact exchange is present in the density functional, while the cluster geometries do not depend on this choice.


Introduction
Plasma-wall interaction (PWI) remains one of the key issues in nuclear fusion research. In nuclear fusion devices, such as the JET tokamak or the International Thermonuclear Experimental Reactor (ITER; presently under construction), first-wall materials will be directly exposed to plasma components. In ITER, the first-wall is envisaged to be coated with beryllium (at the main wall) and tungsten (in the divertor) [1]. JET has demonstrated successful plasma operation using this combination of materials which strongly supports this choice for ITER [2]. Taking a longer term perspective, such as in the fusion program DEMO and beyond it in industrial applications of nuclear fusion, it would be far better, however, to avoid the highly toxic and hence difficult to handle beryllium. The use of stainless steel (i.e. the Eurofer steel envisaged for DEMO [3,4]) for some portions of the main wall may then come into consideration. Erosion of first-wall materials is a consequence of the impact of hydrogen and its isotopes as main constituents of the hot plasma [5,6]. Besides the formation of gas-phase atomic species in various charge states, also di-and polyatomic molecular species are expected to be formed via PWI processes. These compounds may profoundly disturb the fusion plasma and may also lead to unfavourable re-deposition of materials and composites in other areas of the vessel [7][8][9][10]. Detailed knowledge and quantification of interactions between atoms, molecules and the plasma as well as of the transport of impurities is thus of considerable interest for modelling and simulation of fusion plasmas [11]. Collisions of atoms and molecules with electrons are an important example of such processes. They are mainly characterised by the respective electron-impact ionisation cross sections (EICSs). Their understanding is especially important for modelling the plasma energy balance. Apart from magnetic confinement fusion, EICS data also plays a role in astrophysics and in a variety of other applications, such as low-temperature processing plasmas, gas discharges, and in chemical analysis [12].
Generally speaking, FeH n molecules and clusters are rather unstable, compared to other metal hydrides. Iron hydride, FeH, has only been identified in the gaseous state and in matrices. It is of considerable astronomical importance and has been subjected both to experimental and theoretical investigations [13]. Its so-called Wigner-Ford band, an IR peak at 990 nm, has been found in the spectrum of the sun. FeH 2 has a 5 Δ g ground state. The frequency of its vibrational bending mode is 226 cm −1 , characteristic of a very floppy angular bending [14], while the asymmetric Fe-H stretching frequency at 1674 cm −1 indicates bonds of medium strength [15]. In addition, FeH 2 is normally only stable as a diluted gas or in matrices, although recently [16] a solid high-pressure phase with the stoichiometry of FeH 2 has been synthesised. All other FeH n clusters from FeH 3 onwards can be thought to consist of FeH or FeH 2 with H 2 molecules bound to them. In the following, we simply call all FeH n assemblies iron hydrogen clusters. Takahashi et al. [17] report data on a number of these clusters from plane-wave DFT calculations. Their findings on geometries agree mostly with ours, although there are small discrepancies (see Sect. 3.1).
During the past few decades, a number of semiempirical methods that typically use quantum-chemically calculated electronic structure information as input have been developed in order to derive absolute EICSs for various molecules. Their accuracy is usually in the same range as corresponding experimental data. Among them, the most-widely used methods are the binary-encounter-Bethe (BEB) theory of Kim and Rudd [18], Kim et al. [19] and the Deutsch-Märk (DM) formalism [20]. These methods have been successfully applied to atoms, molecules, clusters, ions and radicals [21].
In the context of fusion-relevant species, EICSs were reported earlier for beryllium [22,23], its hydrides [24], tungsten and its oxides [25,26] and beryllium-tungsten clusters [27]. In this work we report calculated EICSs using both the BEB and the DM method for neutral iron-hydrogen clusters, in particular for FeH n with n = 0, 1, 2, . . . , 10 compounds, since iron is the main component of steels. To the best of our knowledge there is no such cross section data available for iron-hydrogen clusters up to now. We also report parameters obtained by fitting the calculated cross sections to an expression commonly used in codes modelling the impurity transport in fusion edge plasmas such as ERO [28][29][30].

The DM formalism
The DM formalism was originally developed as an easy-touse semi-empirical approach for the calculation of EICSs of atoms in their electronic ground-state from threshold to about 100 eV [20]. In its most recent variant [21,31], the total single electron-impact ionisation cross section σ of an atom is expressed as: where r nl is the radius of maximum radial density of the atomic sub-shell characterised by quantum numbers n and l (as listed in column 1 in the tables of Desclaux [32]) and ξ nl is the number of electrons in that sub-shell. The sum extends over all atomic sub-shells labeled by n and l.
The g nl are weighting factors, which were originally determined by a fitting procedure [33,34] using reliable experimental cross section data for a few selected atoms, for which the accuracy of the reported rate is in the range of 7-15%. The reduced energy u is given by u = E/E nl , where E refers to the incident energy of the electrons and E nl denotes the ionisation energy of the sub-shell characterised by n and l. The energy-dependent quantities b (q) nl (u) were introduced in an effort to merge the highenergy region of the ionisation cross section, which follows the Born-Bethe approximation [35], with the DM formula of the cross sections in the regime of low impact energies.
where the four constants A 1 , A 2 , A 3 and p were determined, together with c nl , from reliably measured cross sections for the various values of n and l. The superscript q refers to the number of electrons in the (n, l)-th sub-shell and allows the possibility of using slightly different func- nl depending on the number of electrons in the respective sub-shell. At high impact energies u goes to infinity, the first term in equation (2) goes to zero and b (q) nl (u) becomes a constant ensuring the high-energy behaviour predicted by the Born-Bethe theory [35].
The DM formalism has been extended to the calculation of EICSs of atoms in excited states, molecules and free radicals, atomic and molecular ions, and clusters [21]. For the calculation of the EICS of a molecule, a population analysis [36,37] must be carried out to obtain the weights with which the atomic orbitals of the constituent atoms contribute to each occupied molecular orbital. These weights are obtained from the molecular orbital coefficients after a transformation employing the overlap matrix in order to correct for the non-orthogonality of the atomic basis functions.

The BEB method
The BEB model [19] was derived from the binaryencounter-dipole model [18] by replacing the df /dE term for the continuum dipole oscillator strengths with a simpler form. Thus, a modified form of the Mott cross section, together with the asymptotic form of the Bethe theory describing the electron-impact ionisation of an atom, was combined into an expression for the cross section of each molecular orbital: where t = T /B, u = U/B, S = 4πa 2 0 N R 2 /B 2 , a 0 denotes the Bohr radius (0.5292Å), R is the Rydberg energy (13.6057 eV), and T denotes the incident electron energy. N , B and U are the electron occupation number, the binding energy (ionisation energy), and the average kinetic energy of the respective molecular orbital, respectively. In the BEB model, the total cross section, similarly to the DM method, is then obtained by summation over the cross sections for all molecular orbitals.
The cross section formula given by equation (3) has experienced several modifications over the years for atoms and molecules involving orbitals corresponding to principal quantum numbers larger than two [38][39][40][41][42][43]. The BEB cross section that includes this modification, which became known also as "acceleration correction" [44], reads: where for atoms n = pqn, with pqn denoting the principle quantum number, for orbitals with pqn 3, n = 1 for orbitals with pqn = 1, 2, and for molecules n = pqn if the Mulliken population analysis of the respective molecular orbital yields that the component of a specific atomic orbital with pqn 3 is dominant, i.e. its contribution is larger than 50%, and n = 1 otherwise. Also for singly charged molecular ions n = 2 is taken, which does not, however, apply to this work.
The quantum chemical data needed to calculate EICSs are normally derived from all-electron calculations. For heavy elements and molecules containing them valence-shell-only calculations using effective-core potentials (ECPs) [45] can be used. This facilitates the quantum chemical calculations and allows the incorporation of relativistic effects. Due to the lack of inner radial nodes of the pseudo-valence orbitals, their kinetic energies are lower than normal and equation (3) can be used to determine the BEB cross section [46]. In this work we compare all-electron and valence-only methods to obtain BEB cross sections. When equation (4) is used together with orbital and kinetic energies from all-electron quantum chemical calculations, we refer to the method as the BEB-AE method for the remainder of this work. Due to the use of ECPs in determining the orbital and kinetic energies, we refer to the second method as BEB-ECP method for the remainder of this work. The BEB-ECP method has earlier been recommended over the all-electron BEB-AE method for molecules that contain heavy (atomic number Z > 10) atoms [47].

Quantum chemical calculations
Recently, globally optimised structures of FeH n with n = 0, 1, 2, . . . , 10 were reported [17]. We used these as starting geometries that were further optimised employing the PBE [48,49], PBE0 [50] and B3LYP [51] density functionals, all in conjunction with the Def2-TZVP basis set [52,53]. This allows exploration of the effect of inclusion of exact exchange on the geometry and energetics of these clusters. The binding energies, E BE of the iron hydrogen clusters were determined according to: where E(x) denotes the energy of compound x including the zero-point vibrational energy. Occupation, binding energy and average kinetic energy for each molecular orbital as required for the calculation of the BEB-AE and BEB-ECP cross sections (see Sect. 2.2) were calculated at the HF/Def2-TZVP and HF/CEP-4G levels of theory, respectively, using the geometries obtained with B3LYP/Def2-TZVP (concerning geometric parameters all three employed density functionals yield very similar results, see Sect. 3.1). The orbital populations required for the DM formalism were derived from HF calculations in conjunction with the minimal CEP-4G basis set [54][55][56]. Orbital energies for the outermost valence electrons were calculated with the OVGF method and the Def2-TZVP basis set [57].
All calculations were performed with the Gaussian 09 software [58].

Analytical expression of the EICSs
We fitted the cross sections to an expression that resembles the one used in the ERO code [28][29][30] which is used for impurity transport simulations in fusion edge plasmas. The fitting expression is given by: (6) Here, the cross section σ is expressed in 10 −16 cm 2 , the incident electron energy E and the threshold energy (first ionisation energy) E t are both expressed in eV, and the fit parameter a 1 is expressed in 10 −16 cm 2 eV. The fit parameters a 2 , a 3 and a 4 are dimensionless.

Structure and energetics of the iron hydrogen clusters
The structures of the iron hydrogen clusters considered in this work are depicted in Figure 1. In Table 1, we summarise the binding energies and average Fe-H bond lengths of the various clusters. We also include the values of the respective quantities reported by Takahashi et al. [17]. Takahashi et al. [17] used a plane-wave basis for their calculations, which can explain the slight discrepancies between their and our PBE values. In particular, at the PBE/Def2-TZVP level of theory we obtain a small negative binding energy for the FeH molecule and the increase of binding energies with increasing number of hydrogen atoms in the cluster is less pronounced in our results. In the case of the hybrid functionals, the binding energy hardly increases with the number of hydrogens and most of them are in the range 0.8 to 1.5 eV for PBE0 and 0.6 to 1.0 eV for B3LYP. Exceptions are the less stable FeH, FeH 4 and FeH 5 clusters in both cases. Whereas FeH and FeH 4 are slightly stable, FeH 5 , in contrast, is predicted to be (thermodynamically) unstable by both methods (at zero temperature). The variation of the binding energy with the number of hydrogens follows the same trend for both functionals with quantitative differences within 0.5 eV.
The geometries of the clusters are very similar for the three density functionals, indicating in particular that inclusion of exact exchange plays no significant role. This can also be seen from the average Fe-H bond lengths given in Table 1. They vary by mostly 0.05Å for the three methods. PBE0 yields slightly longer bonds than PBE, with differences in average bond length within 0.02Å. B3LYP yields slightly longer bond lengths than PBE0 and the differences in average bond length are within 0.03Å.

EICS of Fe
In Figure 2, we plot the EICSs of iron calculated in this work, together with those from former theoretical and experimental studies. We note that our DM cross section   (thick black line in Fig. 2) with a maximum of 4.58 × 10 −16 cm 2 at 33 eV is in reasonable agreement with the experimental cross section obtained by Shah et al. [59] (full black circles in Fig. 2; maximum of 4.08 × 10 −16 cm 2 at 35 eV). This experimental cross section appears to be more reliable than the cross section obtained by Freund et al. [60] (open black squares in Fig. 2) as the latter shows a non-zero cross section below the ionisation energy of iron in its electronic ground state, i.e. 7.92 eV [61]. The substantially higher maximum of 5.34 × 10 −16 cm 2 at 29 eV in the latter cross section [61] is also an indication that this cross section includes contributions from the ionisation of excited metastable iron [62]. There is a slight difference between our DM cross section and the one obtained by Deutsch et al. [62], which yields a maximum of 4.59 × 10 −16 cm 2 at 25 eV and was also obtained using the DM formalism. However, it is not completely clear which data was used for the atomic orbital energies in their work [62].
In the case of atomic Fe, we used the experimental ionisation energy of 7.92 eV for the electrons in the 4s orbital, while the energies of the five 3d orbitals were obtained by the OVGF method and were averaged to obtain a single value for the 3d orbital. Contributions from all lower lying orbitals were found to be negligible for the cross section. Data from earlier theoretical studies yield higher  Fig. 2) yields by far the highest value for the cross section (maximum at 6.77 × 10 −16 cm 2 at 42.5 eV). The BEB-ECP cross section (thick dashed red line in Fig. 2) is lower and in significantly better agreement with experimental and other theoretical data (maximum at 4.86 × 10 −16 cm 2 at 46.5 eV). The experimental ground state of the iron atom is [Ar] 3d 6 4s 2 and our cross sections have been constructed for this state. HF calculation of the quintetconfiguration results in the 3d 7 4s 1 electronic configuration which is experimentally the first excited electronic state. Such behaviour is not unusual in the case of transition metals with several electronic states close to each other. In order to calculate the cross section for the [Ar] 3d 6 4s 2 state, we simply removed one of the 3d electrons and inserted it into the 4s orbital, for which the experimental ionisation energy was used. The orbital energies of the six remaining electrons in 3d orbitals were left at values calculated with the OVGF method. This procedure amounts to a slight additional approximation. The BEB-AE, BEB-ECP and DM cross sections for the 3d 7 4s 1 configuration of iron are also depicted in Figure 2 (thin black and red lines) and are substantially lower than the respective cross sections for the Fe ground state, yielding 5.65 × 10 −16 cm 2 at 50 eV, 3.72 × 10 −16 cm 2 at 56 eV and 2.69 × 10 −16 cm 2 at 38.5 eV for BEB-AE, BEB-ECP and DM, respectively.
Convergence of the various cross sections at elevated electron impact energies is observed for our DM cross section and the cross sections obtained by Deutsch et al. [62], McGuire [63] and Shah et al. [59] above 500 eV. The cross section reported by Tsipinyuk et al. [44] approaches the experimental cross section of Freund et al. [60] above 100 eV. The cross section obtained by Goswami et al. [26] seems to approach our BEB cross sections above 1000 eV.

EICSs of the iron hydrogen clusters
In Figure 3, we depict the BEB and DM cross sections obtained for the iron hydrogen clusters FeH n with n = 1, 2, . . . , 10. Tabulated data of the cross sections are also available as supplementary material. In Section S4 of the supplementary information, we give the respective cross section maxima and their respective energies together with the parameters obtained by fitting the cross section data from the DM, BEB-AE and BEB-ECP calculations to equation (6). In summary, the maxima of the cross sections for the iron hydrogen clusters range from: 6.13 × 10 −16 cm 2 at 60 eV to 8.76 × 10 −16 cm 2 at 76 eV, 4.15 × 10 −16 cm 2 at 77 eV to 7.61 × 10 −16 cm 2 at 80 eV, 2.49 × 10 −16 cm 2 at 43.5 eV to 7.04 × 10 −16 cm 2 at 51 eV for the BEB-AE, BEB-ECP and DM methods, respectively. Again we note that the BEB-AE method yields substantially higher cross sections than the DM method, whereas the BEB-ECP yields better agreement with the latter. The BEB-AE cross section maxima for FeH and FeH 2 are higher by more than a factor two than the respective DM values. In addition, the BEB-ECP cross sections are substantially larger than the DM cross sections, up to almost a factor of two for the FeH, FeH 2 and FeH 3 clusters. For the larger clusters the cross section maxima are more similar. For the FeH 10 cluster, BEB-AE and BEB-ECP cross sections are, however, still larger by about 24% and 8%, respectively, than the DM cross section. It should be noted that discrepancies of up to 50% between the results of different methods, as well as similar discrepancies between calculated and experimentally determined ionisation cross sections are not unusual [19,21,34]. Nevertheless, deviations of a factor of two or more appear surprising. To our knowledge, only for atomic tungsten has such a deviation (of about a factor of two) been reported previously [62].
It is interesting to note that especially the cross sections of the small iron hydrogen clusters are substantially smaller than those obtained for atomic iron. This is the   case for both the DM and the BEB methods. Inspection of the orbital populations shows a low population of the Fe 4s orbitals and a significant amount of electron transfer from iron to the hydrogen atoms. This is in line with the results of an earlier study on the stability and bonding mechanisms of iron hydrogen clusters [17]. In fact, the cross sections for FeH are rather similar to those of atomic iron in the 3d 7 4s 1 configuration also indicating the depopulation of the 4s orbital upon Fe-H bonding. Since the contribution of the 4s electrons dominates (by almost 80%) the DM cross section of atomic iron due to the comparably large value of the maximum radial density of the 4s orbital, depopulation of the 4s orbital upon Fe-H bonding can significantly reduce the respective cross section. The BEB formalism captures the contributions from the populations of different atomic orbitals only implicitly via potential and kinetic energies of the molecular orbitals and lacks any further geometric parameters. Due to the large size of the 4s orbital, its kinetic energy is smaller than the kinetic energies of the Fe 3d orbitals and its depopulation causes a lowering of the cross section as can be seen from the BEB formula. As such, both the BEB and DM methods succeed in reproducing the same general trend.
Upon increasing the number of hydrogen atoms in the cluster, the DM cross sections increase continuously. Two BEB cross sections show a different behaviour. The maxima of the BEB-AE and BEB-ECP cross sections for FeH 3 and for FeH 5 are higher than the maxima for FeH 2 , FeH 4 and FeH 6 .

Conclusion
We calculated electron impact ionisation cross sections (EICSs) of iron hydrogen clusters, FeH n with n = 1, 2, . . . , 10, from the ionisation threshold to 10 keV using the Deutsch-Märk (DM) and the binary-encounter-Bethe (BEB) formalisms. We also compared the effect of using either all-electron basis sets or effective-core potentials for inner-core electrons in the quantum chemical calculations necessary to obtain the orbital and kinetic energies required for the BEB method on the resulting BEB-AE or BEB-ECP cross sections, respectively. The maxima of the cross sections for the iron hydrogen clusters range from: 6.13 × 10 −16 cm 2 at 60 eV to 8.76 × 10 −16 cm 2 at 76 eV, 4.15 × 10 −16 cm 2 at 77 eV to 7.61 × 10 −16 cm 2 at 80 eV, 2.49 × 10 −16 cm 2 at 43.5 eV to 7.04 × 10 −16 cm 2 at 51 eV for the BEB-AE, BEB-ECP and DM methods, respectively. Both BEB approaches yield cross sections substantially higher than those obtained via the DM method, most pronounced for the smallest FeH and FeH 2 clusters with deviations of more than a factor of two for BEB-AE and almost a factor of two for BEB-ECP. Generally, the agreement between BEB-ECP and DM is better than between BEB-AE and DM. Together with the comparison of the obtained cross sections for atomic iron with experimental and other theoretical data, this finding agrees with a former recommendation of the use of the BEB-ECP approach for molecules containing heavy elements [47]. Upon Fe-H bonding, we note a substantial depopulation of the iron 4s orbitals underlying significantly lower cross sections for small iron hydrogen cross sections than for the cross sections of atomic iron, but similar to those obtained for excited iron in its 3d 7 4s 1 electronic configuration. Both the DM and BEB cross sections could be fitted perfectly against a simple expression used in modelling and simulation codes in the framework of nuclear fusion research. Moreover, we analysed structures and energies of the iron hydrogen clusters, which revealed a significant effect of the inclusion of exact exchange in the density functional for the energetics while the cluster geometries were insensitive to it.