Correspondence between the surface integral and linear combination of atomic orbitals methods for ionic-covalent interactions in mutual neutralisation processes involving H$^-$/D$^-$

The surface integral method for estimating ionic-covalent interactions in diatomic systems been successful in producing cross sections for mutual neutralisation (MN) in reasonable agreement with experimental results for branching fractions between final states in systems such as O$^+$/O$^-$ and N$^+$/O$^-$. However, for simpler cases of MN involving H$^-$ or D$^-$, such as Li$^+$/D$^-$ and Na$^+$/D$^-$, it has not produced results that are in agreement with experiments and other theoretical calculations; in particular, for Li$^+$/D$^-$ calculations predict the wrong ordering of importance of final channels, including the incorrect most populated channel. The reason for this anomaly is investigated and a leading constant to the asymptotic H$^-$ wavefunction is found that is different by roughly a factor $1/\sqrt{2}$ to that which has been used in previous calculations with the surface integral method involving H$^-$ or D$^-$. With this correction, far better agreement with both experimental results and with calculations with full quantum and LCAO methods is obtained. Further, it is shown that the surface integral method and LCAO methods have the same asymptotic behaviour, in contrast to previous claims. This result suggests the surface integral method, which is comparatively easy to calculate, has greater potential for estimating MN processes than earlier comparisons had suggested.

The surface integral method for estimating ionic-covalent interactions in diatomic systems been successful in producing cross sections for mutual neutralisation (MN) in reasonable agreement with experimental results for branching fractions between final states in systems such as O + /O − and N + /O − . However, for simpler cases of MN involving H − or D − , such as Li + /D − and Na + /D − , it has not produced results that are in agreement with experiments and other theoretical calculations; in particular, for Li + /D − calculations predict the wrong ordering of importance of final channels, including the incorrect most populated channel. The reason for this anomaly is investigated and a leading constant to the asymptotic H − wavefunction is found that is different by roughly a factor 1/ √ 2 to that which has been used in previous calculations with the surface integral method involving H − or D − . With this correction, far better agreement with both experimental results and with calculations with full quantum and LCAO methods is obtained. Further, it is shown that the surface integral method and LCAO methods have the same asymptotic behaviour, in contrast to previous claims. This result suggests the surface integral method, which is comparatively easy to calculate, has greater potential for estimating MN processes than earlier comparisons had suggested.

