Hydrogen molecule-antihydrogen atom potential energy surface and scattering calculations

We have calculated ground state interaction energies for an antihydrogen atom and a hydrogen molecule within the Born–Oppenheimer approximation. Leptonic energies were calculated using a large basis set of explicitly correlated Gaussian functions. Energies were calculated at over 2800 geometries including different proton–proton distances. The energies have been fit to functional forms using a neural network for the short-range interaction which is combined with asymptotic formulas at long range. A two-dimensional rigid rotor and a three-dimensional atom–molecule potential energy surface (PES) have been determined. Rigid-rotor scattering calculations on these surfaces have been carried out using the S-matrix Kohn variational method with a two-dimensional Gaussian basis set. We have calculated cross sections for elastic, rotationally inelastic and annihilation collisions on the two-dimensional PES. This includes the first calculation of leptonic annihilation for this system.


Introduction
Experiments to produce, confine and characterise antihydrogen atoms (H) [1][2][3][4][5][6] have inspired theoretical work on the interactions ofH with normal matter atoms and molecules [7]. These interactions lead to the destruction ofH and so are important from an experimental perspective while also providing a stringent test for molecular dynamics theories.
-H H 2 is the simplest system for the interaction of an antiatom with a neutral molecule and thus is an important benchmark. It is also of interest as H 2 is present within the antimatter traps and is expected to contribute to antihydrogen losses. The latter may occur due toH gaining kinetic energy in elastic collisions with H 2 , particle rearrangement processes (see below), as well as in-flight proton-antiproton and electron-positron annihilation. The relative and absolute rates (or probabilities or cross sections) for these different collision outcomes vary by many orders of magnitude as a function of collision energy. A detailed knowledge of this system is therefore of theoretical and practical interest. Despite this, there has been limited theoretical work carried out on the system, in particular at the low collision energies (corresponding to temperatures from room temperature down to near absolute zero) relevant to the experiments cited above, with only four papers published to date to our knowledge. One of us has calculated a partial potential energy surface (PES) using an explicitly correlated Gaussian (ECG) basis set within the Born-Oppenheimer (BO) approximation [8]. The potential energies were calculated for linear and perpendicular geometries. For some systems containing both matter and antimatter nuclei, the leptons are unbound for certain configurations of the nuclei: for example in H-H when the internuclear distance is below a critical value, R c [9]. This situation does not arise in the present case: for both one proton and the antiproton very close or all three heavy particles nearly coincident, the 'united atom limit' corresponds to {proton, electron, electron, positron} which has one bound state [10], so the three leptons remain bound in the BO approximation.
Only limited scattering calculations have been carried out for the¯-H H 2 system. Gregory and Armour carried out quantum mechanical, Kohn variational, calculations for elastic scattering at very low energies (10 −4 E h , were E h represents Hartree, the atomic unit of energy; in temperature units 30 K) using a quantum Monte Carlo PES and treating the H 2 as a rigid rotor [11]. Prolate spheroidal coordinates were used to describe the position of the antiproton relative to the H 2 and the trial variational wavefunction included terms composed of products of powers and exponentials of those coordinates. An estimate of the in-flight hadronic annihilation cross section was also obtained using a delta potential approach. The latter is very sensitive to the quality of the wavefunction when the antiproton is close to one of the protons and so specially-designed terms were added to the basis set for such configurations. Further quantum scattering calculations were carried out on this system, also using the rigid-rotor approximation, by Sultanov et al [12]. They used a coupled-channel approach with Jacobi coordinates and an approximate PES constructed from the¯-H H potential energy curve (PEC) [9]. In this approach, the radial wavefunction is determined by numerical propagation rather than with a basis set. This gave a low-energy elastic cross section in rough agreement with Gregory and Armour's results. Cross sections and rate constants for various rotational transitions of para-H 2 were also calculated. Neither Gregory and Armour nor Sultanov et al attempted reactive (rearrangement) scattering calculations. The latter has been carried out by Cohen for both¯-H H 2 and¯-+ H H 2 using the quasiclassical fermion molecular dynamics (FMD) method at energies above 0.1 E h [13]. Both systems were found to be highly reactive with cross sections for protonium (Pn) formation processes largest at the lowest energies considered. For¯ At the lowest energies considered the reaction forming Pn, Ps (positronium) and H was found to be most probable followed by the formation of the metastable PnPs complex, although this may be an artefact of the FMD method.
The focus of the present work is the interaction of an antihydrogen atom with a hydrogen molecule. We have carried out ab initio calculations of leptonic energies at a large number of geometries. Two and three-dimensional PESs have been constructed by fitting analytical functions to these energies. These surfaces (as opposed to one-dimensional cuts) are the first published for the¯-H H 2 system. We have also used these PESs to carry out quantum scattering calculations using the S-matrix Kohn variational principle with a twodimensional Gaussian basis set. Scattering calculations have been carried out for low-energy rigid-rotor elastic, inelastic and annihilation collisions, but no account has been taken of reactive processes. A full quantum-mechanical treatment of the latter would be significantly more challenging.
In section 2 below we give details of the atom-diatom scattering method, the Gaussian basis set used and the treatments of leptonic and hadronic in-flight annihilations. In section 3 we report and discuss the two-and threedimensional PESs and the results of elastic and inelastic scattering calculations including annihilation. Concluding remarks are given in section 4. The scattering calculations reported in detail here are restricted to the rigid-rotor model, that is two-dimensional heavy-particle dynamics. For the latter, the two-dimensional PES has been used and, to test the sensitivity of the dynamics calculations to changes in the PES, also the three-dimensional PES with fixed H 2 bond length. The three-dimensional PES should also be of use for three-dimensional heavy-particle dynamics, including preliminary treatments allowing the H-H distance to vary, and for¯-H H 2 reactive scattering calculations. These latter applications are briefly discussed in sections 3.5 and 4.

