Effects of quantum chemistry models for bound electrons on positron annihilation spectra for atoms and small molecules

The Doppler-shift spectra of the γ-rays from positron annihilation in molecules were determined by using the momentum distribution of the annihilation electron–positron pair. The effect of the positron wavefunction on spectra was analysed in a recent paper (Green et al 2012 New J. Phys. 14 035021). In this companion paper, we focus on the dominant contribution to the spectra, which arises from the momenta of the bound electrons. In particular, we use computational quantum chemistry models (Hartree–Fock with two basis sets and density functional theory (DFT)) to calculate the wavefunctions of the bound electrons. Numerical results are presented for noble gases and small molecules such as H2, N2, O2, CH4 and CF4. The calculations reveal relatively small effects on the Doppler-shift spectra from the level of inclusion of electron correlation energy in the models. For atoms, the difference in the full-width at half-maximum of the spectra obtained using the Hartree–Fock and DFT models does not exceed 2%. For molecules the difference can be much larger, reaching 8% for some molecular orbitals. These results indicate that the predicted positron annihilation spectra for molecules are generally more sensitive to inclusion of electron correlation energies in the quantum chemistry model than the spectra for atoms are.

3 correlation accounted for when using different computational quantum chemistry methods. This is particularly important for molecules, where the single-centre quantum mechanical methods developed for atoms are inapplicable, and the overall understanding is poor.
The aims of the study described in this paper are (i) to investigate the accuracy of boundelectron wavefunctions (e.g. the effects of specific computational quantum-mechanical models on annihilation γ-ray spectra), and (ii) to investigate the contributions of individual orbitals of the atoms or molecules to the γ-ray spectra. This paper provides useful information towards the overall goal of developing accurate physical models and an increased understanding of γ-ray annihilation spectra for molecules.