I. INTRODUCTION
Charge transfer, where an electron moves from one atom or ion to another during a collision, is a fundamental atomic process. Mutual neutralisation (MN), with corresponding reverse process ion-pair production is an important special case, and methods for estimating cross sections are important for interpreting experimental results and testing our understanding of the physical mechanisms involved [1][2][3], as well as modelling various plasma environments such as planetary and planetary satellite atmospheres [4] including the earth's ionosphere [5], stellar atmospheres [6,7], as well as being important for diagnostics in and fusion plasma applications and high energy physics experiments employing negative ion sources [8]. The transferred electron is captured into specific excited states A * , and thus affects the distribution of populations of states of atom A and its observed spectrum. These processes occur predominantly through interactions between ionic (A + +B − ) and covalent (A * +B) configurations at avoided crossings in the adiabatic potential energy curves e.g. [9,10]. Any method for estimating MN cross sections thus hinges on the ability to calculate the ionic-covalent interactions. A promising approach is the surface integral method. The general technique was independently developed roughly simultaneously by a number of workers, and thus is known by various names including the Holstein-Herring method, e.g. [11], the Firsov-Landau-Herring method, e.g. [12], and the Landau-Herring method, e.g. [13][14][15]. The original papers on the idea include those by Firsov [16], Holstein [17], and Herring [18], and the method was used to solve problems in Landau and Lifshitz's text Quantum Mechanics [19], published in Russian in the 1950's; the asymptotic (large internuclear distance R) splitting between gerade and ungerade states for H + 2 is solved in §81 of Ref. [19], building on problem 3 in §50. The basic idea is that the interaction energy can be calculated from the electronic current flowing across a surface between the two nuclei, obtained by manipulating Schrödinger equations corresponding to cases where the electron is localised on either nuclei. On the assumption that the wavefunctions are close to zero at the surface, i.e. in the limit of large internuclear separations, the method allows the interaction to be derived in terms of the atomic electronic wavefunctions at the surface.
A key point in the use of the method is the choice of the integration surface to ensure that wavefunctions are small. Smirnov [20] derived expressions for the interaction in the A + + e + B system, assuming the interaction between e and B to be zero, and taking the integration surface as the midplane, halfway between the two nuclei. Janev & Salin [14,15] instead account for the potential due to B, and choose the integration surface as a sphere around B where the potential becomes suitably small; see also Ref. [13] for an earlier application in resonant charge transfer.
Janev & Salin's expressions for the couplings from the surface integral method, which for consistency with previous work [21,22] [26]. In Refs. [21,22] the method was compared with full quantum and linear combination of atomic orbitals (LCAO) calculation results for MN of Li + , Na + , and Mg + with H − , and the results found to differ significantly. Recent experi-mental results for branching fractions in Li + /D − [1,2] and Na + /D − [3] have allowed calculations to be tested, and full quantum and LCAO results [7] have shown good agreement, implying the LHJ method is in disagreement with these experiments. This is puzzling, given that the LHJ method has been shown to produce estimates in reasonable agreement with experiment for MN of N + /O − and O + /O − [27,28].
Asymptotic methods such as LHJ and LCAO are of significant importance for many applications, as estimates can be made rather inexpensively compared to full quantum calculation (quantum chemistry potentials and couplings with quantum scattering), and thus it is valuable to resolve the origin of this discrepancy. The fact that the LHJ method appears to work well for MN involving O − , and less well for the more simple H − , would seem to suggest the following possible explanations: 1) there is a problem with the H − wavefunction used, 2) the agreement for O − is fortuitous, or 3) the full quantum, LCAO and experimental results are in error. In this paper it will be shown that an error appears to have been made in deriving the asymptotic H − wavefunction used in all applications to MN involving H − , and correcting this error brings the LHJ method into reasonable agreement with other theory and experiments. Previous claims that the discrepancy between the LHJ and LCAO methods is due to fundamental problems in the LCAO method are reexamined in light of this.

II. THEORY
Mutual neutralisation occurs predominantly at avoided crossings between adiabatic molecular potential curves, corresponding to real crossings in a diabatic representation, e.g. [9,10]. In a diabatic representation the states of the system can be described in terms of ionic and covalent configurations. The ionic configuration is the case where the active electron is located on core B, that is A + + (e + B), labelled i with corresponding diabatic electronic wavefunction Φ i , and the covalent configuration the case where the active electron is located on the core A + , that is (A + + e) + B, labelled c with wavefunction Φ c ; see Fig. 1 for a sketch of the system. These two configurations have corresponding diabatic states with potentials that cross, H ii = Φ i |H|Φ i and H cc = Φ c |H|Φ c where H is the Hamiltonian, and have an off-diagonal interaction (coupling) H ic = Φ i |H|Φ c . In the adiabatic representation, the potentials avoid crossing, and a key quantity is the splitting between the adiabatic potential curves ∆U , as it enters the Landau-Zener formula for the transition probability [9,29,30]. The transition probability between adiabatic states is (in atomic units, used throughout) where v R is the radial component of velocity of relative motion of the nuclei, and R is the internuclear distance and R x the distance at which the diabatic potentials H ii and H cc cross. The splitting ∆U is related to the interaction H ic = Φ i |H|Φ c in the (non-orthogonal) diabatic basis (e.g. [31,32]) (3) where S ic = Φ i |Φ c is the overlap between states i and c. Note that if the diabatic basis is orthogonal References given above, as well as for example Refs. [10,21,33], can be consulted for more details on the Landau-Zener model, and the relationship between adiabatic and diabatic representations, including orthogonal and nonorthogonal representations.

