Gaussian basis functions for highly oscillatory scattering wavefunctions

We have applied a basis set of distributed Gaussian functions within the S-matrix version of the Kohn variational method to scattering problems involving deep potential energy wells. The Gaussian positions and widths are tailored to the potential using the procedure of Bačić and Light (1986 J. Chem. Phys. 85 4594) which has previously been applied to bound-state problems. The placement procedure is shown to be very efficient and gives scattering wavefunctions and observables in agreement with direct numerical solutions. We demonstrate the basis function placement method with applications to hydrogen atom–hydrogen atom scattering and antihydrogen atom–hydrogen atom scattering.


Introduction
Variational approaches to quantum scattering have a long history and many versions have been devised. Examples include those of Kohn [1], Newton [2] and Schwinger [3]. Variational methods have been used for a diverse range of scattering problems including ones with deep potential energy wells and/or high collision energies, for example H+H 2 reactive scattering [4], electron-molecule elastic and inelastic scattering [5][6][7], positron-H 2 scattering [8], antihydrogen-H 2 scattering [9] and nuclear scattering [10]. The main idea of these methods is to expand the scattering wavefunction using a linear combination of basis functions. By varying the parameters of the basis the correct form of the scattering wavefunction can be obtained. The successful application of a variational method to a given quantum scattering problem thus requires the choice of a suitable basis set.
A particularly simple variational method was developed by Miller et al: the S-matrix Kohn variational principle (SKVP) [11,12]. This makes use of complex boundary conditions of the scattering wavefunction and has been shown to be essentially free of so called 'Kohn anomalies' which can arise for real boundary conditions [13,14]. The SKVP has been widely applied both to model systems [15,16] and reactions such as H+H 2 [4], H+HD [17], F+H 2 [18] and He+H 2 + [19]. One-dimensional scattering problems can of course be solved numerically using a suitable propagation method for the wavefunction and extracting the phase shift or another equivalent parameter to calculate the scattering cross section [20]. In higher dimensions however such as reactive atommolecule collisions, other approaches are required such as coupled-channel methods which typically use hyperspherical coordinates [21,22]. A great advantage of the SKVP is the simplicity by which it can be generalised to inelastic and reactive collisions, the only difference being the dimensions and complexity of the integrals which are required compared with one-dimensional scattering. The method also makes no assumptions about the form of the scattering wavefunction and thus is suitable for antimatter-matter collisions for which the potential energy surfaces (PES) typically have no classically forbidden regions.
Scattering problems involving deep potential wells result in highly oscillatory wavefunctions. Applying variational methods to these problems can lead to difficulties if conventional translational basis sets such as Slater or equally spaced Gaussians are used. If the potential also changes character quickly over a small range (such as going from attractive to repulsive) this further complicates the choice of basis set and can lead to inefficiencies.
In this work we apply a distributed Gaussian basis set within the SKVP to problems involving highly oscillatory scattering wavefunctions. The Gaussian basis set is tailored to the potential using a modified version of the methods of Light et al [23,24]. This will be shown to give very efficient and general basis sets for scattering problems which are even suitable for potentials where conventional basis sets fail. In the next section a brief description of the SKVP will be given and the distributed Gaussian placement method described. In section 3 this method is applied to H-H scattering and compared to direct numerical solutions. In section 4 the method is applied to antihydrogen (H)-H scattering. We discuss the method and present our conclusions in section 5.