Atom-molecule scattering and S-matrix Kohn variational method
The calculations of elastic and inelastic scattering of¯-H H 2 use Jacobi coordinates: R, the distance from the antiproton to the centre of mass of the two protons, r the distance between the two protons and θ the angle between the vectors corresponding to R and r. The present calculations are restricted to total angular momentum J=0 (including the angular momenta of the rotations of R and r but excluding the spins of the nuclei and leptons). The Hamiltonian in this case is [14]  where μ is the reduced mass of the atom-molecule system, μ r is the reduced mass of the molecule and V(R, r, θ) is the potential energy. For the rigid-rotor calculations, r was fixed at the equilibrium value for H 2 , 1.4011 a 0 , and the term involving the second derivative of r was omitted. The scattering calculations use the S-matrix Kohn variational principle (SKVP). This method was derived by Zhang et al [15]. Zhang and Miller also gave a detailed discussion of applying the SKVP to atom-diatom inelastic and reactive scattering [16]. In appendix A we give the working equations for applying this method to J=0 rigid-rotor elastic and inelastic scattering with the Hamiltonian of equation (2.1).
The conventional approach to solving the timeindependent Schrödinger equation with this Hamiltonian is to expand the wavefunction for a given total angular momentum using a basis set consisting of a radial part for R and coupled spherical harmonic functions for the rotation of R and r. The latter reduce to Legendre polynomials, P j , in the case of J=0. In a preliminary attempt, not reported here, rigid-rotor S-matrix Kohn variational elastic scattering calculations on the¯-H H 2 system used a basis set consisting of a direct product of distributed Gaussians along the R coordinate (placed as described below) and rotational functions for the angular coordinate θ. However rotational functions up to and including j=100 were used in the calculation with no sign of convergence. This failure is attributed to the highly anisotropic nature of the PES due to the Coulombic singularities when the antiproton approaches the protons.
For the results reported below, a basis set consisting of two-dimensional Gaussian functions (of R and θ) was used. In a previous paper [17], we have shown that distributed Gaussian functions with widths and placements tailored to the potential are effective for highly oscillatory scattering wavefunctions (such as in the¯-H H system) which have very different characters in different regions. The same problem is present for¯-H H 2 where the potential changes dramatically with the geometry of the particles. For this reason the Gaussian basis set placement method of our earlier paper [17] was extended for rigid-rotor scattering using two-dimensional Gaussian functions.
Multidimensional Gaussian functions have been used in variational calculations for vibrational bound states of atomic complexes [18], atom-molecule complexes [19,20] and the water molecule [18]. In these applications Gaussians were only required in regions of low potential energy, that is, where the wavefunctions have significant amplitude.
It appears that multidimensional Gaussian functions have not been used for scattering calculations, although they were suggested as a possibility by Zhang and Miller [16]. Thē -H H 2 rigid-rotor PES is totally attractive with no classically forbidden regions. In this case placing the Gaussian functions in a suitable way is not as straightforward as for bound-state problems. Nevertheless, the fact that multidimensional Gaussian functions can be tailored to the potential makes them attractive for use if a suitable placement scheme can be devised. Basis functions in R and θ were used here of the form with the radial cut-off function, ( ) -e 1 AR , included to ensure the basis functions are finite for R=0 where A determines how quickly the cut-off function goes to 0 [17].
This form of 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=2a 0 1 . For brevity these basis functions will all be referred to simply as 'Gaussians' henceforth.
An angle, θ R , is chosen to generate a radial Gaussian basis. The angle should be for the most attractive angular cut through the PES which for the¯-H H 2 system is θ R =0 or π.
In practice a small value of θ R is chosen to avoid the Coulomb singularity. A radial Gaussian basis is then generated along this angle using the same placement procedure as described in our earlier paper [17] based on the method of Bačić and Light [21]. An initial Gaussian function, , is placed at R min , a small value close to R=0. The exponent, a i R , 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. The next Gaussian, G j (R), is placed at some distance R j with its exponent calculated in the same way. The overlap, S of these two normalised Gaussians is then calculated and if the value of S is equal to S R (a chosen parameter) to within a set tolerance, S Rtol , 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 the correct overlap is obtained. This procedure is continued until a Gaussian is placed above R max (out of the range of the potential). The basis size can be increased by increasing the C R parameter which makes the Gaussians narrower or by increasing S R which puts the Gaussians closer together. The value of a R min should be as small as possible and will depend on E. Angular Gaussians, ( ) , are then generated between q min and θ max at each radial Gaussian point. These angular Gaussians are placed in the same way as the radial functions but with a q sin term from the volume element included in overlap integrals which have the range 0 to π. Thirteen parameters are required for generating the basis: C R , a R min , S R , S Rtol , R min , R max , C θ , a q min , S θ , q S tol , θ min , θ max and θ R . The parameters for the angular Gaussians in θ are analogous to those for R.
This placement procedure ensures that basis functions are placed over the whole range of space leaving no 'gaps'. A deficiency of the method is the efficiency for the radial range since the radial basis is generated for the most attractive angular cut through the PES. This procedure also makes no use of the symmetry of the¯-H H 2 system. Gaussians in regions where the potential is very low have only a small range. Thus matrix elements between widely separated Gaussians are negligible. This feature is important for making multidimensional Gaussian basis sets efficient so that even for very large basis sets, only a small percentage of the integrals need to be evaluated. Here, if the overlap between two Gaussians was less than =´-S 1 10 tol 12 then the M matrix element (see appendix A) was set to zero.