A. Surface integral method
The surface integral method can be used to estimate the quantity ∆ via where s is the surface of integration; see Refs. [12,14,15,31] for more details. Note the wavefunctions Φ are here only the spatial components. The LHJ theory, choosing the surface s to be a sphere centred on B, derives an analytic expression for ∆(R x ) written in terms of the asymptotic forms of the spatial atomic wavefunction ϕ of the active electron. That is such that R nl (r) is the radial part of the wavefunction, Y lm is the spherical harmonic function, and r is the active electron coordinate with respect to the relevant nucleus ( r A or r B ). When the active electron is on core A + , i.e. the system is in the covalent configuration c, the asymptotic behaviour of the Coulomb wavefunction (i.e. the leading term of the asymptotic expansion) is used where −γ 2 c /2 is the binding energy of the electron, and the quantum numbers nl are now dropped from the radial function R for simplicity. When the electron is on core B, thus forming a negative ion, and the system is in the ionic configuration i, the appropriate asymptotic form is used where −γ 2 i /2 is the binding energy of the electron, and thus E − B = γ 2 i /2 is the electron affinity of the negative ion, k l (x) is the modified spherical Bessel function (e.g. [34]), and l i the orbital angular momentum quantum number. For the specific case of l i = 0, Both wavefunctions contain leading constants, N c and N i , that must be determined. In the case of the N c , this can be estimated from the Coulomb wavefunction normalisation constant from quantum defect theory [35][36][37] (see, e.g., eqn. 19 of Ref. [14]), which can be expected to give accurate results for sufficiently excited states with diffuse wavefunctions well described by quantum defect theory. The leading constant factor for the negative ion function is, however, more problematic to determine. Due to the 1/r factor, the normalisation is dominated by contributions at short distances where the asymptotic form is not valid and overestimates the wavefunction, and the value given by normalisation N i = √ 2γ i (see Ref. [38]) cannot be expected to give an accurate representation of the asymptotic behaviour. N i should therefore be determined by comparison with detailed atomic structure calculations. One expects N i > √ 2γ i , and this has sometimes led to the use of a cutoff value r 0 , below which the wavefunction is zero and can be chosen to give correct normalisation e.g. Ref. [21,38].
Using these wavefunctions, Refs. [15,39] give the expression for ∆(R x ) as which for l i = 0, as in the case of H − ( 1 S), reduces to the rather simple expression The modification due to angular momentum coupling to core electrons in complex systems, including those with equivalent electrons, is given in Refs. [12,15,39]. The major advantage of this approach is that, as can be seen, it reduces to an analytic expression, and does not require the evaluation of any matrix elements of the Hamiltonian.

B. LCAO method
The LCAO approach can be used to calculate ∆(R x ) through calculation of the matrix elements of the Hamiltonian using the same or similar atomic wavefunctions for the active electron. This has been done in Refs. [21,32,40,41] for the case involving H − , considering two electrons. The main differences between these different descriptions are the treatment of the two-electron wavefunctions, leading to slightly different expressions, but they are equivalent asymptotically (at large separations R). For the purposes of our discussion, a much simpler one-electron version is preferable, again equivalent to the above formulations at large R.
Considering the model and coordinate system given in Fig. 1, ignoring any interactions with neutral atom B, the Hamiltonian is H = −1/r A , and we have the rather simple expression in terms of atomic wavefunctions which for large R we have and similarly Using these results in eqn. 3, we obtain where we now define the "post-interaction operator" T of charge transfer theory also found in Refs. [32,40]. The same result is found through a different derivation in § 10.3 of Ref. [42]. Analytic expressions for T ic (R) can be obtained adopting the asymptotic forms of the atomic wavefunctions given above, though for any given l c and l i , they are significantly more complex than those arising from the surface integral method for ∆(R x ) given above. Expressions both for T ic (R) and S ic (R) for the l c = l i = 0 case are given in the Appendix.