SKVP and distributed Gaussian functions
We illustrate the Gaussian basis function placement method with application to one-dimensional elastic scattering. The Schrödinger equation for a particle of mass μ scattered by a central potential V(R) (or equivalently two particles in the centre of mass frame) is where E is the energy of the scattered particle, conserved for an elastic collision. In order to find (axially symmetric) solutions it is convenient to expand the scattering wavefunction as where the A l are expansion coefficients and the P l are Legendre polynomials. This is the partial wave expansion of the total scattering wavefunction [21]. Substituting this expansion into equation (2.1) gives a radial Schrödinger equation for each value of the orbital angular momentum quantum number l   m m y y - Since terms in the full scattering wavefunction Y µ , it follows that ψ l (0)=0 so that the total scattering wavefunction remains finite. If the potential energy, V(R), decays faster than 1/ R then the asymptotic solution of equation (2.3) is where the wavevector  m = k E 2 2 and δ l is the phase shift [21]. The asymptotic solution for ψ l can also be taken to have the complex form 2i l is the scattering or S-matrix, for elastic scattering a 1×1 matrix. This form of boundary condition is used in the S-matrix Kohn variational method. The SKVP was derived by Zhang et al including discussions of inelastic and reactive collisions [12]. Zhang and Miller also give a detailed discussion of applying the SKVP to atom-diatom inelastic and reactive scattering [17]. We direct the reader to these papers for full details of the SKVP.
Here we give only the working equations for one-dimensional, s-wave (l = 0) elastic scattering.
The trial scattering wavefunction is expanded as where the c n are variational parameters. In this expansion u 0 (R) is the incoming wave The function f (R) is a cut-off function for which f (0)=0 and = It is included so that u 0 (0)=0 as required from equation (2.2). In applications we have used where r 0 controls where the function cuts off and p determines how quickly this occurs. In the applications described in the next sections we have used r 0 =8.0 a 0 and p=8.0.
where, following the convention of Zhang et al, the bra is not complex conjugated [12]. The 'M' matrices are given as follows. M 0,0 is the matrix element between incoming waves Similarly M 1,0 is the matrix element between incoming and outgoing waves M 0 is a vector with elements between the N−1 basis functions with n2 and u 0 (R) with each element given by M is a square matrix with elements between each of the N−1 basis functions given by For all of these matrix elements H is given by the terms in square brackets in equation (2.3) with l=0 for s-wave scattering and E is the scattering energy of interest.
Unlike bound-state problems, where an upper bound to the energy can be obtained and systematically lowered by varying the basis set, variational approaches to scattering determine a certain quantity, here the S-matrix, such that it is stationary with respect to variations in the trial wavefunction [25]. Convergence of the S-matrix is checked by adding basis functions and making small changes to nonlinear parameters [12,13] and by checking that the S-matrix is unitary [17].
In applications of the SKVP to reactive scattering, equally spaced distributed Gaussian functions [23] have frequently been used for the translation coordinate [4,[17][18][19]. For the atom-diatomic reactions considered, the PES contain only mildly attractive regions and repulsive regions where atoms are close together. The scattering wavefunctions at low energies are thus smoothly oscillatory, decaying when the potential becomes repulsive. To converge the calculations with respect to the translational basis, the number of Gaussians within a given range is increased and the Gaussians are made narrower.
A problem with equally spaced distributed Gaussian basis sets occurs if the potential has deep wells in a certain region while being repulsive or mildly attractive in other regions. In this case the Gaussian basis must be able to reproduce the most oscillatory part of the scattering wavefunction in the potential well. However, this will be inefficient as the Gaussians will not need to be as narrow in other regions. As will be shown in section 4, for very attractive potentials (such as¯-H H) equally spaced Gaussian functions are completely unsuitable.
Instead of using equally spaced Gaussian functions as a basis set we tailor the positions and widths to the potential. This approach was used by Hamilton and Light in their original distributed Gaussian paper [23]. In applications we have used a placement procedure based on that of Bačić and Light [24]. To the best of our knowledge this approach has not yet been applied to scattering problems.
The placement method is as follows. An initial Gaussian function, is placed at R min . The choice of R min depends on the potential. For fully attractive or mildly repulsive potentials it is set at a small value close to the origin. For highly repulsive potentials it should be set sufficiently far into the classically forbidden region so that the scattering wavefunction is negligible. The exponent, α i , of the Gaussian is related to the de Broglie wavelength at that point via where C R is a parameter to be chosen and E is the scattering energy. This Gaussian is then normalised to unity. The next Gaussian, G j (R), is placed at some distance R j with its exponent calculated in the same way and normalised. The overlap of these two Gaussians is then calculated by the usual integral If the value of S is equal to S R (a chosen parameter) to within a set tolerance, S tol , then the Gaussian at R j is accepted and placed. If the overlap is too low or high then the Gaussian at R j is moved toward or away from R i respectively, until an acceptable overlap is obtained. This procedure is continued until a Gaussian is placed above R max . The value of R max is chosen to be effectively out of the range of the potential. This method of placing the functions is thus tailored to the potential. In regions where the potential is very attractive, narrow Gaussians are placed whereas for repulsive or mildly attractive regions, wide functions are used. This is an efficient placement scheme since more basis functions are placed where they are needed. The placement procedure also ensures nearly equally overlapping basis functions within the entire range. Thus six parameters are used to generate the basis: R min , R max , C R , α min , S R , S tol . R min should be small enough so that the scattering wavefunction is negligible there whilst R max should be large enough that the potential is essentially negligible. R min /R max should not be decreased/increased so much that the number of basis functions is increased unnecessarily. C R and S R should both be increased to increase the basis set size until convergence is achieved: increasing C R gives narrower Gaussians whilst increasing S R puts the Gaussians closer together. S tol is chosen to be small enough that overlaps close to the desired S R are achieved but not so small that the placement process takes a significant time. Convergence is not sensitive to the value of α min provided it is large enough to ensure the least oscillatory part of the wavefunction (within the range of the potential) can be reproduced. This value can be estimated using equation (2.18). For low total scattering energies the value of α min determined in this way might give rise to Gaussians which are very wide relative to most of the others (and the range of the potential) and so increase the range of quadratures required. In this case α min can be arbitrarily increased, keeping in mind that this would again increase the number of basis functions. An advantage of a bigger α min is that it may not need adjusting for each of a series of different scattering energies. For ¹ l 0 scattering calculations, the angular momentum gives rise to an effective repulsive potential which should be included in V(R) before the basis is generated.
An important modification of simple Gaussian functions is required for scattering problems. At R=0, ψ(R)=0 to ensure that the full wavefunction, Ψ=ψ(R)/R, is finite (equation (2.2)). For repulsive potentials ψ(R) is 0 before the origin and so simple Gaussians can be used. For attractive potentials however, ψ(R) is only 0 at the origin with finite amplitude even at close range. To ensure that the condition ψ (0)=0 can be met, the Gaussians are multiplied by a smooth cut-off function which goes to 0 as R goes to 0. In applications we have used basis functions of the form for n2, where A determines how quickly the cut-off function (1−e −AR ) goes to 0 and G n (R) is given by equation (2.17). The cut-off function goes to 1 for sufficiently large R and so normal Gaussian functions are recovered. The value of A is not crucial and in the applications here we set A=1a 0 1 . For brevity these basis functions will all be referred to simply as 'Gaussians' henceforth. The overlap integrals using this form of basis function are still analytical.
To apply these basis functions within the SKVP we have used a combination of analytical and numerical integration: Laguerre quadrature was used for integrals for the M 0,0 and M 1,0 matrices while Legendre quadrature was used for M 0 and M.
Atomic units will be used from this point onwards unless stated otherwise.

Hydrogen atom-hydrogen atom scattering
The first system we apply the distributed Gaussian placement method to is ground state, l=0 H+H elastic scattering treated within the Born-Oppenheimer (BO) approximation using the SKVP. For the Hamiltonian in equation (2.3) we set the reduced mass, μ, to half the proton mass (μ= m p /2=918.0736). The ab initio energies of Kołos et al [26] were used to fit the 1 Σ + g potential energy curve (PEC). The PEC was fit to an analytical form as follows. The ab initio electronic energies of Kołos et al were fit using an analytical form consisting of a sum of sigmoid functions that is, a one-dimensional neural network [27] where d 0 , d p , w p and b p were determined by least-squares fitting. A nine-node network with 28 free parameters was used to fit all 55 ab initio energies between R=0.2−12 a 0 . At long range the interaction energy was described using the dispersion energy formula  [28]. The two ranges are smoothly joined using a switching function. Following Cvitaš et al [29] the switching function is with a=3.0a 0 1 and s=11.0 a 0 . For the combined analytical form a root mean square deviation (rms) of 0.26 mE h was achieved with respect to the ab initio data. For R<0.2 a 0 we set V(R)=V(0.2) as for low energies the scattering wavefunction is negligible here. Thus the total potential is given by H h is the energy of a hydrogen atom. The H-H PEC is shown in figure 1.
The following parameters were required to generate a Gaussian basis set for a converged scattering calculation. R min was set to 0.1 a 0 and R max to 15.0 a 0 . A value of 1.0a 0 2 was used for α min . C R was set to 0.2 and S R to 0.9 with S tol at 0.05. This large value of S R is required due to the H-H PEC becoming steeply repulsive below 1 a 0 and so if a smaller value is used the Gaussians are placed too far apart. For efficiency the value of S R can be changed within different ranges but for simplicity it is kept constant here. Using these parameters generates 73 Gaussians. The narrowest function is placed at R=1.398 a 0 with a » with C=0.75. The scattering calculation is converged by increasing the number of functions. To converge the calculation, 150 equally spaced Gaussians were required with α n ≈57a 0 2 . This is similar to the value of the exponent of the narrowest Gaussian placed using equation (2.18). The inefficiency of equally spaced Gaussians is thus highlighted: sufficiently many functions must be added until the exponents of each are large enough to be able to reproduce the most oscillatory part of the scattering wavefunction.
As an independent check on the accuracy of the SKVP calculations we also solved the radial Schrödinger equation, equation (2.3), numerically by directly integrating the scattering wavefunction using the Numerov method [30]. The wavefunction was propagated from R=0.1 a 0 (where ψ(R) was set to 0) outward using a step size of 1×10 −3 a 0 . The scattering phase shift, δ 0 , is found by matching the propagated solution to the asymptotic form, equation (2.4) [20].
The H-H radial scattering wavefunctions, ψ 0 , generated using both the Numerov method and the SKVP with the tailored Gaussian basis set are shown in figure 2 for E=1×10 −10 E h . For the SKVP the complex function is multiplied by  [32][33][34] and potential [35] used greatly affecting its value. Since the emphasis here is on the basis set we did not attempt to improve our potential fit any further.

Antihydrogen atom-hydrogen atom scattering
The¯-H H system provides a demanding test for variational approaches to scattering due to the strongly attractive nature of the PEC. The potential is attractive at all distances and as  R 0 it becomes Coulombic when the antiproton and proton are close together. This results in highly oscillatory scattering wavefunctions. Previous work using Kohn methods for this system did not converge the calculation of the elastic cross section using a basis set for the radial coordinate [36]. For this reason a numerically propagated BO radial scattering wavefunction was incorporated in a subsequent reactive calculation [37]. Finding suitable scattering basis functions for this and related systems is thus an important problem.
Elastic¯-H H scattering will be treated here within the BO approximation. For this system the approximation is unreliable at internuclear distances below around 0.744 a 0 [28]. At and below this critical distance, R c , the leptons (electron and positron) can dissociate from the nuclei as positronium (Ps, a bound state of an electron and positron), rendering the BO approximation invalid. This makes carrying out physically meaningful scattering calculations challenging. Despite this the¯-H H system treated within the BO approximation is an ideal model potential for developing methods since different groups using different approaches have reported similar results [36,38]. The reduced mass of the system, μ, was again set to m p /2 here (and also to be consistent with the literature cited).
An analytical PEC of the¯-H H system was fit in a similar way to that described above for H-H. Ab initio leptonic energies computed by Strasburger [28] were fit using a onedimensional neural network. A six node network with 19 free parameters was used to fit all 45 ab initio energies between R=0.744 and 20 a 0 to give V lep (R). At longer range the same dispersion energy formula, equation (3.2), was used with the same coefficients as those for H-H. The neural network fit and dispersion energy functions were smoothly combined using equation (3.3) with = a a 3.0 0 1 and s=7.0 a 0 . Using this form a rms error of 2.90 mE h was achieved with respect to the ab initio data. For RR c the leptonic energy was set equal to that of the ground state of positronium, −0.25 E h and V(R) obtained by adding this to the Coulomb interaction of the nuclei, −1/R. Thus the total potential is given by The¯-H H PEC is shown in figure 3. Despite the very attractive PEC for¯-H H, generating a basis set for converged scattering calculations using the SKVP was still possible with the placement procedure described above. The following parameters were used. R min was set to 1×10 −5 a 0 and R max to 12 a 0 . Such a small value of R min is required as the scattering wavefunction, ψ 0 (R) is finite at all ranges with only ψ 0 (0)=0 (as required from   [31] 0.5425 equation (2.2)). C R was set to 0.3 with S R =0.7 and S tol = 0.05. α min was again set to 1.0a 0 2 . Using these parameters generates 80 Gaussians. The narrowest function is placed at R=R min with a »´a 5.5 10 We also attempted the SKVP calculation using a basis set of equally spaced Gaussians but convergence was not reached despite using hundreds of functions-see table 2. This is not surprising since, as discussed in the previous section, the exponent of the equally spaced Gaussians must be around the same value as the narrowest function generated using the placement method. Using this exponent (a »´a 5.5 10 2 ) would required the Gaussians to be spaced by around 10 −4 a 0 . For the range 0−12 a 0 required this would mean using about 10 5 functions. Clearly equally spaced Gaussians are totally inadequate for this system.
We again numerically solved the radial Schrödinger equation using the Numerov method to check our SKVP results. The wavefunction was propagated from R=1×10 −7 a 0 outward using a step size of 1×10 −4 a 0 .
The¯-H H radial scattering wavefunctions, ψ 0 , generated using the Numerov method and the SKVP using the tailored Gaussian basis set are shown in figure 4 for E=1×10 −8 E h . The tailored Gaussian basis is capable of reproducing the scattering wavefunction accurately despite its highly oscillatory nature. The scattering lengths given by the Numerov and SKVP methods are shown in table 2 along with those of Jonsell et al [38] and Armour and Chamberlain [36]. The Numerov and Kohn values calculated here are in good agreement with each other and with the literature values which were calculated using different leptonic ab initio data.
For¯-H H scattering a further check on the accuracy of the SKVP elastic scattering wavefunction can be carried out by using it to compute annihilation cross sections. For this system both hadronic (antiproton-proton) and leptonic (positron-electron) annihilation can occur. Both annihilation cross sections were computed using the 'delta' potential method following Jonsell et al [38,39]. This approach uses the elastic scattering wavefunction to compute annihilation cross sections and is classed as a first order perturbation type method, ignoring any changes in the elastic scattering wavefunction caused by the strong force.
For the hadrons the annihilation cross section is computed using the value of the scattering wavefunction (Ψ(R)=ψ 0 (R)/R) at R=0 where the antiproton and proton's location coincides. Details of this calculation are given by Jonsell et al [38], but briefly, the antiproton-proton annhilation cross section is whereĀ pp is the rate constant for the annihilation. The value of A pp is taken to be 1.7×10 −7 E h a 0 3 [38]. At low energies the annihilation cross section is related to the scattering energy by where C a is a constant. Using the elastic scattering wavefunction generated from the SKVP we compute a value of 0.11 a E 0 2 h 1 2 for C a . This is in reasonable agreement with Jonsell et alʼs value of 0.14 [38] given that the latter authors used a different potential.
Leptonic annihilation is calculated in a similar way but must take into account the finite probability of annihilation over the whole range of R. Froelich et al have given a detailed discussion of this calculation [39], but briefly the positronelectron annhilation cross section is given by   where

+ -
A e e is the rate constant for the annihilation and P(R) is the positron-electron coalescence probability at separation R.

+ -
A e e is taken to be 4.86×10 −6 E h a 0 3 for singlet, two-photon, annihilation collisions and 4.28×10 −9 E h a 0 3 for triplet, three-photon, annihilation collisions [39]. We used Strasburger's positron-electron coalescence probabilities for  R R c [40]. Below the critical distance we use Strasburger's suggestion of simply setting the coalescence probability to p From the annihilation cross sections calculated we can conclude that the Gaussian basis set implemented within the SKVP gives accurate scattering wavefunctions which can be used to compute scattering observables. The calculation of the hadronic annihilation cross section shows that the basis set gives an accurate representation of the wavefunction close to and at R=0 while the leptonic annihilation cross section requires accurate values of the wavefunction over a wide range of R values.

Discussion and conclusions
We have applied a modified version of the Gaussian placement procedure of Bačić and Light to generate basis sets for use in variational scattering calculations, specifically using the SKVP. We have shown that the basis sets generated using this method are especially suitable for potentials with deep wells resulting in highly oscillatory wavefunctions. Another application could be to high-energy collisions where the quantity [E−V(R)] is also large and so the scattering wavefunction is again highly oscillatory. The basis sets are also efficient, placing many narrow functions in regions of low potential energy while placing significantly fewer, wider functions in repulsive or mildly attractive regions. Such potentials are difficult or impossible to treat using equally spaced Gaussian basis sets which have commonly been applied to variational approaches to atom-diatom reactive collisions. The wavefunctions generated using the placement method are accurate and can be used to calculate scattering observables.
There are other approaches which can be used to deal with highly oscillatory one-dimensional wavefunctions. Basis functions such as high-energy harmonic oscillator functions or hydrogen atom-like radial wavefunctions with large N quantum number are intrinsically oscillatory. Such basis sets would be able to cope with the oscillatory nature of scattering wavefunctions in attractive potentials but are awkward to use and difficult to generalise to multidimensional systems. The use of such basis sets would also require tests to assess which functions were suitable for a given problem. Another approach would be to use the Gaussian expansion method (GEM) of Hiyama et al [41]. This method uses a basis set consisting of Gaussians centred at R=0 with different decay rates. A modification for oscillatory wavefunctions involves multiplying these functions by sine or cosine terms The GEM has been widely applied to multidimensional atomic and nuclear problems [41] and is clearly a useful method. A deficiency of the method however is the requirement of choosing parameters of the basis. The authors recommend using a geometric series which fixes the exponents of the Gaussians and frequency of the sine and cosine functions. The series still requires a choice of three parameters which need to be optimised for each system apparently in a non-systematic way.
The discussion in the previous paragraph should be contrasted with the Gaussian basis set placement method used here. The method of placing the Gaussians is straightforward. The range and character of the potential determines the range of where the Gaussians are placed. The basis set can be systematically improved by increasing C R and S R until scattering calculation convergence is obtained. Gaussian functions are also simple to manipulate in order to find derivatives or integrals [23].
We stress that one-dimensional systems have been used here to demonstrate the use of the placement method and ideas but that higher-dimensional problems are of course the main interest. Multidimensional Gaussian basis sets have been successfully applied for bound-state problems [23,[42][43][44]. For bound-state problems the basis set is only required over a small range for a given energy. In a future publication we will expand the Gaussian basis placement procedure to rigid-rotor scattering using a multidimensional Gaussian basis set.