Annihilation
Both hadronic and leptonic in-flight annihilation cross sections are calculated using the 'delta potential' method. This assumes that annihilation occurs when a particle and corresponding antiparticle's position coincide. The method also assumes that annihilation has no effect on the elastic scattering wavefunction and treats the process as a first order perturbation.
The delta potential treatment of hadronic annihilation in mixed matter/antimatter systems has been applied to¯-H H [22],¯-H He [23] and¯-H H 2 [11] scattering. For the¯-H H system the simple delta potential approach works well and gives hadronic annihilation cross sections in good agreement with a more sophisticated scattering length approach [23] and complex potentials [24]. The elastic scattering cross section is only corrected by around 10% when hadronic annihilation is taken into account. For the¯-H He system the delta potential method underestimates the hadronic annihilation cross section compared with the scattering length approach [23] and optical potentials [25] giving around twice as large a value. Accounting for hadronic annihilation in¯-H He scattering also drastically changes the elastic cross section. For the¯-H H 2 system a more sophisticated treatment of hadronic annihilation is also expected to affect the elastic cross section [11]. Despite this the delta potential treatment of hadronic annihilation is the simplest approach and easy to generalise to systems with more than two hadrons. It should give at least the order of magnitude of the hadronic annihilation cross section.
Here¯-H H 2 hadronic in-flight annihilation cross sections were calculated from the Kohn rigid-rotor elastic scattering wavefunctions (determined using the twodimensional PES) as 3 from the width of the protonium 1s state [22]. See appendix B for our derivation of equation (2.4) following Jonsell et al [22].
Leptonic annihilation cross sections were calculated using an extension of the method described by Froelich et al [26] (see appendix B) as The positron-electron coalescence probability density, P(R, θ), is a function of R and θ. One of us has previously calculated the values of P(R, 0) and ( ) p P R, 2 using an ECG basis set [8]. To calculate the leptonic annihilation rate the averages of the linear and perpendicular P(R, θ) values are used. The values of P(R) were least-squares fit using a sixnode neural network (NN) giving a standard deviation of 3.1×10 −6 . Figure 1 shows the calculated values of P(R, 0) and ( ) p P R, 2 , the average value and the fitted function.
Although the values of P(R, 0) and ( ) p P R, 2 differ by tens of percent at a given R, both decrease at essentially the same rate and are negligible by R=10 a 0 . Using the average of the values for θ=0 and p 2 should give the correct order of magnitude of the leptonic annihilation rate. The fitted function to P(R) has some deficiencies at small R values but the form of P(R) at these distances is not important as small values of R do not significantly contribute to the integral in equation (2.5) as discussed in the next section.

Ab initio calculations
The leptonic contribution, E lep , to the PES has been computed in the BO (clamped nuclei) approximation using a similar method to that published by one of us for linear and perpendicular geometries [8]. See appendix C for details of the present calculations. The energies of 2805 unique geometries of the¯-H H 2 system were calculated. The geometries are generated using prolate spheroidal coordinates for a range of p-p distances. These coordinates are defined as r 3 , m = + r r r 1 2 3 and n = r r r where r 3 is the p-p distance and r 1/2 are the two¯-p p distances. Fourteen r 3 values between 0.5 and 20 a 0 were used including the H 2 equilibrium bond length of 1.4011 a 0 . For each value of r 3 ten values of ν were used between 0 and 1. ν=0 and ν=1 correspond to perpendicular and linear geometries respectively. The values of μ varied for each r 3 and ν.

Analytical representation of PESs
The PESs have been constructed by combining analytical functions least-squares fit to the ab initio energies, E lep , with the nuclear Coulombic interaction energies and long-range terms as described in detail in appendix C. Very little work has been done on fitting analytical forms to the short-range potentials of mixed matter/antimatter systems. For this reason a NN approach [28,29] which is capable of fitting any shape of potential was used for the short-range interactions. For asymptotic geometries, functional forms based on long-range dispersion energy formulas were used.
3.2.1. Two-dimensional PES. The leptonic energies of 293 geometries with the p-p distance fixed at 1.4011 a 0 were fit using a two-dimensional NN. All geometries were equally weighted. The fit was carried out using 12 symmetrically unique nodes with 49 independent parameters. For the best fit a standard deviation of 0.036 E m h was obtained although achieving similar accuracy with different starting parameters is straightforward. The largest discrepancy between an ab initio energy and the fitted surface was 0.12 E m h . The complete functional form of the two-dimensional -H H 2 PES, including the nuclear Coulombic interaction energies and long-range terms (equation (4.25) in appendix C) gives a root mean square (rms) error to the ab initio energies of 0.031 mE h with a maximum deviation of 0.12 E m h . For the linear geometry the potential is attractive at all distances. For the perpendicular geometry there is a barrier at around 2.5 a 0 followed by a local minimum of ≈1 E m h depth at 3.8 a 0 , as shown in figure 2. The maximum of this barrier is still below the separated atom-molecule energy however. The reason for the barrier at perpendicular geometries and lack of one for linear geometries was discussed previously [8]. Cuts through the full PES are given in figure 3 showing how the barrier decreases as the angle is lowered. There are no classically forbidden regions in the two-dimensional PES.