C. H − asymptotic wavefunction
The general theory of asymptotic wavefunctions has been well expounded, for example in § 2.1 of Ref. [43], as well as Ref. [44]. The derivation of the asymptotic wavefunction for H − is particularly simple, and so it is instructive for the following discussion to present the basic equations; see also e.g. Refs. [41,42,44] for details. The two-electron non-relativistic Hamiltonian, in atomic units, is where r 1 and r 2 are distances of the electrons with labels 1 and 2 to the proton, and r 12 is the distance between electrons. Taking the asymptotic case r 1 r 2 , then r 12 ≈ r 1 , the asymptotic form then allows separation into terms related to each electron such that H 1 is the Hamiltonian for the loosely bound distant electron, and H 2 is the Hamiltonian for the tightly bound electron close to the proton. Then the timeindependent Schrödinger equation becomes where ψ is the electronic wavefunction, E 1s H is the energy of the hydrogen ground state (−0.5 au) and E H − is the electron affinity of H (a positive value, see below). Separating the total wavefunction ψ into spatial and spin functions ψ = ϕ( r 1 , r 2 )χ S=0 , where χ S=0 is the singlet spin function, then the asymptotic spatial wavefunction becomes where H 2 ϕ H 1s ( r 2 ) = E 1s H ϕ H 1s ( r 2 ) and thus ϕ H 1s ( r) = exp (−r)/ √ π is the hydrogen ground state function, and the long range (LR) electron function is the solution of The spherically symmetric (l = 0) solution is where γ = √ 2E H − , and A is an arbitrary constant related to N i above in eqn. 9, which we need to determine. The electron affinity of H from modern measurements is E H − = 0.754195(19) eV or 0.0277162 atomic units [45], which gives γ = 0.235441 in atomic units. The validity condition γr 1 thus requires r 1/γ = 4.26 au. As will be shown below, A = N i 1 √ 4π and correct normalisation would imply N i = √ 2γ ≈ 0.686 and A ≈ 0.194. Unfortunately a range of different definitions and notations for the asymptotic function and the leading constant have been used in various sources. Part of this difference arises often from separation into radial and angular components, for this case (l = 0) In early papers by Smirnov [20], the following definition was used R(r) = A S65 2γ exp (−γr)/r (24) where the √ 2γ arises since the function is correctly normalised if A S65 = 1 (see § II A). In Ref. [20], the value A 2 S65 = 2.65 was derived, which means A S65 = 1.63. This implies A = A S65 √ 2γ 1 √ 4π = 0.316. Later work by Smirnov and others [14,23,38,42,46] have adopted the definition R(r) = A S exp (−γr)/r = B S 2γ exp (−γr)/r (25) such that B in Ref. [42], here denoted B S , is equal to A S65 . Note, A S = N i and is thus the value required for use in the equations of Janev and Salin [14,15], that is, equations 10 and 11. Smirnov [38,42,46] has derived A S from the simple Chandrasekhar 2-and 3-parameter wavefunctions [47], finding a value of A S = 1.13 ± 0.06, and thus B S = A S65 = 1.64 ± 0.08. A value of A S in this range has been used in all applications of the LHJ method to processes involving H − [21][22][23][24][25][26]. This gives In LCAO applications [40,41] the value A = 0.223106 has been used, determined from the 203-parameter variational wavefunction of Pekeris [48], asymptotic values of the wavefunction being presented in Ref. [49]. This derivation will be discussed below in § II C 1. In Ref. [22] the LCAO method was applied, and the same value is essentially used, but a factor of √ 2 must be applied to adjust for the antisymmetrized version of the asymptotic wavefunction that was employed.
Clearly, there is a discrepancy between the leading constant of the asymptotic wavefunction derived by Smirnov and used in all LHJ theory applications, A ≈ 0.319, and the value that has been applied in LCAO calculations, i.e. A ≈ 0.223. The two values of A differ by a factor of roughly √ 2. In the following we repeat both derivations in detail to show the source of the problem.