Methods and computational details
In an independent-particle approximation, the photon spectrum is determined by the annihilation amplitude (see, e.g., [4,18]) where ψ i (r) is the wavefunction of the electron in orbital i of the target in the ground state, φ k (r) is the wavefunction of the incident positron with momentum k, and P is the total momentum of the annihilation photons. These electron and positron orbitals, which are calculated quantum mechanically in the coordinate representation, are then Fourier transformed using equation (1). The probability distribution function of the photon momentum P in two-photon annihilation is given by where r 0 is the classical electron radius, c is the light speed and the subscript i refers to a particular hole in the final state of the electronic system. If the total momentum of the two photons is not zero, then the energy of each of the photons is Doppler-shifted. The shift of the photon energy from the centre of the spectral line (energy mc 2 = 511 keV) is given by ε = c P cos θ/2, where θ is the angle between the direction of the photon and the velocity of the electron-positron pair. As a result, the photon energy spectrum for a spherically averaged target is given by At large separations the positron is described by a plane wave, and for the low positron momentum k in the vicinity of the target φ k ≈ e ik·r ≈ 1, which we term the LEPWP approximation [12]. In this approximation, the γ-ray spectrum is determined by the bound electron contribution, which can be obtained by solving the electron Schrödinger equation for the target quantum mechanically. Equation (3) gives the annihilation spectrum for a particular orbital i. The total spectrum is found by summing over all electronic orbitals w(ε) = i w i (ε).
As mentioned in the introduction, the annihilation spectra obtained in the LEPWP approximation are broader than those measured experimentally [12][13][14]. As an example, shown in figure 1 are the annihilation spectra obtained from equations (1)-(3) in the LEPWP approximation for the valence orbitals of N 2 using the Hartree-Fock/triplet zeta with valence polarized (HF/TZVP) model (see below), together with the total spectrum for these valence orbitals. As shown in the figure, when the measured spectrum is represented by a Gaussian fit (see [8] for details) and broadened artificially by a factor of 1.44, it compares well with the model calculation. The factor of 1.44 arises from a systematic effect due to the high-momentum electron contributions that come from a range of distances close to the nuclei and would be suppressed if a proper positron wavefunction were used. Although this suppression is important for explaining the experimental spectra, it is relatively featureless [13]. Leaving this effect aside, we test here how sensitive the calculated spectra are to the models used for computing the electronic wavefunctions.
If the γ-ray annihilation spectra for bound-electron systems such as atoms and molecules are largely determined by the atomic or molecular wavefunctions, a key issue is understanding the role that the accuracy of the electronic wavefunctions plays in determining the γ-ray annihilation spectra. This paper focuses on exploring the effects of the bound-electron wavefunctions in the LEPWP approximation. The accuracy of the atomic or molecular wavefunctions depends on the quantum mechanical model employed, i.e. the theory and basis set. The most widely used modern ab initio quantum chemistry methods include the HF and post-HF (e.g. many-body theory) methods, as well as density functional theory (DFT) methods with models such as Becke's three-parameter fit for the exchange energy [19][20][21][22][23] and the Lee-Yang-Parr form for the correlation energy (B3LYP) [22].
The HF model contains the exact exchange energy (V x ) but no electron correlation energy (V c ). In the ab initio models, if the basis set is sufficiently complete, more electron correlation energies are included in the post-HF models, and more accurate total energies are produced. Alternatively, the DFT functionals, such as B3LYP, contain various exchange-correlation functionals (V xc ). However, the degree of inclusion of exchange and correlation energies in many DFT-based methods is not precisely known. A significant difference between the post-HF models and the DFT models is that the post-HF models are energy focused, whereas the 5 DFT models are more wavefunction (density) focused. As a result, the post-HF models provide more accurate energetic properties, such as orbitals energies, without taking into consideration wavefunctions, so that the wavefunctions of the target remain at the HF level. The DFT models are significantly different: their improvement relative to the HF model arises from a more accurate treatment of the density (wavefunction) rather than the energetics. As a result, the simulated γ-ray annihilation spectra of a target, which depends on the electron wavefunctions, will be the same when using either HF or post-HF models, but will be different from those simulated using the DFT models.
The choice of basis set is also important for the development of an accurate model. Usually, basis sets with larger sizes lead to more accurate results but at a significant computational cost [25]. The present study uses two popular basis sets, namely the TZVP basis set [12] and the 6-311+ + G * * basis set [20], to examine the effects of the basis set on the calculated Dopplershift spectra. These basis functions are linear combinations of atomic wavefunctions with the coefficients and exponents obtained by optimization. As the positron and electron interaction process can result in strong polarization in the valence space, both basis sets employed in this study include polarization functions that would allow the orbitals to change shape as a result of the positron-electron interactions. (This will be especially important in the future for going beyond the LEPWP approximation.) Diffuse functions allow one to describe orbitals with large radial extent. The small exponents in these basis set functions ensure that the electrons are still associated with the nucleus when they are far away from the atom or molecule. As in Gaussiantype basis functions, the small exponents in the basis functions represent properties far from the nuclei and the larger exponents in the basis functions represent properties in the core region. The TZVP basis set is augmented using polarization functions, but includes no diffuse functions, whereas the 6-311+ + G * * basis set includes both polarization functions and diffuse functions.
The quantum mechanical models HF/6−311+ + G * * , HF/TZVP and B3LYP/TZVP (i.e. 6-311+ + G * * versus TZVPs and HF versus B3LYP) are used here to calculate the γ-ray spectra of positron annihilation in noble gases and small molecules such as H 2 , N 2 , O 2 , CF 4 and CH 4 . All calculations of the electronic wavefunctions are performed using the standard computational chemistry package GAUSSIAN03 [26]. Numerical calculations of the annihilation spectra are performed using a large but finite two-photon cut-off momentum P max = 10 au in equation (3) [12]. This value is sufficient for the valence orbitals of the atoms and molecules studied. The core electron orbitals are characterized by large electron momenta, but in practice they contribute very little to the annihilation spectra. Rather than comparing the shapes of the spectra obtained using different models, we focus here on one parameter that characterizes each spectrum, namely its full-width at half-maximum (FWHM). Table 1 compares the experimental FWHM values ε of the annihilation γ-ray spectra for noble gases with those calculated using different methods in the LEPWP approximation [8]. Also shown in the table are the orbital ionization energies −ε i . The electron wavefunctions of the noble gases, He, Ne, Ar and Kr, were obtained using the HF/6−311+ + G * * , HF/TZVP and B3LYP/TZVP models. The 6-311+ + G * * is a large basis set with both polarized and diffuse functions and is therefore more accurate but more computationally expensive. In contrast, the TZVP basis has only polarized functions, but is less expensive and typically has sufficient accuracy based on our previous electron momentum spectroscopy studies [24]. We investigate Table 1. Comparison of the annihilation γ-ray FWHM spectral linewidths ( ε in keV), calculated in the LEPWP approximation for noble gases using different computational quantum chemistry models for the atomic electrons. Also shown are the orbital energies ε i (in eV). the basis set effect for the HF approximation by comparing the HF/6−311+ + G * * and HF/TZVP results. Note that the results for Xe in table 1 are only given for comparison, since the basis set used for Xe was DGDZVP, as the TZVP basis set is not available for Xe [12]. In addition, Xe 7 is a rather heavy atom, so that wavefunctions obtained using the nonrelativistic Schrödinger equation may suffer from neglect of relativistic effects. As a result, the present study does not include Xe. Moreover, we focus here predominantly on the contributions of the valence electrons since they are expected to provide the dominant contribution to the annihilation process and hence the linewidths. Basis sets affect the accuracy of the wavefunctions of the electrons, particularly for electrons in the outermost valence orbitals, namely the 1s, 2p, 3p and 4p shells for He, Ne, Ar and Kr, respectively. Using the same HF method, the FWHM widths of the Doppler-shift spectra obtained using the 6-311+ + G * * basis set are slightly smaller (and are more in accord with experiments) than those obtained using the TZVP basis set, as shown in table 1. This may be due to the diffuse functions included in the 6-311+ + G * * basis set, which cover a larger region of space to accommodate the valence electrons. However, the improvement using the HF/6−311+ + G * * basis set, as compared to the HF/TZVP basis, is small (not exceeding 0.03 keV). As shown in table 1, for the outermost valence electrons, the difference between the results obtained using the two basis sets decreases as the atomic number increases. For heavy atoms, the two basis sets yield essentially the same results. For this reason, in the present study, only the computationally less expensive TZVP basis set is employed for polyatomic molecules.