3.2.2.
Three-dimensional PES. The leptonic energies of all 2805 ab initio data points were fit using a NN. All geometries were again equally weighted. The best 24-node fit (with 121 independent parameters) determined an analytical function ( ) V r r r , , NN 1 2 3 which gave a standard deviation of 0.44 E m h with a maximum deviation of 4.63 E m h . The worst agreement between the fit and data was found to occur for high leptonic energies where the antiproton is close to one or both of the protons. In this case the potential energy is dominated by the Coulombic interaction of the nuclei and so the disagreement is less important at these geometries.
After carrying out this fit, 79 additional ab initio data points became available. Rather than include these in the fit it was decided to use them as a check of the fitted NN. The additional data contained 77 points with all the nuclei within 6 a 0 of each other. The other two geometries were Pn-H type with thep close to one of the protons (0.5 and 0.6 a 0 ) but over 20 a 0 away from the other proton. The rms error for the NN fit to these points was 1.98 E m h . The largest discrepancy was for one of the Pn-H geometries at 16.01 E m h . Such a large difference between the fit and ab initio energy for this geometry is expected since so few of the initial data set contained this type of geometry. The next largest difference between the fit and data was then 2.25 E m h at a geometry with high leptonic energy.
The combined function, including the nuclear Coulombic interaction energies and long-range terms, equation (4.31) in appendix C, gives a rms error of 0.43 mE h for the 2805 data points used in the fit.
The barrier for perpendicular geometries is shown in figure 2. The three-dimensional fit has too high a barrier with the potential energy higher than the separated atom and molecule. The barrier is only too high by a few tenths of a millihartree however and since it decreases as the¯-H H 2 angle becomes more linear, it should not seriously affect scattering calculations. Figure 4 shows the barrier for different p-p distances. To fit all of the curves into one plot, H was subtracted from the total energy for each p-p distance shown. As the p-p distance is decreased the barrier increases. This is expected since for the¯-H He system there is a barrier of around 0.3 E m h at approximately the same distance [30]. The leptonic energy for H 2 approaches that of He as the p-p distance is shortened. It follows that for small p-p distances the barrier should be of a similar height. For small values of r the potential energy increases as R decreases from 7 to around 5.5 a 0 . This is an artefact of the fitted PES due to the limited accuracy of the NN fit. This is also responsible for the shape of the potential for all r between R=7-9 a 0 where the shortrange NN fit is switched to the long-range¯-H H 2 form. This  was not present for the two-dimensional fit (see figure 3) as the short-range NN function is more accurate. This is not a serious problem however as for these r values the total energy is higher than the isolated H 2 +H equilibrium energy by at least 36 E m h (for r=2.0 a 0 ) and thus will not affect lowenergy scattering calculations.
Another comparison of the two-and three-dimensional PES fits is shown in figure 5 which shows the differences between the surfaces in Jacobi coordinates. The largest differences coincide with the Coulombic singularities where, as discussed, the NN fit is least accurate. The region around the barrier for perpendicular geometries also shows larger differences. For small and large R distances the agreement between the surfaces is better.