Direct matching of the asymptotic wavefunction
The most straightforward approach to determining the appropriate constant in the asymptotic wavefunction is through direct matching to a calculated wavefunction for H − . This is the approach that has been used in the LCAO case, and the method stems from Ref. [49]. They define a function equal to the two-electron spatial wavefunction where either r 1 = 0 or r 2 = 0, ψ(r) in their notation. Here we denote this same function: ϕ 0 (r) ≡ ϕ( r 1 = 0, | r 2 | = r) = ϕ(| r 1 | = r, r 2 = 0) = C(r) exp (−γr)/r.
such that C(r) is the same as in Ref. [49], and converges to a constant value C(∞) at large r [50]. Equating to the asymptotic wavefunction from above, eqn. 20, taking r 2 = 0, we obtain which means A = C(∞) √ π. In Ref. [49], C(∞) = 0.125874 and thus A = 0.223106 was found by matching points at r = 14-16 atomic units from the calculated values of ϕ 0 (r) extracted from the 203-parameter wavefunction of Ref. [48] given in their Table 1. This is the value used in LCAO calculations. Figure 2 compares ϕ 0 (r) extracted from various detailed calculations of the wavefunction ϕ( r 1 , r 2 ), with asymptotic forms of the wavefunction ϕ asymp 0 (r) for the two different values of A discussed above. It is clearly seen that A ∼ 0.223 matches all the functions better than the higher value of A ∼ 0.319. Note, among detailed calculations, only the Pekeris wavefunction has the expected asymptotic form even at very long range, since the variational form of the wavefunction used ensures correct asymptotic behaviour. The 2-and 3-parameter wavefunctions derived by Chandrasekhar [47] are of particular interest as they were used by Smirnov, but also because they are so simple they can be used to exemplify various issues. The Chandrasekhar wavefunction is: ϕ( r 1 , r 2 ) = N (exp (−ar 1 − br 2 ) + exp (−ar 2 − br 1 )) × (1 + cr 12 ) (28) where for the two parameter function (c = 0) he finds through variational calculations a = 1.03925 and b = 0.28309, and it can be found from the normalisation condition that N = 0.031443, and for the three parameter function introducing correlation he finds a = 1.07478, b = 0.47758, c = 0.31214, and it can be found that N = 0.031226. That a ≈ 1 in both cases shows the expected result of an unscreened hydrogen 1s-like function for the electron close to the proton, with a more distant long range electron. Having said that, it is clear that the functions do not have the correct asymptotic behaviour when the distant electron is very far from the nucleus. The Chandrasekhar 2-parameter function does not have asymptotic behaviour at any distance, and the 3-parameter and Hart and Herzberg 20-parameter function [51] have only roughly asymptotic behaviour in the region r ≈ 4-10 au. This is due to the fact that the asymptotic region has little influence on the state energies and thus the wavefunction determined by variational methods.

Electronic density matching method
An alternative approach, and that which has been used by Smirnov [38,42,46], is to match the electron density The asymptotic H − wavefunction ϕ asymp (r) as defined in eqn. 26, extracted from various detailed calculations of the wavefunction ϕ( r1, r2): the 203-parameter function from Pekeris [48] calculated in Ref. [49], the 19-parameter function of Hart and Herzberg [51], the 2-and 3-parameter functions of Chandrasekhar [47]. Asymptotic forms (eqn. 27) with the two values of A discussed in the text are also shown.
function. Generally for an N -electron system, and ignoring spin, the electronic density is and thus for H − and two electrons ρ( r) = 2 d r 2 |ϕ( r, r 2 )| 2 .
For the asymptotic form this gives: As done by Smirnov, one can derive analytical expressions for the asymptotic constant A or some related version of it, e.g. A S or A S65 , as a function of r, by equating the asymptotic form with the electronic density calculated from Chandrasekhar's wavefunctions, and the expressions found by Smirnov are given in Refs. [38,42]. Using Mathematica this analysis was repeated here, making use of Hylleraas coordinates for the 3-parameter case [52,53], and the results are plotted in Fig. 3. factor of 1/ √ 2 lower than found by Smirnov A ≈ 0.32 (A S ≈ 1.13). The difference is not precisely a factor of 1/ √ 2 since different wavefunctions are used, the direct method permitting the easy use of the accurate 203parameter function of Pekeris.
Note, the criterion for validity of the asymptotic wavefunction, r 1/γ = 4.26 a.u., is not satisfied across this entire region r = 2-5 a.u. Further, it is clear that the Chandrasekhar wavefunctions do not have correct asymptotic behaviour, as is seen from application of the direct method, which gives rapidly varying values of A (see also Fig. 3). That the density matching method using these simple wavefunctions gives values that agree reasonably with the values from direct matching of the accurate and asymptotically correct wavefunction of Pekeris, is thus perhaps somewhat surprising. It may reflect that the density method, by accounting for contributions from both electrons at intermediate distances, reaches asymptotic behaviour at smaller internuclear distances than the direct method where the second electron is fixed at the nucleus. Contributions from both electrons to the wavefunction are equal, ϕ LR ( r) = ϕ H 1s ( r), at roughly r = 2.3 a.u.