Results and discussion
When combined with the TZVP basis set, the HF method, which includes the full electron exchange energy but no electron correlation energy, and the B3LYP method, which includes some electron exchange energy and some electron correlation energy, exhibit only small differences in the simulated line shapes and the FWHM values ε of the γ-ray annihilation spectra. For example, for the FWHM of the 4s shell of Kr, the ab initio HF model produces ε = 2.15 keV, which is only slightly smaller than ε = 2.17 keV given by the B3LYP model. This observation suggests that, at least for atoms, the effects of the electron correlation energy on the widths of the γ-ray annihilation spectra are quite small. As a result, even the simplest HF/TZVP model is able to produce sufficiently accurate orbitals of bound electrons to describe the γ-ray spectra of the noble gas atoms. This was seen in the earlier many-body theory studies [18], which accounted for the positron wavefunction and electron-positron correlation effects, while using the HF electron wavefunctions. For atoms heavier than He, table 1 shows that the FWHM values of the γ-ray spectra for the valence electrons, calculated in the LEPWP approximation, are about 40% greater than the experimental values. Most of this discrepancy is due to the neglect of the effect of positron repulsion from the nuclei as discussed above. At the same time, it is interesting to note that the experimental FWHM values are quite close to those of the γ-ray spectra for the outermost ns orbitals in all the atoms, in agreement with our previous study [12]. Figure 2 shows the relative differences between the calculated FWHM values ( ε) of the γ-ray spectra for the outermost ns orbitals with respect to the HF/TZVP model for the noble gas atoms (except for Xe, which has a different basis set). Here the relative differences shown are ε HF/6−311++G * * / ε HF/TZVP , which compares the 6-311+ + G * * and TZVP basis sets at the HF level of theory, and ε B3LYP/TZVP / ε HF/TZVP , which compares the HF and the B3LYP models using the same basis set of TZVP. For each noble gas atom, the ratios are very close to unity for the quantum mechanical models considered here. However, figure 2 clearly indicates that the effects of the basis sets and models (although small) are the opposite: starting from He, the relative magnitude of ratio due to the basis set effect for the HF model reaches a peak for Ne and then decreases for Ar and decreases further for Kr. At the same time, for the same basis set (TZVP), the model effect (B3LYP relative to HF) is minimum for Ne, then increases for Turning now to molecules, we note that even small molecules possess significantly more complex electronic structures than atoms, due to their multiple centres and chemical bonds. Table 2 reports the calculated orbital-based FWHM values ε of the γ-ray spectra for the diatomic molecules H 2 , N 2 and O 2 , and table 3 shows the results for the small polyatomic molecules CH 4 and CF 4 . For the hydrogen molecule, the models overestimate the experimental ε by approximately 16%. This difference is smaller than that seen for the noble gas atoms, which can be related to the relatively small effect of the nuclear repulsion compared with electron-positron correlations (see [13,14] for a detailed analysis). It should be noted that the inclusion of electron correlation energy via the DFT-based B3LYP model has a weaker effect on the results than the choice of basis set. For example, the ε value for H 2 produced by the HF/TZVP model is 2.02 keV, compared with 2.01 keV from the B3LYP/TZVP model. However, using the same HF theory with different basis sets, we obtain ε = 2.10 keV with the 6-311+ + G * * basis and 2.02 keV with the TZVP basis.
For other diatomic molecules, N 2 and O 2 , different valence orbitals produce quite different ε values. Similar to the atomic cases, it is the molecular orbital (MO) that bonds by s-electrons that gives the ε value closest to the experimental FWHM. As shown in table 2, it is the 3σ g orbital in N 2 and the 2σ g orbital in O 2 . However, unlike the atomic cases, the 2σ g orbital is not the highest-occupied s-electron-dominant MO in O 2 . Perhaps the more obvious difference between a molecule with a few electrons, such as H 2 , and molecules with many electrons, such as O 2 , is that the electron correlation energy plays a more important role for the latter. Inclusion Table 2. Comparison of the orbital-based annihilation γ-ray FWHM spectral linewidths, ε (in keV) for H 2 , N 2 and O 2 , obtained using different models for the molecular electrons. Also shown are the orbital energies ε i (in eV). HF/6−311+ + G * * HF/TZVP B3LYP/TZVP Expt. of the electron correlation energy via the DFT-based B3LYP model reduces the FWHM values ε, whereas the choice of basis set plays a relatively small role in many-electron molecules. For example, the ε value of the 3σ g MO of N 2 produced by the HF model is 2.25 keV when the 6-311+ + G * * basis is employed but 2.30 keV when the TZVP basis is used; whereas the B3LYP/TZVP model gives 2.24 keV. The common features of the results for N 2 and O 2 are that (i) not all valence electrons contribute equally to the γ-ray Doppler shifts; and (ii) the most chemically active electrons in the frontier orbitals do not dominate the annihilation process. Table 3 compares the simulated orbital-based FWHM values for methane (CH 4 ) and fluoromethane (CF 4 ). Note that, due to the computational costs, methane and fluoromethane are calculated only using the HF/TZVP and B3LYP/TZVP models. The 2a 1 electrons of methane give a ε value closest to that of the experimental γ-ray line shape, whereas in fluoromethane, the 2t 2 electron contribution agrees best with the experimental result. The HF/TZVP model gives 1.97 keV for the 2a 1 MO of methane, which is almost the same as that produced by the B3LYP/TZVP model. In fluoromethane, the HF and B3LYP models also produce almost the same ε values of 3.76 and 3.77 keV, respectively, for the same 2t 2 orbital.
It is possible in these molecules that more than one valence orbital determines the shape of the Doppler-shift spectrum. For example, in fluoromethane, the 3a 1 orbital (2.00 keV) and 2t 2 orbital (3.76 keV) are the most probable contributors to the measured Doppler width of 3.04 keV. Determining whether all the valence electrons, a specific group or even a single valence electron dominates the annihilation in a specific molecule requires the inclusion of both the effects of the positron wavefunction and electron-positron correlations. As a partial answer to this question, table 3 shows that some orbitals of the target molecule may be particularly sensitive to the inclusion of electron correlation energies. For example, the HF/TZVP model (5.85 keV FWHM) and the B3LYP/TZVP model (5.37 keV) produce very different Doppler widths for the 4t 2 orbital of CF 4 . It is worth noting that, for both the diatomic and polyatomic molecules examined, the FWHM values for the spectra calculated in the LEPWP approximation for the total valence orbitals are again about 40% greater than the experimental values. Comparison of the FWHM ( ε) values obtained using different models for atoms and molecule reveals that the electron correlation energy effects matter more for the γ-ray spectra in molecules than they do for atoms. The valence-orbital-based relative FWHM values of the γ-ray spectra calculated for N 2 and O 2 and for CH 4 and CF 4 are given in figures 3 and 4, respectively. As shown in figure 3 for N 2 (upper panel) and O 2 (lower panel), the choice of the basis set and the quantum-mechanical model affect the FWHM values for the orbitals differently. While the relative changes in the ε values of the γ-ray spectra of N 2 are both orbital dependent and model dependent, it is clear that in N 2 the values agree to within 1% for the valence orbitals except for the outermost orbital, 3σ g . In O 2 , on the other hand, the electron correlation effects play a more Comparison of basis set effects (triangles) and electron correlation effects (circles) on the FWHM widths ε of γ-ray spectra of positron annihilation for each shell in N 2 (upper panel) and O 2 (lower panel). Triangles represent relative basis set effects: ε HF/6−311++G * * / ε HF/TZVP , and circles represent the relative electron correlation effects: ε B3LYP/TZVP / ε HF/TZVP . important role for the 2σ u and 1 u orbitals. Figure 4 compares the orbital-based linewidths ε calculated using HF and B3LYP models for methane (upper panel) and fluoromethane (lower panel). Interestingly, the model difference is significant in the 1t 1 electrons of methane, whereas in CF 4 , most of the valence electrons do not exhibit significant model-dependent differences except for the 4t 2 electrons. Apart from the model and the orbital-based effects, the numerical scheme [12] implemented in the calculations also plays a role. Although the core shells play a less important role in the γ-ray spectra [8,[12][13][14]18], the assumed numerical cut-off at P max = 10 au in momentum space (used as the integration upper limit in equation (3)) led to a noticeable removal of the contribution to the electron density from the inner shells. For example, the ε values of the total core-shell LEPWP spectra of N 2 and O 2 are approximately 11.98 and 13.56 keV (HF), respectively, which correspond to the momentum P ∼ 7 au. In contrast, the electron density loss for valence electrons due to such a numerical scheme is less significant, as valence orbitals are less extended in momentum space compared with their core counterparts.