Elastic and inelastic scattering
A suitable two-dimensional Gaussian basis set for the¯-H H 2 S-matrix Kohn variational scattering calculations was found by checking the convergence of low-energy elastic scattering using the two-dimensional PES. The parameters for the converged basis are given in table 1. The total energy of the scattering calculation is also required to place the basis set. However it was found that the basis generated with the lowest and highest energies considered for elastic scattering were almost identical and gave elastic cross sections at a given energy in agreement to four significant figures. This is due to the energies of interest being low and the potential being very attractive so that the optimum basis is essentially insensitive to the scattering energy.
The parameters of table 1 gave a basis consisting of 12617 two-dimensional Gaussian functions with 98% of these having R i < 2 a 0 . The basis function locations are shown in figures S1 and S2 in the supplementary information available online at stacks.iop.org/JPB/52/185201/mmedia. The density of functions is greatest around the two Coulomb singularities at R = r e /2 = 0.700 55 a 0 , θ = 0 and π. At and above R = 4 a 0 the exponents of the functions in both coordinates take their minimum values and the Gaussians are then equally spaced.
The SKVP elastic scattering calculations were carried out for total angular momentum J=0 using the Hamiltonian of equation (2.1). This should give a good approximation to the total elastic cross section at the low energies considered. The reduced mass of the atom and molecule, μ in equation With these parameters the incoming and outgoing functions are cut-off at around R=8 a 0 but the details of the cut-off function are not important. For example using r 0 =7.0 a 0 and p=7.0 changed the elastic scattering phase shift by less than 1%. This is also indirect evidence that the basis set is sufficiently large to give a converged calculation [15].
The elastic scattering phase shifts and cross sections for -H H 2 at selected energies are shown in table 2. A plot of the   phase shifts (modulo π adjusted) at different values of k is shown in figure S3 in the Supplementary Information where a fifth order polynomial has been fit to the data. The calculated phase shifts vary smoothly as k increases. A plot of the elastic scattering cross sections is shown in figure 6. At low energies the elastic cross section becomes constant as expected from threshold laws [31]. The cross section falls to zero when the phase shift passes through zero. The limiting elastic cross section for low energies is 1776 a 0 2 corresponding to a scattering length of 11.89 a 0 .
The results here can be compared with previous work. Gregory and Armour calculated the limiting elastic cross section to be 4801 a 0 2 giving a scattering length of 19.5 a 0 . They state that there was a 16% variation in the phase shift when nonlinear parameters in their basis set were changed. The limiting cross section calculated here is two or three times smaller. The difference in values could be due to the PES used. Gregory and Armour used ab initio energies calculated using a quantum Monte-Carlo method and are expected to be correct to 10-100 μE h . The PES used here was fit to ab initio energies expected to be accurate to a few μE h . Despite this the cross sections are in fair agreement, especially compared to¯-H He elastic scattering where small changes in the PEC resulted in very large differences in the elastic cross section [32]. The calculations of Sultanov et al gave a limiting elastic cross section of 9470 a 0 2 giving a scattering length of 27.5 a 0 . This is a larger disagreement with the calculations here. The PES used by Sultanov et al was constructed by adding two¯-H H PECs which is less accurate (by millihartrees) than the PESs used here and by Armour and Gregory.
Elastic rigid-rotor scattering was also carried out using the three-dimensional¯-H H 2 PES with r 3 =r fixed at the equilibrium value for H 2 in V(r 1 , r 2 , r 3 ). This allows a test of the sensitivity of the computed cross sections to changes in the PES by comparing the results with those of the twodimensional PES. The Gaussian basis set constructed for the two-dimensional surface was used with the same parameters for the Kohn calculation. The elastic scattering phase shifts for the two-and three-dimensional PESs are compared in figure S4 in the Supplementary Information where again the phase shifts were fitted to fifth order polynomials. The elastic scattering cross sections for the two-and three-dimensional PES are compared in figure 7. The limiting elastic cross section for the full three-dimensional surface is 3067 a 0 2 giving a scattering length of 15.62 a 0 . From figure 7 it can be seen that the behaviour of both cross sections with increasing energy is qualitatively similar. The agreement between the two-and three-dimensional surface is reasonable. The threedimensional PES gives a rms error of 0.29 mE h for the ab initio energies used to fit the two-dimensional PES which has a rms error of 0.031 mE h . Thus an order of magnitude increase in the error in the PES changes the computed limiting elastic cross section by less than a factor of two.
The sensitivity of the elastic cross sections for the threedimensional surface with respect to small changes in the PES was also investigated. It was found that small changes, such as changing the parameters of the switching function between the short-range and long-range forms, did not cause significant changes to the cross sections. For example changing s from 8 to 15 a 0 in equation (4.31) of appendix C so that the NN is used out to a longer distance gave a less than 1% change in the elastic cross section. This is in contrast to¯-H He scattering where small changes in how the BO potential energy curve was interpolated led to large changes in the elastic and annihilation cross sections [32].
For the higher energies considered rotationally inelastic scattering channels (that is, collisions where the molecule's rotational state can change) are open. As the H 2 molecule is homonuclear, only Δj=±2 transitions are possible. For total angular momentum, J=0, only even j states can be accessed from j=0, corresponding to para-hydrogen molecules. The j=2 rotational state has an energy of 1.66×10 −3 E h and is open for scattering energies at and above this. For two rotational states the S-matrix is a 2×2 matrix. The S-matrix for inelastic scattering also contains results for =  ¢ = j j 0 0 and =  ¢ = j j 2 2 elastic scattering.  The same Gaussian basis set and parameters of the elastic calculation were used for the inelastic calculations. The twodimensional PES was used. The cross sections are calculated using equation (4.11) in appendix A accounting for the degeneracy of the rotational state of the incoming scattering function. Table 3 shows the probabilities and cross sections for inelastic scattering. Note that the values of 1.01 and 1.02 for | | ¬ S 0 0 2 and | | ¬ S 2 2 2 respectively are greater than 1.00 due to numerical errors.
From table 3 it appears the rotational excitation and relaxation cross sections vary considerably with energy. For a scattering energy of 4.0 E m h the excitation and relaxation cross sections are very small. Calculations with smaller energy increments are needed to fully investigate these features.
Inelastic¯-H H 2 scattering was also considered by Sultanov et al using a coupled-channel approach [12]. Their results for rotational excitation from j=0 to j=2 and the corresponding relaxation process were reported for energies much closer to the j=2 threshold energy. For the energy which overlaps their range (1.7 E m h ) the s ¬ in2 0 cross section here is around twenty times lower than their result. As discussed in the previous sections they used an approximate PES and H 2 rotational basis which seems to be inadequate for thē -H H 2 system.

Annihilation cross sections
An initial calculation of the hadronic annihilation cross section using the basis set generated with the parameters given in table 1 gave a very small cross section. Gregory also found that using a basis set which was suitable to converge the elastic scattering cross section gave very small values for hadronic annihilation cross sections [33]. The delta potential method relies on the Kohn wavefunction being able to accurately represent the short-range interaction between the antiproton and the protons, especially for geometries where the particles coincide. To improve the trial basis Gregory and Armour added further functions which were designed to model the analytical solutions for protonium in the spherical prolate coordinates used for the scattering calculation [11]. This gave very large hadronic annihilation cross sections which at low energies behaved ass = -E 4.68 a pp 1 2 . The low hadronic annihilation cross section calculated using the Gaussian basis set described above was suspected to be caused by the inability of the basis to accurately describe the scattering wavefunction where the antiproton coincides with one of the protons. The parameters used to place the Gaussians only partially took this into account. They were sufficient to generate a basis for a converged elastic cross section but not to give an accurate wavefunction where the antiproton is very close to the protons.
To improve the basis set, Gaussians were generated using different parameters. As the singularities occur at θ=0 and π, (and = = R 0.700 55 re 2 a 0 ) the effect of decreasing θ R on the hadronic annihilation cross section was investigated. As θ R is reduced the radial basis is generated closer to the singularity giving Gaussians with larger radial exponents and thus more basis functions. This was found to increase the hadronic annihilation cross section. The effect of increasing the C R parameter was also investigated. Since the basis set was only required to be improved around the singularities, the value of C R was only changed for Gaussians in the range R=0.69-0.71 a 0 and set to 0.3 for all other R. Increasing C R was also found to increase the hadronic annihilation cross section. Changing other parameters such as increasing S R or C θ was not found to significantly increase the annihilation cross section. Around a 5% variation in the elastic cross section occurs as the basis is changed but the result is fairly stable.
The basis which gave the largest value of the annihilation cross section with q =´-1 10 R 6 and C R =0.9 and all other parameters as given in table 1 was used to calculate the annihilation cross section at different energies.  There are a number of reasons that suggest the hadronic annihilation cross sections calculated here are significantly too low. From test calculations the annihilation cross section does not appear to be converged with respect to the Gaussian basis set. Increasing the values of C R continued to increase the cross section significantly. Increasing C R further lead to difficulties in placing the Gaussian functions. The cross sections are also around a third of those of¯-H H scattering which at low energies behave ass » -E 0.14 a pp 1 2 [23]. From the appearance of the variationally calculated scattering wavefunctions it was seen that there is an increase in the magnitude of the wavefunction at positions where the antiproton coincides with the protons. Thus it appears that the Gaussian basis is reproducing approximately the correct shape of scattering wavefunction. The reason for the small hadronic annihilation cross sections calculated here is likely due to the delta potential method requiring very accurate wavefunctions at a single point. Although the Gaussian basis gives the correct overall shape, the actual amplitude of the wavefunction is too small at the Coulomb singularities.
The larger annihilation cross sections of Gregory and Armour seem likely to be more accurate than those calculated here and further work is needed to confirm their result. It appears that more specialist basis functions are required to account for the hadronic annihilation if treated using the Kohn method. Additional basis functions (such as those used by Gregory and Armour) are more difficult to implement than the Gaussian functions used here and are not as general. Further work is required to develop basis sets and methods to account for hadronic annihilation in a more general way for multidimensional systems such as¯-H H 2 .
The scattering wavefunction used to compute the hadronic annihilation cross sections in figure 8 was also used to compute leptonic annihilation cross sections. It was found that small values of R<1.5 a 0 did not significantly contribute to the leptonic cross section computed using equation (2.5). Similarly the contribution from angles near θ=0 and π was negligible. This suggests that unlike for hadronic annihilation, the leptonic cross section is much less dependent on the details of the scattering wavefunction when the antiproton is close to the protons. To test this, the basis set used for elastic scattering was also used to compute the leptonic annihilation cross section. The values for the two bases differed by 8%. The leptonic cross sections are thus much less sensitive to the details of the scattering wavefunction where the antiproton coincides with one of the protons.
The leptonic annihilation cross sections for two-photon and three-photon annihilation are shown in figure 9. At low energies the cross sections are given by 7.3×10 −5 -E 1 2 for two-photon annihilation and 6.4×10 −8 -E 1 2 for three-photon annihilation. These values are around twice as large as those calculated for the¯-H H system [26]. The reason for the larger leptonic annihilation cross sections is the larger values of P(R) for¯-H H 2 with two electrons compared to¯-H H. Repeating the calculation using P(R) for the¯-H H system [27] gave leptonic annihilation cross sections to within 10% of those for¯-H H.
It was pointed out by Froelich et al that if the leptonic annihilation is not a negligible channel then the BO approximation will break down as the leptons will no longer provide the potential for the hadrons [26]. However¯-H H 2 leptonic annihilation occurs with a very small probability and so the BO approximation should still be reliable.

Three-dimensional elastic scattering calculations
We attempted to carry out low-energy three-dimensional -H H 2 elastic scattering using the three-dimensional PES, that is with the r coordinate free to vary and H 2 vibrational basis functions up to v=2. However no indications of convergence were observed. There are no barriers in the¯-H H 2 PES and the adiabatic reaction to form PsH + Pn can occur even in the limit of zero collision energy. It is thus unlikely that a calculation not accounting for this reaction can be converged. The situation is similar to that described by Soldán et al for vibrationally inelastic collisions in the Na + Na 2 system [34].
In this case the Na atom can insert between the atoms in Na 2 . A full reactive calculation was required to calculate the inelastic cross-sections.

Conclusions
We have carried out ab initio calculations of the leptonic energies for the¯-H H 2 system within the BO approximation using a basis set of ECG functions. Energies were calculated at over 2800 geometries of the system. We have fit analytical functional forms to the ab initio energies to construct twodimensional and three-dimensional PESs, the first to be published for this system. The two-dimensional surface should be useful to develop and test scattering methods for the¯-H H 2 system. The three-dimensional PES should be of use for preliminary treatments of¯-H H 2 reactive scattering calculations. Further ab initio calculations for Pn-PsH type geometries would be useful to include in the future. We have carried out scattering calculations using the S-matrix Kohn variational method, treating the molecule as a rigid rotor and employing a basis set of two-dimensional Gaussian functions. These functions were placed to ensure sufficient overlap with neighbouring functions and with exponents tailored to the potential.
Elastic scattering calculations using the two-dimensional PES gave low-energy limiting elastic cross sections in fair agreement with those calculated by Gregory and Armour [11] but with a large discrepancy with those of Sultanov et al [12]. Elastic scattering cross sections computed using the threedimensional PES (with fixed H-H distance) were less than a factor of two different than the two-dimensional surface. Cross sections for rotational excitation and quenching seem to show a complicated structure.
Hadronic annihilation cross sections were computed from the variational scattering wavefunctions but were not converged despite improving the basis set. The cross sections for this process computed here are two orders of magnitude smaller than those of Gregory and Armour [11]. This appears to be due to the Gaussian basis functions not properly accounting for the proton-antiproton interaction at close distances. We have carried out the first calculation of leptonic annihilation for the¯-H H 2 system. The cross section for this process is around twice that for the¯-H H system but still remains small.
Further scattering calculations for the¯-H H 2 system are needed to take into account reactive channels and it appears that a three-dimensional scattering calculation can only be carried out with reaction included. The reaction to form Pn + PsH could possibly be treated on the ground state BO PES but other reactions such as Pn + Ps + H are open even at zero collision energy, so ultimately a non-adiabatic method may be required to obtain accurate cross sections.
The ab initio energies and fitted two-dimensional and three-dimensional PESs (the latter in the form of FORTRAN subroutines) and figures S1-4 mentioned in the text are available as electronic supplementary information.

*
Each 'M' matrix involves elements of the Hamiltonian and the basis functions used. The general form of the integrals for the Kohn method applied to rigid-rotor scattering is where j indicates the rotational channel which has internal energy ò j and orbital angular momentum l.
n j n j 0 , 0 The M matrix is a 'large' matrix consisting of elements between the basis functions: The integrals required were evaluated either analytically or numerically using Legendre and Laguerre-Gaussian quadrature.
From the S-matrix the scattering cross section is given by [36] Appendix B. Annihilation Following Jonsell et al [22] the hadronic annihilation rate is given by The flux is related to the normalisation of the asymptotic scattering wavefunction and is derived by making use of the plane wave boundary conditions. From equation (4.5) the asymptotic form of the incoming, J=0 elastic scattering wavefunction is given by is used giving a further factor of 2.
Within the BO approximation the leptonic annihilation rate is 9 E h a 0 3 for orthopositronium [26].
The full expression for the leptonic annihilation rate is given by where the factor of 2π comes from angular integration over the polar azimuthal angle f. The cross section is calculated by dividing the annihilation rate by the flux, F, as for hadronic annihilation scattering.
Appendix C. Potential energy surfaces

C.1. Ab initio calculations
The leptonic Schrödinger equation whereĤ lep is the leptonic Hamiltonian, has been solved variationally for a set of configurations of two protons and an antiproton. The wavefunction of the light particles has been expressed as the linear combination of K ECG functionŝ whereÂ is the antisymmetrizer, Θ is the spin singlet function of two electrons,P is the spatial symmetry projector, the vectors r i represent the lepton positions and C I , a i I , , b ij I , and R i I , are all variational parameters. The nonlinear basis function parameters, exponents a i I , and b ij I , and vectors R i I , , were optimised with the method of conjugated directions [37][38][39], independently in each point in the space of nuclear configurations. The optimisation of the nonlinear parameters was continued until the energy in a cycle was lowered by at most 30 E n h . This procedure is expected to give an accuracy of a few μE h . The initial form of the basis set is the Cartesian product of the basis consisting of 75 symmetrized ECG functions for the electrons in the H 2 molecule, and the positronic basis consisting of 10 Gaussian orbitals, centred on the antiproton.
The symmetry projectorP is a linear combination of symmetry operationsp k , that belong to a chosen finite point group.ˆˆ( a p 4.20 k N k k 1 sym N sym is the order of this group and the a k coefficients depend on the chosen irreducible representation, which-for practical reasons-has to be one-dimensional. For the problem considered here of interacting antihydrogen atom with hydrogen molecule both in their ground states, the leptonic wavefunction is always fully symmetric, resulting in a k =1 for each k. The hydrogen molecule and the linear configurations of all three nuclei have either ¥ D h or ¥ C v symmetries. They are represented by the C s and C 1 groups respectively, with all R i I , vectors fixed on the axis containing the nuclei. Triangular structures have the C 2v (T-shape) or C s point group symmetries. The basis set sizes are equal to 750 for the nuclear configurations having the symmetry plane perpendicular to the line connecting the protons ( ¥ D h and C 2v ) and 1500 for remaining configurations.

C.2. Analytical representation of PESs
The NN method requires no prior knowledge about the shape of the potential and fits are systematically improvable. These properties make the NN approach very useful for the unusual shapes of antimatter-matter potentials. The specific form of NN functions used here is based on the implementation of Manzhos et al [40] which was subsequently used by Law et al [41]. The functional form is a sum of sigmoid functions where w p are the NN weights, b p are biases and x is a vector of coordinates. Each term in the sum is called a node. The parameters w p , d p , b p and d 0 are determined by least-squares fitting to ab initio data using the I-NoLLS programme [42].

C.2.1. Two-dimensional PES
The internuclear distances r 1 , r 2 and r 3 (where r 3 is the p-p distance, which is fixed in the two-dimensional PES) were chosen as the coordinates in x in equation (4.21). The potential energy is identical for permutation of the two protons and these coordinates make this straightforward to take into account. To ensure the permutational symmetry, each node of the NN was used to generate another node with r 1 and r 2 exchanged. That is, equation ( [43]. Here we have fixed the proton-proton distance in our two-dimensional surface below to the equilibrium value for H 2 , 1.4011 a 0 , but not adjusted the C coefficients. In this formula R and θ are Jacobi coordinates with R the distance from the antiproton to the centre of mass of the two protons and θ the angle between the line joining the two protons and R. P n , are the Legendre polynomials [35]. Although the accuracy of the NN fit is good (see section 3.2.1), at large¯-H H 2 distances there are small oscillations of the potential energy. This deficiency in the NN fit is not important however as at large¯-H H 2 distances the NN is smoothly combined with the dispersion energy formula of equation (4.23) using a switching function (following Cvitaš et al [44]): where a affects how quickly the function switches from 0 to 1 and s affects where it switches. In this case a was set to 3.0a 0 1 and s to 7.0 a 0 so that the switch from the NN to longrange dispersion formula takes place around 7 a 0 where both functional forms are accurate. To avoid numerical problems of using equation (4.23) for small values of R, for R below 1.0 a 0 only the NN function is used.
The combined functional form of the two-dimensional -H H 2 PES is given in Jacobi coordinates by is the nuclear potential energy.
The energy of the separatedH and H 2 ,¯+ E H H 2 , is subtracted to give an interaction energy.

C.2.2. Three-dimensional PES
The isolated hydrogen molecule potential energy curve was fit to a functional form as it is used in the long-range three-dimensional¯-H H 2 function (see below). The PEC is also required for calculating the rotational-vibrational channel functions of the molecule for use in scattering calculations where the H 2 bond length is free to vary (along with the threedimensional PES of course). The p-p distance was varied with theH atom at R=20 a 0 and q = p 2 . A sum of powers of Morse variables [45] was used for the PEC fit: The nonlinear parameters c and d were fixed for each leastsquares fit attempted and so the optimum linear coefficients, C i , could be obtained in a single iteration of least-squares fitting.
The long-range functional form for the three-dimensional PES was based on equation (4.23) for the dispersion energy. For the three-dimensional surface the dispersion energy needs to vary also as a function of the p-p distance, r. It was found that simply multiplying the terms by r (and adjusting the C coefficients) was effective so that the long-range function was The energy of the H 2 molecule for a given p-p distance was calculated from the fitted PEC and the hydrogen atom energy added to this value. The dispersion terms then account for the interaction between theH atom and the hydrogen molecule. The full long-range functional fit in Jacobi coordinates is therefore The energy of interaction for geometries where the three atoms are considered separated was taken into account by adding the dispersion energies for each pair of atoms. This is given by the dispersion energy formula for two H atoms:  [43]. This approach ignores three body effects between the atoms but still gives very good accuracy at long range: the rms error is 0.69 mE h for the ab initio data for geometries where all atoms are over 8 a 0 apart. The accuracy of this approach improves the further apart the atoms are with the largest deviation of 2.97 mE h occurring for the geometry where the atoms were closest.
For fitting the H 2 molecule PEC the parameter d in equation (4.27) was set to 1.4011 a 0 , the molecule's equilibrium bond length, while the c parameter was optimised by manually changing its value for a given number of Morse functions to obtain the best fit. The fit was also constrained so that the sum of linear coefficients was equal to −1 E h so that for infinite separations the function reproduced the sum of two isolated H atom energies. The best fit was obtained with an eight-term sum of Morse functions with c set to 0.78a 0 1 . This gave a standard deviation of 0.073 E m h to the data with the largest difference of 0.13 E m h for the r=6 a 0 data point. As an additional check on the accuracy of the fitted PEC, a variational calculation of the vibrational energy levels was carried out. Using 8 harmonic oscillator basis functions the v=0 to v=1 vibrational transition was calculated as 4171 cm −1 . This is within 10 cm −1 of the experimental value of 4161 cm −1 [46].
The long-range functional form for the¯-H H 2 interaction energy was fit using equation (4.28). To determine the coefficients, this function was fit to 1214 ab initio energies with geometries where the p-p distance  r 6 3 a 0 and the Jacobi coordinate R between thep and centre of mass of the two protons >7 a 0 . This gave a standard deviation of 0.051 E m h with a maximum difference between the fit and data of 0.13 mE h .
The NN and asymptotic functions were combined in a similar way to the two-dimensional surface. In this case additional conditions are imposed depending on which geometry is required.
The potential energy is given as a function of the three internuclear distances. If all three of these distances are below 7 a 0 then the potential is given as a sum of the nuclear energy and the NN leptonic energy. This avoids numerical problems of trying to evaluate the asymptotic functions for geometries with small r i . If all internuclear distances are above 8 a 0 then the atoms are considered separated and the energy is given by the sum of dispersion energies equation (4.30). For¯-H H 2 geometries with the p-p distance less or equal to 6 a 0 and the two¯-p p distances greater than 7 a 0 , the potential is given as a sum of the NN energies and¯-H H 2 asymptotic form, equation (4.29), with a switching function to smoothly join these together. The switching function, S(R), is the same as used for the twodimensional surface, equation (4.24), with a set to 3.0a 0 1 and s to 8.0 a 0 . For all other geometries the NN energy is used.
The full expression is given by