III. COMPARISON OF THEORETICAL CALCULATIONS AND EXPERIMENTS
Experiments have recently been performed in Louvain on MN of Li + /D − [1] and at the DESIREE facility in Stockholm on Li + /D − [2] and Na + /D − [3], resolving final states and thus measuring branching fractions for the neutral products. D − is preferred over H − in the experiments for practical reasons, but is basically identical to H − in terms of binding energy (electron affinity) and thus electronic structure. The different mass leads to trajectory (Coulomb focussing) effects, but is easily accounted for in the dynamical calculations. In a recent paper [7], existing full quantum calculations [54][55][56] and (two-electron) LCAO method calculations [21,22] were compared with experiment, generally finding good agreement, with the full quantum calculations performing best, in line with expectations. The reader is referred to Ref. [7] for a discussion regarding these comparisons, including the discussion of Coulomb focussing effects caused by H − versus D − .
Here, these comparisons are supplemented with calculations from the LHJ theory using the two A values discussed, and these calculations use the same code as the LCAO calculations, described in Refs. [21,22]. The results are shown for Li + /D − in Fig. 4, and for Na + /D − in Fig. 5. To aid comparison (see Ref. [7]), the plots are presented on the reduced energy scale, which is defined as where E CM is the collision energy in the centre-of-mass frame, µ is the reduced mass of the system, and v is the relative velocity. It is clear that the comparisons with experiment are significantly improved for both Li + /D − and Na + /D − if the low value of A = 0.223 is adopted, compared to the high value of A = 0.319. In particular we note that for Li + /D − (Fig. 4), the high value leads to the prediction of 3p being the dominant channel at low energy, rather than 3s as clearly seen to a high degree of certainty in experiments. This discrepancy is resolved for the LHJ theory with the lower value of A. Further the low value of A leads to much better consistency between full quantum, LCAO and LHJ predictions in both cases. Note also that LHJ with the high value of A leads to variation of the branching fractions at lower collision energies. Janev and Radulovic [23] also performed calculations with the LHJ method, for Li + /H − and Na + /H − , though state-resolved results are only shown for Na + /H − (Fig. 1 of Ref. [23]). Our results for Na + /H − partial cross sections are significantly different from theirs; in particular at low energy they find 3p to be the second most populated channel, only an order of magnitude lower than 4s. This is in strong disagreement with experiment where the 3p channel is not observed, and with all other calculations, including our LHJ calculations with the old, high value of A, where we find the cross section for 3p to be more than five orders of magnitude lower than 4s; this may suggest a labelling error. In Ref. [56] substantial differences were also noted compared to full quantum calculations. This was investigated further by comparing our LHJ theory avoided crossing parameters with those in Table 1 of Ref. [23], in particular ∆(R x ). However, for both Li and Na we were unable to reproduce all of their results using the equations and numbers in the paper. Using their equations, their value for the leading constant in the negative ion radial wavefunction, and their values for R x we obtain ∆(R x ) in agreement (within 10%) for a few crossings (Li 2s and 2p, Na 3p), but for the majority there are significant differences (up to a factor of 2.5 for Na 4p), usually our result being larger (only for Li 3d do we obtain a roughly 20% smaller value). Our LHJ theory values using A = 0.319 agree well with what we calculate from the equations in Ref. [23], which are just simpler versions of the more general theory. As such, detailed comparison of the avoided crossing parameters is not particularly enlightening on the more general problem of interest here, namely the relationship between LHJ theory, LCAO theory, and experiment. Branching fractions for the MN reaction Na + + D − → Na(nl) + D as a function of reduced collision energy. The 4s, 3d and 4p channels are shown in separate panels. Experimental results from DESIREE [3] are shown, with estimated errors (1σ). Theoretical results are shown from full quantum (FQ) calculations, LCAO and LHJ asymptotic methods with the old (high) value for A, and the new (low) value for A; see text for further details. Theoretical results are also shown for the case where D − is replaced by H − , marked by (H) in the legend. semi-empirically found a similar need for a factor 1/ √ 2 correction to the couplings in LHJ theory, and in doing so obtain good matches to experimental results for MN of D − with He + , Li + , Na + , C + , N + , though with some discrepancies for O + .
Finally, it is worth noting that calculations with the Smirnov formulation of the surface integral method [20], taking the surface of integration as the midplane between the two atoms, while improved with the lower value of A, still perform substantially worse than the Janev and Salin formulation; see also Ref. [26].