Concluding remarks
Orbital-and model-dependent γ-ray Doppler-shift spectra of positron annihilation in noble gases and small diatomic and polyatomic molecules have been computed using modern computational quantum chemistry tools in the LEPWP approximation [12]. The results were obtained using different quantum mechanical models for the bound electrons in the atoms and molecules, namely the HF/6−311+ + G * * , HF/TZVP and B3LYP/TZVP models, which include varying levels of electron exchange and correlation. As expected, the results overestimate the widths (FWHM, ε) of the experimentally measured Doppler-shift spectra. The FWHM values support the understanding that a low-energy positron likely annihilates with electrons in valence orbitals, rather than with core electrons in the cases of atoms as well as small molecules. The present study indicates that modern computational chemistry models are capable of producing accurate bound-electronic orbitals that can be used to predict the γ-ray spectra of positron annihilation in atoms or molecules. Inclusion of electron exchange and correlation energies (such as the DFT-based B3LYP) affects molecules more than atoms, and can change the FWHM value by up to 8% for some orbitals. At the same time, the choice of basis set typically did not change the linewidth by more than 1-2%.
In spite of the increased accuracy of the electronic model, the calculated total FWHM values are about 40% greater than the experimental ones. The largest remaining error in the calculations is not due to the accuracy of the bound-electron wavefunctions, but rather due to the neglect of the details of the positron wavefunction and positron-electron interactions [13,14]. As expected, this will require going beyond the LEPWP approximation used here. The numerical scheme employed in the simulations also contains a small error due to the electron density loss in the core shells. Apart from the positron wavefunction, an important message from the shell-dependent results of the present study is that the Doppler-shift spectral widths ε may strongly depend on the shape of particular valence orbitals in the atoms and especially in the molecules, and this warrants further study.