Mutual neutralization in low energy H^+ + H^- collisions
aa r X i v : . [ phy s i c s . c h e m - ph ] F e b Mutual neutralization in low energy H + + H − collisions:a quantum ab initio study Michael Stenrup ∗ Department of Theoretical Chemistry, School of Biotechnology,Royal Institute of Technology, SE-106 91 Stockholm, Sweden ˚Asa Larson and Nils Elander
Department of Physics, Stockholm University, AlbaNova University Center, SE-106 91 Stockholm, Sweden (Dated: May 29, 2018)The mutual neutralization of H + and H − at low collision energies is studied by means of amolecular close-coupling approach. All degrees of freedom are treated at the full quantum level alsotaking into account the identity of the nuclei. The relevant Σ + g and Σ + u electronic states as well asthe associated non-adiabatic radial couplings are calculated for internuclear distances between 0.5and 50 a . Following a transformation into a strictly diabatic basis, these quantities enter into aset of coupled equations for the motion of the nuclei. Numerical solution of these equations allowsthe cross sections for neutralization into the H(1) + H( n ), n = 1 , , PACS numbers: 34.70.+e, 31.50.Df, 31.50.Gh
I. INTRODUCTION
The mutual neutralization of H + and H − is the pro-totype reaction for electron transfer between two oppo-sitely charged ions. With only two protons and two elec-trons involved, it is the simplest ion-ion reaction possi-ble. Here, calculations should be able to provide accurateresults which could serve as benchmarks in the develop-ment of new theoretical methods and models. Further-more, the reactants have well defined quantum states andare relatively easy to prepare in the laboratory. A de-tailed comparison between experiment and theory shouldthus be possible and should provide a natural startingpoint to understanding more complicated ion-ion colli-sion phenomena. However, the simplicity of this reactionis deceptive, a fact which has stimulated a large num-ber of studies in the past. Nevertheless, surprisingly fewof these have investigated the neutralization cross sec-tion for collision energies below a few eV, which playsan important role in determining the H − abundance inlow temperature ionized environments. A compelling ex-ample from astrophysics is its influence on the gas phaseformation of H , in particular in the chemistry of theprimordial gas, where one of the main reaction paths in-volves H − as an intermediate [1, 2, 3].Bates and Lewis [4] were the first to study the mutualneutralization of H + and H − , which is written schemat-ically as H + + H − → H(1) + H( n ) , ∗ [email protected] where n is the main quantum number of the excited hy-drogen atom. Their calculations were based on the semi-classical Landau-Zener model [5, 6] in which the reactionis thought to proceed along a system of horizontal co-valent potential energy curves crossed by an attractiveion-pair potential. This pioneering study was later fol-lowed by others, among them by Olson et al. [7] and morerecently by Eerden et al. [8]. Fussen and Kubach [9] havegone beyond the Landau-Zener concept and investigatedthe reaction at low collision energies using a one electronclose-coupling model.The only measurement of the neutralization cross-section for collision energies below 3 eV is the mergedbeam experiment of Moseley et al. [10]. Their resultsseem to indicate a cross section that is approximatelya factor of three higher than obtained in most of thecited theoretical studies. However, the theoretical resultsare mutually consistent and, for higher energies, also inmore or less good agreement with the measurements re-ported by Szucs et al. [11] and Peart and Hayton [12].This intriguing situation has led to the suspicion thatthe low energy cross section of Moseley et al. is overesti-mated [11, 12, 13, 14].In this paper we report results of a fully quantum me-chanical ab initio study of the reaction at hand. We aremainly concerned with the collision energy region 0.001to 10 eV, but the calculations have been extended up to100 eV in order to facilitate comparison with the citedexperiments. The present study is based on a molecu-lar close-coupling expansion of the total wave function.Accurate electronic structure methods have been used tocalculate adiabatic states and non-adiabatic couplings ofboth gerade and ungerade symmetry over a wide rangeof internuclear distances. Subsequently, a transforma-tion into a strictly diabatic representation has been per-formed and the diabatic quantities entered into a set ofcoupled equations for the motion of the nuclei. By nu-merically solving these equations, we have been able tocalculate the neutralization cross sections for scatteringwithin both the gerade and ungerade inversion symme-tries. Finally, the identity of the nuclei has been ac-counted for by properly combining the gerade and unger-ade cross sections.The paper is organized as follows. Section II describesthe theoretical framework underlying the present study.Section III describes the electronic structure calculationsand discusses the resulting potential energy curves andcoupling matrix elements. Section IV gives a brief ac-count of the numerical treatment of the scattering prob-lem. Section V presents and discusses the calculated neu-tralization cross sections and the reaction rate coefficient.Section VI concludes the paper and summarizes the mainresults. Unless stated otherwise, vector and matrix quan-tities are written in boldface and atomic units are used. II. THEORYA. The molecular expansion: adiabatic anddiabatic formulations
Consider a diatomic system composed of two nucleiand N electrons. The masses and charges of the nucleiare denoted by M i and Z i , respectively. Let R = ( R, θ, ϕ )be the internuclear separation vector and r i a vectorpointing from the center of mass of the nuclei to the i -th electron. All vectors are expressed with respect toa space fixed axis system. For convenience, the notation r = ( r , r , . . . , r N ) will occasionally be used. In thesecoordinates, the center of mass motion of the total sys-tem can be separated out, and the Hamiltonian takes theform [15] H = − µ ∇ R + H el , (1)where µ is the reduced mass of the nuclei and H el isthe electronic Hamiltonian. Neglecting terms of the form[2( M + M )] − ∇ r i · ∇ r j , it follows that H el = − ∇ r + Z Z R + N X i>j =1 | r i − r j |− N X i =1 Z (cid:12)(cid:12)(cid:12) µM R + r i (cid:12)(cid:12)(cid:12) + Z (cid:12)(cid:12)(cid:12) µM R − r i (cid:12)(cid:12)(cid:12) . (2)The total wave function ψ ( R , r ) satisfies the time-independent Schr¨odinger equation Hψ ( R , r ) = Eψ ( R , r ) , (3) where E is the total energy in the center of mass system.We expand ψ ( R , r ) in terms of electronic and nuclearstates, φ a i ( R , r ) and χ a i ( R ), according to ψ ( R , r ) = ∞ X i =1 φ a i ( R , r ) χ a i ( R ) = φ a χ a , (4)where φ a i ( R , r ) are normalized solutions to the electronicSchr¨odinger equation H el φ a i ( R , r ) = ǫ a i ( R ) φ a i ( R , r ) . (5)The electronic energies ǫ a i ( R ) depend parametrically on R and obey the non-crossing rule [16]. The particular choiceof basis functions defined by Eq. (5) is called the adia-batic representation and is denoted by the superscript a.Inserting the expansion (4) into Eq. (3) and using the or-thonormality of the electronic states yields a set of cou-pled differential equations for the motion of the nuclei: (cid:18) − µ ∇ R − µ F a · ∇ R − µ G a + V a (cid:19) χ a = E χ a , (6)with F a ij ( R ) = (cid:10) φ a i |∇ R | φ a j (cid:11) , (7) G a ij ( R ) = (cid:10) φ a i (cid:12)(cid:12) ∇ R (cid:12)(cid:12) φ a j (cid:11) , (8) V a ij ( R ) = (cid:10) φ a i (cid:12)(cid:12) H el (cid:12)(cid:12) φ a j (cid:11) . (9)By definition, the matrix V a is diagonal and containsthe adiabatic potential energy curves ǫ a i ( R ). The matrixvector F a and the matrix G a both have off-diagonal el-ements, which are the non-adiabatic couplings betweendifferent electronic states. In addition, G a also has diag-onal elements commonly referred to as adiabatic energycorrections. If the electronic wave functions are taken tobe real, the chain rule for derivatives can be used to showthat F a is antisymmetric, i.e., F a ij ( R ) = − F a ji ( R ) . (10)The electronic states approximately correlate withthe separated atomic limits and therefore each term inEq. (4) defines a scattering channel. To determine theoutcome of any collision process, we need to solve thecoupled equations out to the asymptotic region wherethese limits are reached. However, it has been known fora long time that a product of the form χ a i ( R ) φ a i ( R , r )does not properly describe the translational motion ofthe electrons along with the nuclei [17, 18, 19]. Thiscan lead to difficulties of both formal and practical na-ture. One of the most prominent is the appearance ofnon-vanishing couplings at large or even infinite internu-clear distances [18, 19, 20]. However, it is demonstratedin Sec. V B that the influence of these couplings on thepresent results is expected to be small.The non-adiabatic coupling terms in Eq. (6) often varyrapidly with R and may cause problems in the numeri-cal solution procedure. Furthermore, the presence of F a leads to differential equations which contain both firstand second order derivatives. It is therefore desirable totransform Eq. (6) into an electronic basis where some orall of the non-adiabatic couplings disappear. FollowingMead and Truhlar [21], we shall refer to a representationin which all three components of F a vanish as strictlydiabatic. From this point on, it is assumed that the to-tal wave function can be represented in a finite set of M adiabatic basis functions. The transformation may thenbe written in terms of an orthogonal matrix T ij ( R ) [22]as φ d i ( R , r ) = M X j =1 φ a j ( R , r ) T T ij ( R ) . (11)In order for the total wave function to be preserved, thenuclear states must transform according to χ d i ( R ) = M X j =1 χ a j ( R ) T T ij ( R ) . (12)Using the relation (12) to replace the adiabatic states inEq. (6) yields after some manipulation (cid:18) − µ ∇ R − µ F d · ∇ R − µ G d + V d (cid:19) χ d = E χ d , (13)with F d = T T ( ∇ R + F a ) T , (14) G d = T T (cid:0) ∇ R + 2 F a · ∇ R + G a (cid:1) T , (15) V d = T T V a T . (16)The first derivative coupling F d expressed in the newbasis is seen to vanish provided that T is a solution tothe equation ( ∇ R + F a ) T = . (17)In the special case where all electronic states have thesame angular momenta projection onto the molecularaxis, only the radial component of F a can be non-zero [21]. Equation (17) then reduces to (cid:20) dd R + τ ( R ) (cid:21) T ( R ) = , (18)where τ ij ( R ) = (cid:28) φ a i (cid:12)(cid:12)(cid:12)(cid:12) ∂∂R (cid:12)(cid:12)(cid:12)(cid:12) φ a j (cid:29) . (19) Assuming that the first derivative couplings between theadiabatic states with i ≤ M and those with i > M arezero, also G d can be shown to vanish [22]. The coupledequations then take the form (cid:18) − µ ∇ R + V d (cid:19) χ d = E χ d . (20)In the strictly diabatic representation, the non-adiabatic couplings transform into off-diagonal elementsof the potential matrix V d (electronic couplings) and asa consequence the diagonal elements of this matrix (di-abatic potentials) no longer obey the non-crossing rule.As will be illustrated in Sec. III C, the strictly diabaticrepresentation may show little resemblance to the simpleand intuitive curve crossing picture often referred to asdiabatic. B. Asymptotic boundary conditions and theintegral cross section
The interpretation of χ d i ( R ) is usually simplified byrequiring that lim R →∞ T ( R ) = , (21)which makes the adiabatic and the diabatic representa-tions coincide in the region where the scattering wavefunction is evaluated. Thus, in any of the representa-tions, the physical situation is that of an incoming planewave in some entrance channel j , and outgoing spheri-cal waves in all channels which are energetically allowed.The direction of the incoming plane wave can be chosenalong the space fixed z-axis. If the potential falls off rea-sonably fast [23], the appropriate boundary conditionsare given by χ d i ( R ) ∼ R →∞ δ ij e ik j z + f ij ( E, θ ) e ik i R R , (22)where f ij ( E, θ ) is the scattering amplitude for enteringthe collision in channel j and leaving it in channel i . Sincethe potential matrix is spherically symmetric, the scat-tering amplitude is independent of the azimuthal angle ϕ .The asymptotic wave number k i is defined as k i = (cid:2) µ (cid:0) E − E th i (cid:1)(cid:3) / , (23)where E th i = lim R →∞ V d ii ( R ) is the corresponding thresh-old energy. In the case of a Coulomb potential, the ex-ponents in Eq. (22) have to be modified with logarithmicphase factors which arise due to the long range nature ofthis interaction [23].Evaluating the probability fluxes associated with thetwo terms in Eq. (22) leads to the familiar multichannelexpression for the integral cross section: σ ij ( E ) = 2 πk i k j Z π | f ij ( E, θ ) | sin θ d θ . (24)If the adiabatic states are connected by non-vanishingcouplings at infinity, neither the boundary condition (21)nor the scattering formalism developed below will be ap-propriate. In this case, the couplings can be manuallycut at some large but finite internuclear distance. Thisprocedure is usually justified if the collision energy is lowenough, but its validity should be considered from caseto case. C. The scattering formalism
The conventional approach to solving Eq. (20) is toexpress χ d ( R ) in terms of a partial wave expansion, χ d i ( R ) = 1 R ∞ X l =0 A l u i,l ( R ) P l (cos θ ) , (25)where P l (cos θ ) are the well known Legendre polynomi-als [24] and A l are constants to be chosen so that theboundary conditions (22) are fulfilled. The radial wavefunctions u i,l ( R ) can be shown to satisfy the equations (cid:20) − µ d d R + l ( l + 1)2 µR (cid:21) u i,l ( R ) + M X j =1 V d ij u j,l ( R )= Eu i,l ( R ) , (26)subject to the physical boundary conditions u i,l (0) = 0.If the corresponding derivatives u ′ i,l (0) are chosen appro-priately, these equations have M linearly independentsolutions which can be combined into a square matrix e u .It follows that (cid:18) d d R + Q l (cid:19) e u l = , (27)where Q l = 2 µ (cid:0) E − V d (cid:1) − l ( l + 1) R . (28)It is now assumed that the potential matrix is diagonalbeyond some radius ρ and here has elements correspond-ing to either pure Coulomb or pure short-range interac-tions. In this region, Eq. (27) is analytically solvable andthe radial wave functions may be expressed in terms ofthe incoming and outgoing wave solutions, α ij,l ( R ) and β ij,l ( R ), as defined in the Appendix. These are essen-tially the same as adopted by Johnson [25] but have beenextended to cover the Coulomb case. Thus, e u l ( R ) = R ≥ ρ α l ( R ) − β l ( R ) S l , (29)which defines the scattering matrix S l . Inserting this ex-pression into the partial wave expansion (25) and com-paring with the boundary conditions (22) leads to thewell known formula f ij ( E, θ ) = 12 i ( k i k j ) / ∞ X l =0 (2 l + 1) (cid:0) S oo ij,l − δ ij (cid:1) P l (cos θ ) , (30) for the scattering amplitude in terms of the open-openor physical partition of S l . The integral cross section forscattering from channel j to channel i is readily evaluatedto be σ ij ( E ) = ∞ X l =0 σ ij,l ( E ) , (31)where σ ij,l ( E ) = πk j (2 l + 1) (cid:12)(cid:12) S oo ij,l − δ ij (cid:12)(cid:12) . (32)Instead of working with S l directly, it is convenient tointroduce the reactance matrix K l [26] which is definedby e u l ( R ) = R ≥ ρ a l ( R ) − b l ( R ) K l , (33)where a ij,l ( R ) and b ij,l ( R ) are regular and irregular solu-tions to Eq. (27) explicitly given in the Appendix. Com-bining the definitions above together with those from theAppendix leads to the relationship [25] S oo l = (1 + i K oo l ) − (1 − i K oo l ) , (34)between the open-open partitions of S l and K l . This ex-pression resembles the Cayley transform and, since S oo l isboth symmetric and unitary, implies that K oo l is real [23].Finally, we introduce the logarithmic derivative y l = e u ′ l e u − l , (35)which transforms the radial Schr¨odinger equation intothe numerically more stable matrix Riccati equation [27,28] y ′ l + Q l + y l = , (36)while the boundary condition at the origin becomes thatof a diagonal matrix with infinite elements. This is theset of coupled equations that we ultimately have to solve.The relationship between y l and K l is easily seen to be K l = [ y l ( ρ ) b l ( ρ ) − b ′ l ( ρ )] − × [ y l ( ρ ) a l ( ρ ) − a ′ l ( ρ )] . (37) D. Modifications in the case of identical nuclei
Up to this point in the discussion complications thatmight arise if the nuclei are identical have been ignored.These may be considered in two steps. Firstly, since thecharges of the nuclei that enter H el are equal, the elec-tronic states will be either symmetric (gerade) or anti-symmetric (ungerade) under inversion of r . If the nucleiare labeled A and B, and α, β, γ, ... represent groups ofelectrons in some definite quantum states, each asymp-totic channel will be a linear combination of two config-urations, say A( α ) + B( β ) and A( β ) + B( α ). Such linearcombinations can not describe the situation where a nu-cleus is known to carry a certain number of electrons in acertain quantum state (unless α = β ). Localization of theelectron cloud can be achieved by forming proper linearcombinations of the gerade and ungerade solutions. The g / u symmetry is explicitly indicated as a superscript andthe channel index i is allowed to include all other quan-tum numbers. The appropriate scattering amplitudes arethen given by [23] f di ij ( E, θ ) = 12 (cid:2) f gij ( E, θ ) + f uij ( E, θ ) (cid:3) (38)and f ex ij ( E, θ ) = 12 (cid:2) f gij ( E, θ ) − f uij ( E, θ ) (cid:3) , (39)corresponding to the direct and exchange reactionsA( α ) + B( β ) → A( γ ) + B( δ ) (40)and A( α ) + B( β ) → A( δ ) + B( γ ) . (41)Here, the amplitudes f gij and f uij are given by the usualformula (30) applied to the gerade and ungerade mani-folds separately.A second complication arises from the fact that alsothe masses of the nuclei are equal. Consequently, it isimpossible to distinguish direct scattering in the direc-tion θ from exchange scattering in the direction π − θ .This effect can be accounted for by adding the direct andexchange amplitudes coherently. The fully symmetrizedscattering amplitudes are [29] f ± ij ( E, θ ) = f di ij ( E, θ ) ± f ex ij ( E, π − θ ) , (42)where the sign depends on whether the spatial part ofthe wave function should be symmetric (+) or antisym-metric (-) under exchange of the nuclei. To illustrate thispoint we consider a hydrogen quasi molecule where thetwo protons are known to be in a singlet spin state. Ac-cording to the Pauli principle, the total wave function(spin part included) must be antisymmetric under ex-change of the nuclei. Since the singlet spin state itself isantisymmetric under such an operation, the spatial partof the wave function is forced to be symmetric. Thus, f + ij is the appropriate choice of amplitude.The integral cross section is obtained by taking themodulus square of Eq. (42) and integrating the resultover the unit sphere: σ + ij ( E ) = X l even σ gij,l ( E ) + X l odd σ uij,l ( E ) (43)and σ − ij ( E ) = X l odd σ gij,l ( E ) + X l even σ uij,l ( E ) , (44) where σ gij,l and σ uij,l are given by Eq. (32) applied to thegerade and ungerade manifolds separately. The subscripteven (odd) indicates that only even (odd) partial wavesshould be summed over. In the particular case of the hy-drogen quasi molecule, the total spin state of the nucleiis either a singlet (para hydrogen) or a triplet (ortho hy-drogen). Taking into account the degeneracy factors (1and 3) and the exchange symmetry of these spin states,the symmetrized integral cross section for an ensemble ofparticles with randomly oriented spins is σ sym ij ( E ) = 14 X l even σ gij,l ( E ) + 34 X l odd σ gij,l ( E )+ 34 X l even σ uij,l ( E ) + 14 X l odd σ uij,l ( E ) . (45)This can be compared with the expression σ dist ij ( E ) = 12 ∞ X l =0 σ gij,l ( E ) + 12 ∞ X l =0 σ uij,l ( E ) , (46)obtained by treating the nuclei as distinguishable, andthus adding the direct and exchange amplitudes incoher-ently.In deriving Eqs. (45) and (46) we have assumed that allfinal states exist in both gerade and ungerade versions.This is not always the case. Of particular interest to thepresent study are processes in which the final electronicstate corresponds to γ = δ and exists only in a geradeversion. In this case, Eqs. (45) and (46) are still valid(setting σ uij,l = 0), although the arguments leading tothese formulas are somewhat modified. III. ELECTRONIC STRUCTURECALCULATIONS
The electronic states relevant to the present study arethose of H . The amount of data published for this sys-tem is of course substantial. Of particular relevance hereare the very accurate electronic structure calculations re-ported by Wolniewicz and co-workers [30, 31, 32, 33, 34,35, 36] using explicitly correlated basis functions and thecalculations by Detmer et al. [37] using a non-sphericalGaussian basis set. Reactions like the present one, how-ever, require the calculation of not only a large numberof excited states but also the associated non-adiabaticcouplings for a wide range of internuclear distances. Tothe best of our knowledge, a complete and accurate setof data for these quantities does not exist in the litera-ture. We have therefore performed ab initio electronicstructure calculations of all the relevant potential energycurves and non-adiabatic couplings. Where such data isavailable, our results have been compared with the bench-mark calculations of Refs. [30, 31, 32, 33, 34, 35, 36]. TABLE I: Specification of the spherical Gaussian basis setused in the present study.Type Exponent Type Exponent s a - p s p s p s p s d s d s d s d s d p d p d p f p f a Contraction of six primitive s functions with exponents (82.64,12.41, 2.824, 0.7977, 0.2581, 0.08989) and coefficients (0.002006,0.015343, 0.075579, 0.256875, 0.497368, 0.296133). A. Potential energy curves
We have calculated the adiabatic potential energycurves and the associated electronic wave functions corre-sponding to the seven lowest Σ + g states and the six low-est Σ + u states, here denoted as (1-7) Σ + g and (1-6) Σ + u ,respectively. All calculations have been performed at thefull configuration interaction (FCI) level using a modi-fied version of the Dalton 2.0 program package [38]. Themolecular orbitals have been obtained in a (11 s ,8 p ,7 d ,2 f )spherical Gaussian basis set contracted to [9 s ,8 p ,7 d ,2 f ].The compact basis functions of this set have been takenfrom the augmented correlation-consistent polarized va-lence quadruple zeta (aug-cc-pVQZ) basis set of Dunningand co-workers [39, 40], while the more diffuse ones havebeen approximately optimized to obtain a good repre-sentation of the excited states. The basis set is given inTable I.Figure 1 shows the (2-7) Σ + g and (1-6) Σ + u potentialenergy curves in the range 0.5 to 50 a . As is well known,only one bound state of the negative hydrogen ion exists;the S e ground state [41]. Accordingly, there is one ger-ade and one ungerade electronic state which asymptoti-cally correlate with the H + + H − ion-pair limit. As theinternuclear distance is decreased from infinity, the ion-pair configuration is successively carried on through a se-ries of avoided crossings involving states correlated withthe H(1) + H( n ) covalent limits. In the neighborhood ofthese crossings, the first derivative radial couplings (19)are significant and are able to cause transitions betweendifferent adiabatic states.Beginning with the gerade symmetry, the (1-6) Σ + g manifold consists of states correlated with the n = 1 , , n )there are n of these states in total. The 7 Σ + g state hasthe ion-pair configuration from 36 a up to approximately280 a . Here, the adiabatic wave function changes char- -0.75-0.70-0.65-0.60-0.55 0 10 20 30 40 50 P o t en t i a l ene r g y ( un i t s o f E h ) Internuclear distance (units of a ) H(1) + H(2)H(1) + H(3)H + + H - (a)-0.75-0.70-0.65-0.60-0.55 0 10 20 30 40 50 P o t en t i a l ene r g y ( un i t s o f E h ) Internuclear distance (units of a ) H(1) + H(2)H(1) + H(3)H + + H - (b) FIG. 1: Adiabatic potential energy curves for (a) the (2-7) Σ + g states and (b) the (1-6) Σ + u states of H . The separatedatomic limits are indicated at the far right of each subfig-ure. For a discussion of the particular ion-pair limit, see themain text. acter into the n = 4 covalent configuration while the ion-pair configuration is carried on to another higher excitedstate. At such large distances, however, the motion ofthe nuclei is highly diabatic and the change of characterof the adiabatic wave functions will not induce any sig-nificant ionic-covalent transitions [4, 42, 43]. With theseconsiderations in mind, we will refer to the 7 Σ + g stateas that which is correlated with the ion-pair limit.In the ungerade symmetry, the (1-5) Σ + u manifold con-sists of states correlated with the n = 2 , Σ + u state, with the same considerations asabove, can be said to correlate with the ion-pair limit.Because of symmetry restrictions, there can be no sin-glet ungerade state dissociating into the n = 1 covalentlimit. Note that while the gerade and ungerade statesbehave very differently at small internuclear distances,these differences vanish quickly as the separation of thenuclei is increased. In fact, at 10 a the energy differencefor all excited states is less than 9 × − E h (Hartree)and at 30 a it is less than 5 × − E h .Figure 1 illustrates the two avoided crossings near 12and 36 a , which are related to a change from ionic to n = 2 and 3 covalent character (or vice versa) of theadiabatic wave functions, respectively. In what followswe will refer to them simply as the n = 2 and 3 curvecrossings. Figure 1 also shows a rich pattern of avoidedcrossings taking place at smaller internuclear distances.This region, however, is to a large extent shielded bythe centrifugal term in Eq. (26) and the influence of theresulting radial couplings upon the total neutralizationcross section is rather small (see Sec. V B). The higherexcited Σ + g,u states that are not considered in the presentstudy are only accessible to the system through avoidedcrossings at small internuclear distances [37]. The errorintroduced by excluding these states is therefore expectedto be of less significance. In principle, also rotationalcouplings to states of Π symmetry should be considered.However, due to the presence of strong radial couplingsin the vicinity of the n = 2 and 3 curve crossings, anddue to the low collision energies considered in the presentstudy, these are assumed to be negligible [9, 44].To evaluate the quality of our calculated potentialenergy curves, we have compared them with those re-ported by Wolniewicz and co-workers. Their calcula-tions cover the 1 Σ + g state out to 12 a [32, 34], the (2-6) Σ + g states out to 20 a [33] and the (1-6) Σ + u statesout to 150 a [36]. The potential energy curve for the4 Σ + g state has been recalculated and extended out to80 a [35]. With the energy scale used in Fig. 1, thesepotentials would entirely overlap with ours. A more de-tailed comparison shows that the largest energy differ-ence is obtained for the inner part ( R ≤ . ) of the1 Σ + g ground state potential, which is of little impor-tance to the present study anyway. Here, our computedvalues are approximately 5 . × − to 1 . × − E h higher in energy compared to those in Ref. [34]. Forall of the excited states, our calculated energies are ap-proximately 1 . × − to 5 . × − E h higher than inRefs. [33, 35, 36]. It is of particular interest to estimatethe quality of the (4-7) Σ + g and (3-6) Σ + u states near thecritical n = 3 curve crossing. Concerning the unger-ade states, the covalent parts of our potentials are within1 . × − E h and the ionic parts within 3 . × − E h tothose reported in Ref. [36]. Due to the near degeneracy ofthe gerade and ungerade states, we expect this accuracyto be common to both inversion symmetries. The qual-ity of the adiabatic states and in particular their energysplitting at the n = 3 curve crossing will be discussedfurther in the next section. B. Radial couplings
We have calculated the first derivative radial cou-plings (19) between all of the adiabatic states consideredin the previous section. The accurate evaluation of thesequantities is not trivial. The magnitude and shape of theradial couplings are usually very sensitive to the qualityof the electronic wave functions and it is essential that not only the individual states are good in a variationalsense, but also that the relative energies of these statesare well represented. In regions where two or more adia-batic states become nearly degenerate, further complex-ity is added by the fact that a minor change in the atomicbasis can cause the potential energy curves to artificiallypseudo cross, giving rise to dramatic effects in the radialcouplings. These effects do not change the outcome ofthe dynamics but can obscure both interpretation andcomparison with other results.To evaluate the radial couplings we have implementeda three point version of the finite difference method de-scribed in Ref. [45]. This method provides a systematicway to converge the calculations toward the exact resultin the particular type of basis and wave function beingconsidered. We have examined carefully how the cal-culated radial couplings depend on the derivative steplength ∆ R , as well as the numerical accuracies, δ MO and δ CI , with which the molecular orbitals and the CI wavefunctions are obtained. Stable and converged results wereobserved when ∆ R = 5 × − a , δ MO = 10 − E h and δ CI = 10 − E h . The antisymmetry condition (10) canbe used as a simple consistency check of the results. Inthe present case, most of the calculated coupling elementssatisfied this condition to at least four significant digits.For practical reasons, it is still desirable to have a cou-pling matrix that is fully antisymmetric. This has beenaccomplished by taking each element to be the geomet-rical mean of τ ij and − τ ji .The radial couplings are too many to discuss in fulldetail and so we will focus on those relevant for the im-portant n = 2 and 3 curve crossings. The notation ( i, j ) g is used to label the radial coupling between the i -th and j -th Σ + g state (and similarly for the ungerade states).Due to symmetry, gerade and ungerade states are notcoupled to each other.In Fig. 2 we show the radial couplings among the (2-4) Σ + g and (1-3) Σ + u states in the neighborhood of the n = 2 curve crossing. As can be expected from theappearance of the potential energy curves, each of the(2 , g and (1 , u couplings show a rather broad peaknear the crossing point at 12 a resulting from the ex-change of ionic and n = 2 covalent character betweenthe adiabatic wave functions. Similar shapes and mag-nitudes, although displaced to larger internuclear dis-tances, are also exhibited by the (2 , g and (1 , u cou-plings, which connect the adiabatic states tending to the n = 2 covalent limit. The (3 , g and (2 , u couplingsare, on the other hand, much less pronounced in thecrossing region. For comparison, also plotted in Fig. 2 arethe results obtained by Wolniewicz and Dressler [30, 33],which include the radial couplings among all of the (2-4) Σ + g states and between the 1 and 2 Σ + u states. Theoverall agreement with the present results is very good,in particular in the region up to approximately 16 a .For larger internuclear distances some deviations can beobserved. These are most likely due to a slight unbal-ance in our description of the nearly degenerate n = 2 -0.25-0.20-0.15-0.10-0.050.000.050.100.15 4 6 8 10 12 14 16 18 20 C oup li ng s t r eng t h ( un i t s o f a - ) Internuclear distance (units of a )(2,3) g (2,4) g (3,4) g (a)-0.25-0.20-0.15-0.10-0.050.000.050.100.15 4 6 8 10 12 14 16 18 20 C oup li ng s t r eng t h ( un i t s o f a - ) Internuclear distance (units of a )(1,2) u (1,3) u (2,3) u (b) FIG. 2: Radial couplings among (a) the (2-4) Σ + g states and(b) the (1-3) Σ + u states in the neighborhood of the n = 2curve crossing. The notation ( i, j ) g label the radial couplingbetween the i -th and j -th Σ + g state (and similarly for theungerade states). The present results (solid lines) are com-pared with the results obtained by Wolniewicz and Dressler(crosses) [30, 33]. adiabatic states. However, as mentioned at the very be-ginning of this section, the physical relevance of thesetypes of effects is expected to be small.Figure 3 illustrates the radial couplings among the adi-abatic states associated with the n = 3 curve crossing.Here, the gerade and ungerade states are almost com-pletely degenerate and so only the results for the for-mer of these, i.e., the (4-7) Σ + g manifold of states, areshown. All the couplings appear as peak like structureswith their maxima located in the region 35.5 to 37 . a .The most prominent of these is the (4 , g coupling whichreaches its maximum value at 35 . a . This couplingarises from the exchange of ionic and n = 3 covalentcharacter between the 4 and 7 Σ + g adiabatic states. Wenote that the location of its maximum agrees well withthe crossing point of the pure ionic and covalent energiesat 36 . a obtained from the ground state energy [46]and the polarizability [47] of the H − ion. The calcula-tions of Wolniewicz and Dressler do not extend out to -1.2-1.0-0.8-0.6-0.4-0.20.00.20.4 28 30 32 34 36 38 40 42 44 C oup li ng s t r eng t h ( un i t s o f a - ) Internuclear distance (units of a )(4,5) g (4,6) g (4,7) g (5,6) g (5,7) g (6,7) g FIG. 3: Radial couplings among the (4-7) Σ + g states in theneighborhood of the n = 3 curve crossing. The notation ( i, j ) g label the radial coupling between the i -th and j -th Σ + g state. the n = 3 curve crossing and thus no direct comparisonof the radial couplings can be made. However, in theungerade symmetry we can compare the energy splittingof the 3 and 6 Σ + u states at 35 . a . Here, we obtainthe value 3 . × − E h which is in very good agreementwith the value 3 . × − E h calculated from the resultsof Ref. [36]. C. Adiabatic to diabatic transformation
In order to obtain the adiabatic to diabatic transforma-tion (ATDT) matrix, we have numerically solved Eq. (18)using the radial couplings described in the preceding sec-tion as input data. The solution has been obtained with amatrix version of the Runge-Kutta Fehlberg method [48].The boundary condition T = has been imposed at R = 50 a and the integration has been performed in-wards. As hinted above, several of the non-adiabaticcouplings are non-vanishing as R tends to infinity. Inpractice, our choice of boundary condition correspondsto setting these couplings to zero beyond 50 a . Sincethis kind of assumption has to be made in the calcula-tion of the scattering matrix anyway, no additional erroris introduced due to the diabatization procedure itself.In the present case, the transformation to a strictly dia-batic basis is merely for computational reasons and is notdiscussed further. It is worth noting, however, that thediabatic potential energy curves we obtain here are farfrom the intuitive ones conventionally used in connectionwith the H system, see Figure 4. IV. NUMERICAL PROCEDURES
The logarithmic derivative (35) has been calculated bynumerically solving the matrix Riccati equation (36) from -1.0-0.9-0.8-0.7-0.6-0.5 0 10 20 30 40 50 P o t en t i a l ene r g y ( un i t s o f E h ) Internuclear distance (units of a ) H(1) + H(2)H(1) + H(3)H + + H - FIG. 4: Strictly diabatic potential energy curves obtainedfrom a transformation of the (1-7) Σ + g adiabatic states. Thediabatic potential energy curve associated with the H(1) +H(1) limit is not shown in the figure. .
54 a out to 50 a . For this purpose we have used analgorithm developed by Johnson [27, 28] with the gridsize set to 5 × − a . Knowledge of the logarithmicderivative in the final integration point has allowed usto calculate the partial cross sections (32) for scatteringwithin the gerade and ungerade inversion symmetries.The fully symmetrized integral cross sections have beenobtained by combining the partial cross sections accord-ing to Eq. (45). To determine at which point the seriesin Eq. (45) could be truncated, a simple convergence cri-teria was introduced. This was set to terminate the sum-mation if the ratios of the partial cross sections and theaccumulated integral cross sections remained less than10 − for 25 terms in succession. Over the energy rangeconsidered here, this led to the inclusion of approximately250 to 3500 partial waves. V. RESULTS AND DISCUSSIONA. State dependent cross sections
The calculated cross sections for scattering into theH(1)+H( n ), n = 1 , , n = 3 process dominates the others overalmost the entire energy range, which is reasonable con-sidering the favorable location of the n = 3 curve crossingat around 36 a . In contrast, the n = 2 transitions oc-cur much further in and are more easily suppressed bythe centrifugal barrier. Only for collision energies above5 eV, where the n = 2 cross section exhibits a mini-mum, is this process able to compete with neutralizationinto the n = 3 final states. This phenomena was firstnoted by Bates and Lewis [4] and has been confirmed, atleast on a qualitative level, by almost every study sincethen [7, 8, 9]. Neutralization into the n = 1 final state -19 -18 -17 -16 -15 -14 -13 -12 -11 C r o ss s e c t i on ( c m ) Collision energy (eV)Present: n = 1n = 2n = 3Fussen and Kubach: n = 2n = 3
FIG. 5: Calculated cross sections for scattering into theH(1) + H( n ), n = 1 , , is insignificant at all energies and is included only forcompleteness. As the collision energy approaches zero,all cross sections gain the characteristic E − dependencethat can be expected for reactions which are governed bythe Coulomb interaction. Here, the n = 3 cross section isapproximately a factor of 50 larger than the n = 2 crosssection.Some weaker oscillations can be observed over a largepart of the energy range. The reason for this structureis not completely clear, but the most likely explanationis in terms of St¨uckelberg oscillations, i.e., in terms ofquantum interference arising because there are severalcompeting routes through the potential landscape lead-ing to the same asymptotic limit.The close-coupling results of Fussen and Kubach [9]have also been plotted in Fig. 5 for comparison. Theagreement between their n = 3 cross section and thepresent one is good. The difference is about 20 percentat 1 eV and increases slightly with the energy. Com-paring the n = 2 cross sections, the two calculations areseen to agree well above the minimum at 5 eV, whereasat lower energies the differences are more pronounced. Arationale for this discrepancy may be found by consider-ing how the inversion symmetry of the electronic statesis approached. For this purpose we show in Fig. 6 a com-parison of our fully symmetrized cross sections (same asin Fig. 5) together with the pure gerade and ungeradecross sections calculated in the conventional way. Wefirst note that the n = 3 cross section is fairly insensi-tive to the choice of inversion symmetry over the wholeenergy range. Consequently, the effect of averaging overthese symmetries is expected to be small. A similar ob-servation can be made regarding the n = 2 cross sectionabove 5 eV. However, at lower energies the gerade andungerade results differ strongly. This indicates not onlythat our way of combining the gerade and ungerade cross0 -19 -18 -17 -16 -15 -14 -13 -12 -11 C r o ss s e c t i on ( c m ) Collision energy (eV) n = 1n = 2n = 3SymmetrizedGeradeUngerade
FIG. 6: Fully symmetrized vs. pure gerade and ungeradecross sections for scattering into the H(1) + H( n ), n = 1 , , sections will affect the final result, but also that any twoapproaches that treat the inversion symmetry differentlymay lead to different outcomes. In the one electron modelof Fussen and Kubach, this inversion symmetry is actu-ally broken and as a consequence their results are likelyto differ from ours. It should be emphasized though, thatat the energies where this happens the contribution fromthe n = 2 channels to the overall reaction is small. Itshould also be noted that the effects of treating the nu-clei as indistinguishable are very small, and it is equallygood to use the simpler formula (46) in favor of (45) toweight the gerade and ungerade partial cross sections. B. Total neutralization cross section
In Fig. 7 we show the calculated total neutralizationcross section, i.e., the sum of the n = 1 , , E − behavior discussed above.The cross section then continues to fall off until it reachesa broad minimum at around 20 eV. This minimum obvi-ously results from the nearly constant n = 3 cross sectionin combination with the rapid increase in the formationof n = 2 final products.To examine the influence of the non-vanishing asymp-totic couplings upon the cross section, the point at whichthe logarithmic derivative is evaluated has been variedover the range 50 to 80 a . The boundary condition onthe ATDT matrix has been varied accordingly. At thecollision energies 0.001, 1 and 100 eV, the relative varia-tion in the cross section was seen to be 4 × − , 7 × − and 1 × − , respectively. We conclude that the effects ofthe non-vanishing asymptotic couplings over the energyrange considered here are negligible. This is consistentwith the conclusions drawn by Borondo et al. [43]. Additional tests have been performed in order to esti-mate the relevance of the radial couplings at small inter-nuclear distances (as discussed in Sec. III A). This hasbeen done by setting all radial couplings up to a certaindistance R to zero before the diabatization procedure iscarried out. Choosing R = 6 a reduces the total neu-tralization cross section at 0.001, 1 and 100 eV collisionenergy by 8.6, 9.9 and 5.9 percent, respectively. Thesevariations in the cross section can be considered small(although not insignificant) and may be taken as an up-per limit for the influence of the higher lying Σ + g,u adia-batic states not included in the present study, which arecoupled to the lower lying states mainly at small inter-nuclear distances. A more detailed study of these effectsis currently being pursued.In Fig. 7 we also show the results of earlier studiesrelevant for the mutual neutralization of H + and H − atlow collision energies. Up to the present date, no mea-surements of this cross section have been reported forcollision energies below 0.15 eV. In the region 0.15 to3 eV, the only experimental results that are available arethose of Moseley et al. [10]. A comparison shows thatthese are a factor of two to three larger than the presentdata. At energies above a few eV, experiments have beenconducted by Szucs et al. [11] and Peart and Hayton [12],among others. In the region 3 to 10 eV, their results pointtowards a cross section that is slightly lower than ours,which is most evident for the first few data points of Peartand Hayton, whereas for collision energies above 10 eVall three studies agree favorably. Comparing with earliertheoretical studies, our results match within 30 percentthe close-coupling calculations of Fussen and Kubach [9]and within 40 percent the Landau-Zener calculations ofBates and Lewis [4]. Our results also more or less over-lap with the Landau-Zener results obtained by Eerden etal. [8].It may be worthwhile to review the accumulated dataon the mutual neutralization of H + and H − at low col-lision energies (below a few eV). There are now severalcalculations of this cross section that are mutually con-sistent within a few tens of percent. These employ suchdiverse methods as semi-classical Landau-Zener theory aswell as model Hamiltonian and ab initio close-couplingschemes. In the upper part of the energy region, thetheoretical results are in more or less agreement withseveral of the available experiments. One of these mighteven suggest a cross section that is slightly lower thanpredicted by theory. It is a fact, however, that the onlymeasurement that has been published for collision ener-gies below 3 eV shows a cross section that is consider-ably higher than others reported. To clarify this situ-ation, new experimental efforts are obviously required.Recently such measurements have been undertaken bythe group of Urbain [49] at the Universit´e Catholique deLouvain, Belgium, using a merged beam apparatus in theenergy range 0.001 to 0.2 eV. Although uncertainties ofthe order of 50 percent still exist, preliminary results fromthis experiment seem to support the neutralization cross1 -14 -13 -12 -11 C r o ss s e c t i on ( c m ) Collision energy (eV) PresentBates and Lewis (theor.)Fussen and Kubach (theor.)Eerden et al. (theor.)Moseley et al. (exp.)Szucs et al. (exp.)Peart and Hayton (exp.) 10 -14 -13
1 10 100 C r o ss s e c t i on ( c m ) Collision energy (eV)(a) (b)
FIG. 7: Total cross section for the mutual neutralization of H + and H − . The present results are compared with the measure-ments by Moseley et al. [10], Szucs et al. [11], and Peart and Hayton [12], as well as the calculations by Bates and Lewis [4],Fussen and Kubach [9], and Eerden et al. [8]. (a) Cross sections plotted over the full energy range 0.001 to 100 eV. (b)Enlargement of the energy range 1 to 100 eV. section obtained from the various theoretical treatments. C. Rate coefficient
At low collision energies, the total neutralization crosssection may be parametrized in terms of the relative ve-locity v as [10] σ ( v ) = Av − + Bv − + C + Dv . (47)To determine the unknown coefficients of this expression,we have performed a least square fit to our calculatedcross section in the region 0.001 to 4 eV. The result is A = 4 . × − cm s − ,B = − . × − cm s − ,C = 1 . × − cm ,D = − . × − cm s , corresponding to a maximum error of six percent in theparametrized cross section. By integrating Eq. (47) overa Maxwellian velocity distribution, the associated ratecoefficient α ( T ) is obtained: α ( T ) = A (cid:18) µπkT (cid:19) / + B +2 C (cid:18) kTπµ (cid:19) / + 3 DkTµ , (48)where T is the ion temperature and k the Boltzmann con-stant. We estimate the above results to be valid over thetemperature range 10 to 10,000 K. For comparison, therate coefficient we obtain at 300 K is 1 . × − cm s − . VI. SUMMARY
In this paper we have considered the mutual neutral-ization of H + +H − into H(1)+H( n ), n = 1 , , n = 3 final states. Furthermore,a proper weighting of the gerade and ungerade cross sec-tions turns out to be important only for the n = 2 forma-tion at low collision energies, where these cross sectionsanyway are small. In the collision energy region below afew eV, our total neutralization cross section is in goodagreement with previous theoretical studies [4, 8, 9], butis a factor of two to three lower than that measured byMoseley et al. [10]. Acknowledgments
This work was supported by the Swedish ResearchCouncil (VR). We thank X. Urbain and A. E. Orel foruseful discussions, and R. D. Thomas for reading and2commenting on the manuscript.
APPENDIX
Here, we provide analytical expressions for the solutionmatrices a l , b l , α l and β l which are used in the defini-tions of the reactance and scattering matrices of Sec. II C.Open and closed channels are distinguished depending onwhether the total energy is greater or less than the cor-responding threshold energy. The channels are labeledshort-range (covalent) if V d ii ( R ) = R ≥ ρ E th i (A.1)and Coulombic if V d ii ( R ) = R ≥ ρ E th i + Q Q R , (A.2)where Q and Q are the net charges situated on eachof the nuclei. The asymptotic wave number k i for openchannels has been defined in Eq. (23). For closed chan-nels, k i = (cid:2) µ (cid:0) E th i − E (cid:1)(cid:3) / . (A.3)It is further convenient to introduce the scaled coordinate ξ i = k i R and the Coulomb parameter η i = µQ Q k i . (A.4)We look for pairs of linearly independent solutions tothe radial Schr¨odinger equation (27) in the region R ≥ ρ .These are chosen such that a ii,l ( R ) and b ii,l ( R ) behave asregular and irregular solutions, and α ii,l ( R ) and β ii,l ( R )as incoming and outgoing wave solutions, respectively.For open short-range channels, two linearly independentsolutions are the Riccati-Bessel functions of the first andsecond kind, e S l ( ξ i ) and e C l ( ξ i ), respectively [24]. Weadopt the definitions [25] a ij,l = δ ij k − / i e S l , (A.5) b ij,l = δ ij k − / i e C l , (A.6) α ij,l = δ ij k − / i (cid:16) − i e S l + e C l (cid:17) , (A.7) β ij,l = δ ij k − / i (cid:16) i e S l + e C l (cid:17) . (A.8)For open Coulomb channels, the radial Schr¨odinger equa-tion is solved by the regular and irregular Coulomb func-tions, e F l ( η i , ξ i ) and e G l ( η i , ξ i ), respectively [24]. Accord-ingly, a ij,l = δ ij k − / i e F l , (A.9) b ij,l = δ ij k − / i e G l , (A.10) α ij,l = δ ij k − / i (cid:16) − i e F l + e G l (cid:17) , (A.11) β ij,l = δ ij k − / i (cid:16) i e F l + e G l (cid:17) . (A.12)Finally, in the case of closed short-range channels, pos-sible solutions are ξ i e i l ( ξ i ) and ξ i e k l ( ξ i ), where e i l ( ξ i ) and e k l ( ξ i ) are the modified spherical Bessel functions of thefirst and second kind, respectively [24] (first and thirdkind in the cited reference). Appropriate definitionsare [25] a ij,l = δ ij / π − / ξ i e i l , (A.13) b ij,l = δ ij / π − / ξ i e k l , (A.14) α ij,l = δ ij (2 π ) / ξ i k − i he i l + ( − l π − e k l i , (A.15) β ij,l = δ ij / π − / ξ i k − i e k l . (A.16)Closed Coulomb channels play no role in the presentstudy and are therefore not considered. [1] D. Galli and F. Palla, Astron. Astrophys. , 403(1998).[2] S. C. O. Glover, ApJ , 331 (2003).[3] S. C. Glover, D. W. Savin, and A.-K. Jappsen, ApJ ,553 (2006).[4] D. R. Bates and J. T. Lewis, Proc. Phys. Soc. A , 173(1955).[5] L. Landau, Phys. Z. Sowjetunion , 46 (1932).[6] C. Zener, Proc. R. Soc. Lond. A , 696 (1932). [7] R. E. Olson, J. R. Peterson, and J. Moseley, J. Chem.Phys. , 3391 (1970).[8] M. J. J. Eerden, M. C. M. van de Sanden, D. K. Otor-baev, and D. C. Schram, Phys. Rev. A , 3362 (1995).[9] D. Fussen and C. Kubach, J. Phys. B: At. Mol. Phys. ,L31 (1986).[10] J. Moseley, W. Aberth, and J. R. Peterson, Phys. Rev.Lett. , 435 (1970).[11] S. Szucs, M. Karemera, M. Terao, and F. Brouillard, J. Phys. B: At. Mol. Phys. , 1613 (1984).[12] B. Peart and D. A. Hayton, J. Phys. B: At. Mol. Opt.Phys. , 5109 (1992).[13] B. Peart, M. A. Bennett, and K. Dolder, J. Phys. B: At.Mol. Phys. , L439 (1985).[14] B. Peart, S. J. Foster, and K. Dolder, J. Phys. B: At.Mol. Opt. Phys. , 1035 (1989).[15] S. Geltman, Topics in Atomic Collision Theory (KriegerPublishing Co., Malabar, 1997).[16] J. von Neumann and E. Wigner, Phys. Z. , 467 (1929).[17] D. R. Bates and R. McCarroll, Proc. R. Soc. Lond. A , 175 (1958).[18] J. B. Delos, Rev. Mod. Phys. , 287 (1981).[19] A. K. Belyaev, D. Egorova, J. Grosser, and T. Menzel,Phys. Rev. A. , 052701 (2001).[20] D. W. Jepsen and J. O. Hirschfelder, J. Chem. Phys. ,1323 (1960).[21] C. A. Mead and D. G. Truhlar, J. Chem. Phys. , 6090(1982).[22] M. Baer, Phys. Rep. , 75 (2002).[23] N. F. Mott and H. S. W. Massey, The Theory of AtomicCollisions (Oxford University Press, London, 1965).[24] M. Abramowitz and I. A. Stegun,