IV. DISCUSSION
It has been claimed in Ref. [39], that the LCAO method is "essentially incorrect for the problem of calculation of ∆(R)", with reference to the work of Herring [18]. It is, however, unclear that the criticisms of the LCAO method for the cases discussed by Herring, the exchange interaction [57] in H + 2 and H 2 , apply also to ionic-covalent interactions. First, the problem of the 1/r 12 term described for the exchange interaction in H 2 is not relevant, as it does not enter the one-electron LCAO method (see S II B) and in the two-electron methods only enters empirically via the H − electron affinity and wavefunction (e.g. eqn 10 of Ref. [21]). Second, in both cases, a basis with a single tightly bound 1s wavefunction is used in the LCAO description. In the case studied in detail by Herring of the exchange splitting in H + 2 , the surface integral approach permits a modification factor which increases the value of the wavefunction between the two nuclei, capturing departures of the molecular wavefunction from the atomic ones, which cannot be modelled in a LCAO basis with a single fixed orbital. In contrast, the LCAO method for ionic-covalent interactions employs two wavefunctions that capture the main components of the molecular wavefunction at internuclear distances corresponding to the crossing.
More quantitatively, in Ref. [39], it is claimed that there is a discrepancy between the asymptotic behaviour of the interactions in the LHJ and LCAO formulations. For LHJ the interaction has the asymptotic behaviour at large R x , see eqns. 10 and 7, while for LCAO at large R, Ref. [39] claims the asymptotic behaviour is however, this relationship is asserted without proof, and misses the important fact that γ i and γ c are related at the crossing point R x . Assuming a pure Coulomb interaction for the ionic potential, and a null potential for the covalent state, the crossing point is given by Using the one-electron LCAO expression for T ic (R) for l c = l i = 0 given in the Appendix and making the substitution γ i → γ 2 c /2 − 2/R x , an analytic expression for ∆ LCAO (R x ) = T ic (R = R x ) can be obtained. This results in a complicated expression, which is no longer dependent on γ i , and which can be shown to have the same behaviour as the LHJ theory for large R x ; i.e., for l c = 0 This agreement is most convincingly shown by numerical results for the ratio ∆ LCAO /∆ LHJ , which is dependent only on γ c and R x , shown in Fig. 6. The condition for validity of the LHJ method is R x γ −2 i (Ref. [14]), which, noting that for large R x we have γ i ≈ γ c − 1/(γ c R x ), is roughly equivalent to R x γ −2 c . The ratio is seen to be exactly one in the region of validity of the LHJ theory. Note the effective principle quantum number of the neutral state on A is n * = 1/γ c , thus the validity criterion can also be written R x n * 2 . Outside the region of validity for the LHJ theory, at crossings at shorter internuclear distances, the LCAO theory gives larger interactions, and the disagreement increases very rapidly.
In summary, previous claims that LCAO theory is incorrect for ionic-covalent interactions at large internuclear distance seem to be unfounded, and correction of the leading constant for the H − asymptotic wavefunction in the LHJ theory bring the two theories into reasonable quantitative agreement for the two cases of MN of Li + and Na + with H − and D − . It has also been shown that the LCAO and LHJ theory have the same behaviour for large R x . The LCAO theory has only been developed for cases involving H − and D − , treating two electrons to varying degrees of complexity. The surface integral method, LHJ, is however easily applied to any given situation if the asymptotic wavefunction for the negative ion is known. The results shown here suggest that if correct wavefunctions for the active electron on the negative ion can be obtained, LHJ theory can be used to obtain reasonable estimates for MN processes occurring at large internuclear distance in practically any system, perhaps even those involving molecules, and thus is of potentially great utility for situations such as astrophysical modelling, where often a large number of rate estimates are needed, and completeness can be more important than accuracy. It would be useful also to further develop and investigate the utility of the one-electron LCAO method as outlined in this paper, enabling application to complex atoms beyond H − and D − . ACKNOWLEDGMENTS I thank Xavier Urbain, Arnaud Dochain, Jon Grumer, Anish Amarsi, and Gustav Eklund for useful discussions, sharing preliminary results, and comments on the manuscript. The work relies on experimental results obtained at the Swedish National Infrastructure, DESIREE (Swedish Research Council Contract No. 2017-00621), originally published elsewhere, and I thank the DESIREE group at Stockholm University for collaboration on the experiments. This work is part of the project "Probing charge-and mass-transfer reactions on the atomic level", supported by the Knut and Alice Wallenberg Foundation (2018.0028), and supported by the Swedish Research Council through individual project grants with contract Nos. 2016-03765 and 2020-03404.