Analytical Results for the Three-Body Radiative Attachment Rate Coefficient, with Application to the Positive Antihydrogen Ion H+

To overcome the numerical difficulties inherent in the Maxwell–Boltzmann integral of the velocity-weighted cross section that gives the radiative attachment rate coefficient α R A for producing the negative hydrogen ion H − or its antimatter equivalent, the positive antihydrogen ion H ¯ + , we found the analytic form for this integral. This procedure is useful for temperatures below 700 K, the region for which the production of H ¯ + has potential use as an intermediate stage in the cooling of antihydrogen to ultra-cold (sub-mK) temperatures for spectroscopic studies and probing the gravitational interaction of the anti-atom. Our results, utilizing a 50-term explicitly correlated exponential wave function, confirm our prior numerical results.


Introduction
The Antiproton Decelerator (AD) facility at CERN [1] has provided the foundation for a variety of experiments (e.g., [2][3][4]) over more than a decade. Small numbers of these anti-atoms are trapped by the ALPHA and ATRAP collaborations using specialized magnetic minimum neutral atom traps [5][6][7], with confinement times of many minutes being routine at ALPHA [8]. They have done spectroscopic [9] measurements for H in their quest to investigate possible violations of CPT symmetry, experimental limits on its charge [10], and preliminary limits on the gravitational interaction of the anti-atom [11].
Building on the latter idea, the GBAR collaboration [12][13][14] means to measure the gravitational attraction of matter versus antimatter using neutral H atoms, but cooling them sufficiently is difficult because of their neutrality. They intend to form the antihydrogen ion H + as an intermediate step because its net charge would allow for sympathetic cooling with a mixture of positively charged ions of ordinary matter such as Be + , and, after they are cooled, the extra positron would be stripped off prior to studies of the gravitational interaction of the anti-atom [12][13][14]. The authors of [15,16] calculated the cross section and rate coefficient for the radiative attachment of a second positron to create the H + ion, H (1s) + e + → H + 1s 2 1 S e +hω . (1) We first [15] used the effective range wave function of Ohmura and Ohmura [17] and then [16] a fully two-positron 200-term wave function [18] composed of explicitly correlated exponentials of the kind introduced by Thakkar and Smith [19]. These extend to temperatures lower than Bhatia's [20] results. Calculating the radiative attachment rate coefficient α RA for producing the negative hydrogen ion H − or its antimatter equivalent, the positive antihydrogen ion H + , requires the evaluation of a Maxwell-Boltzmann integral of the velocity-weighted cross section whose integrand is akin to a slightly-rounded Heaviside step function that is difficult to handle numerically, particularly for temperatures below 1000 K. Evaluating this integral analytically would, then, be ideal, perhaps using the analytical results for the underlying six-dimensional photoionization integral for the cross section itself given in Keating's master's thesis [21]. However, integrating squares of sums of the large variety of terms in that final cross section is a daunting task. This variety of terms arises in Keating's work as the Lth derivatives of the Laplace transform of the spherical Bessel function j 1 (kr), where the "(2)" on the summation sign in the last line indicates steps of two, p = α + γ or α + β of Equation (8), and k is the wave number. One might wonder, then, if one could back up to the final radial integral of the cross section that has a consistent analytic form r 3/2+h e −σr j 1 (kr) for the direct and cross terms. It indeed turns out to be possible to perform the Maxwell-Boltzmann integral of products of that analytic form first, and then integrate over each of the radial integrals in the product r 3/2+h e −σr j 1 (kr) R 3/2+j e −τr j 1 (kR), and finally sum over all such product terms and the terms of the explicitly correlated exponential wave function. We give a synopsis of how one finds the radiative attachment cross section, explicate the integrals one needs to calculate, and show how a fully analytical rate coefficient may be found. Tests of the new form confirm the numerical integrals of prior work.

The Radiative Attachment cross Section
Since this approach relies on finding the radiative attachment cross section from the photodetachment cross section via the principle of detailed balance, we give a short history. Photodetachment from the hydrogen ion H − , for instance, is known to be responsible for the opacity of the sun [22,23], garnering much attention in the 1940s-1980s [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] and more recently [43][44][45][46][47][48][49][50][51]. Ward, McDowell, and Humberston [52] described the parallel Ps − case as calculating an allowed dipole transition to the continuum of the two-electron (or two-positron) Hamiltonian Following Ghoshal and and Ho [50], we put the (length gauge 1 ) cross section for photodetachment (or photoionization), σ PI , laid out by Chandrasekhar [24,25] in atomic units, and obtain where α is the fine structure constant, a 0 is the Bohr radius, and the magnitude of the momentum of the detached electron or positron p =hk may be related to the energy ω (and the photon wavelength λ ) and the H − electron or positron affinity, I, by Since the cross section differences between velocity and length gauge formulations (due to the approximate nature of the two-positron wave functions used) are small, we will present only length gauge results in this work. In the matrix element for photodetachment in Equation (4), ψ f is a continuum state wave function for the outgoing positron represented by (the dipole term of) a plane wave multiplied by a hydrogen ground state wave function, Accurate two-positron initial-state wave functions may take the form ψ H (r 1 , r 2 , r 12 ) = 1 √ 2 1 −P 12 e −αr 1 −βr 2 −γr 12 ∑ l,m,n c lmn s l t 2m u n , whereP 12 is the permutation operator for the two identical positrons α ↔ β, with Hylleraas coordinates [53] given by s = r 1 + r 2 , t = r 1 − r 2 , and u = r 12 ≡ |r 1 − r 2 |. One may also express this as sums of powers of r 1 and r 2 instead of powers of s and t, via the binomial theorem. Often, the difficulty of finding the nonlinear parameters in the exponential is reduced by setting β = α and γ = 0. Alternatively, Thakkar and Smith [19] introduced a set of wave functions involving solely exponentials, with the nonlinear inter-electron (inter-positron) correlation parameter γ retained, where the parameters in the exponentials are generated in a quasi-random fashion, where x denotes the fractional part of x. The downside of having to find six nonlinear parameters that minimize the energy, rather than the single nonlinear parameter one varies in a many Hylleraas expansions, is sufficiently compensated for in that the wave function has a consistent form and is generally easier for evaluating integrals. For the fifty-term wave function we use, these parameters are [54]: A 1 = 0.2380, A 2 = 1.3240, B 1 = 0.9800, B 2 = 1.3290, G 1 = −0.0720, G 2 = 0.288, and η = 1 − 2.458 × 10 −7 . The quasi-random assignment of the 50 values for each of α k , β k , and γ k in Equation (10) means that we do not have to vary these 150 parameters directly. Because the optimization algorithm is not perfect, one must scale the wave function with η very slightly different from one so that it satisfies the virial theorem. The coefficients c k are found by diagonalizing the Hamiltonian matrix in order to minimize the ground state energy and then normalized. Using this wave function for a given triplet of α = α k , β = β k , and γ = γ k in the sum in Equation (9), the matrix element for ionizing either positron under the influence of the length dipole operator (z 1 + z 2 ) is the sum of four terms: × e −λr 2 e ikz 1 + e −λr 1 e ikz 2 where the first factor of two comes from I 11 = I 22 and so on, whose subscripts j refer to z j in the dipole operator and in the plane wave, respectively, and we have factored out the coefficient that is common to all terms. We keep λ, the magnitude of the charge of the hydrogen nucleus, in symbolic form rather than setting it to one so that we can take derivatives of it to represent powers of r j . The cross section for radiatively attaching a second positron to H (1s) to create the (1s 2 1 S e ) state of the H + ion, via the reaction in Equation (1), can be obtained from the principle of detailed balance (see, e.g., Landau and Lifshitz [55]) following the lead of Drake [56] and then Jacobs, Bhatia, and Temkin [57], who applied the principle of detailed balance to obtain the radiative attachment coefficient (for an electron) to form the (2p 2 3 P e ) metastable H − state from H (2s, 2p). For the (1s 2 1 S e ) case, we have [16], where g 1 /g 2 = 6/12 is the statistical weight ratio. Here, the photon momentum relative to the ion is given by p ω =hω/c = (k 2 + γ 2 0 )/2c, and p e is the positron momentum k. Note that c in atomic units is the inverse of the fine structure constant α.
To estimate formation rates of H + , it is helpful to calculate the positron attachment to H as a function of temperature rather than energy, as is common in astrophysical applications. This rate coefficient α RA is formed as the expectation value of vσ RA with the normalized Maxwell-Boltzmann distribution f (v) as, whose overall coefficient is where the temperature T is given in Kelvins K.

Evaluating Integrals
The algebra in each case is greatly reduced by making the replacement e −γr 12 r 12 (15) in each term, with the notation where Atoms 2020, 8, 13

of 13
For the cross-terms z 2 e ikz 1 = r 2 cos θ 2 e ikr 1 cos θ 1 , we need the conversion given in Ley-Koo and Bunge [58] where from Edmonds [59] we have so that in this case we recover the law of cosines. Ley-Koo and Bunge [58] note that the only contribution comes from the term with M = 0, "because the point of interest is on the polar axis" so that only the cosine-product term remains, which we confirmed by calculating both terms.
Introducing an addition theorem for the plane wave ( [60], p. 671, Equation (B.44)) helps us to do the first angular integral in each of the terms in Equation (17), All expressions that follow should in principle be multiplied by this factor of i, but, since we take the absolute square of sums of these transition amplitudes to get the cross section, we ignore this factor. We follow Ley-Koo and Bunge [58] in replacing dΩ 2 -the differential solid angle aroundr 2 in a frame of reference in which r 1 is taken as the polar axis-by dΩ 12 = sin θ 12 dθ 12 dφ 12 . One can immediately integrate over dφ 12 . They change variables to cos θ 12 = r 2 1 + r 2 2 − r 2 12 / (2r 1 r 2 ) giving sin θ 12 dθ 12 = (−2r 12 ) dr 12 / (2r 1 r 2 ), but one may also change variables to cos θ 12 = u 12 giving sin θ 12 dθ 12 = du 12 and simply do that integral using integrals we have not found in the literature: 1 The last integrals to do for the cross section are These are easily done and the derivatives in Equation (16) taken, providing the core of the numerical Maxwell-Boltzmann integral in Equation (13) for the rate coefficient α RA for producing the negative hydrogen ion H − or its antimatter equivalent, the positive antihydrogen ion H + [16].

Doing the r 1 Integrals Last
The conventional path to crafting an analytical rate coefficient α RA for producing the positive antihydrogen ion H + , or its matter equivalent, would be to integrate pair-products of terms in the analytical results for the cross section found from Equation (23) above, or as given in Keating's master's thesis [21], reproduced in Equation (2). This latter form organizes the sums of terms in the cross section most compactly, but integrating every pair-product of every term in even this set (Equation (2)) would require some two-dozen analytical integrals such as few of which are easily done. Instead of this obvious approach, we take the road less traveled and take these integrals in reverse order because of the uniformity of the integrands in Equation (23). The downside of this novel approach is that we must form the product of distinct radial integrals rather than squaring the analytical result of the result of a single integration, and there are many dead ends on a path to integrating over both of these radial variables after integrating over the equivalent of x in Equation (13). We did, however, find a means to do so, as follows.
We first take the derivatives in Equation (16) of Equation (23) to obtain the requisite terms of Equation (11), after substituting the conventional Bessel function for the spherical Bessel function , p. 673, Equation (C.2)): and The first term of I 11 may be combined with third-to-last term of I 21 , as may the third term of I 11 , with the last term of I 21 . The second term of I 11 and the second-to-last of I 21 cancel, thus one is left with seven terms comprising three powers of r 1 with two kinds of exponentials, a considerably uniform set of analytic functions to be integrated. If one wishes to do the Maxwell-Boltzmann integral of products of such functions first, the products must be written as unique integrals. That is, the rate coefficient will be sums of terms B (T, a, σ, h, τ, j, β, γ, λ, α, γ where we replace β in the denominators with d and f to allow for multiplication of sums of terms whose positrons are exchanged by theP 12 permutation operator in Equation (9) for the two identical positrons α ↔ β, as well as different values of α and β arising from the various terms in the adjoining sum over terms in the wave function. The latter also explains the need to distinguish γ b from γ in the second denominator. It turned out that the a in the exponential was unneeded for the present problem, but we have left it in for those who might need this sort of integral for another problem and later set it to equal one. One may perform the x and then the r 1 integrals as they stand, but the resulting expression did not allow for integration over R 1 . This is, of course, why one normally would never do the integrals in this order if it could be avoided. However, the hope of an overall simpler summing of products of terms if this unusual and difficult integration order is successful eventually found fruit in a series approach. We first express the Bessel functions in terms of the hypergeometric function [61-63] ; − 1 4 (kr 1 ) 2 (29) and combine their product as [64,65] 1 After expanding The r 1 integral, including coefficients from the second line of Equation (28), for each term in the m-sum follows as [67] We ignored the restriction against R 1 being an element of the reals under the assumption that this result could likely be considered as a distribution whose integral would smooth out any singularities arising from this restriction.
whose ratios have the following values for m an integer: Integrating the one remaining term (under the restriction (j) a, σ, h, τ, j, β, γ, λ, α, γ b , ps, pt) = 6.811556×10 −20 1.10114 × 10 6 √ T √ π (2 2 c 2 ) The rate coefficient is then a quadruple sum. First, we have a sum over the seven terms in I 11 + I 21 of Equations (26) and (27) plus their seven permutations α ↔ β, all taking on appropriate values of σ, h, and ps as used in Equation (28), and being multiplied by the appropriate numerators in Equations (26) and (27) such as 48γ(β + λ). Second, we have a sum over the corresponding 14 terms of the second factor that takes on the values of τ, j, and pt, multiplied by the corresponding numerators. Third, we have a sum over the n terms in the wave function having differing values for the parameters α, β, and γ, and that term's coefficient c whose value is found by diagonalizing the Hamiltonian matrix in order to minimize the ground state energy and then normalized. We display the first and last elements of the fifty-term version we used for both analytical and numerical calculations in case readers wish to check it against their own derivations: α(1) = 0.6878357596671101, · · · , α(50) = 0.37080904876118337 β(1) = 1.2354854281591452, · · · , β(50) = 1.1073078257847453 γ(1) = 0.012984468708341133, · · · , γ(50) = 0.28320160279252 c(1) = −2.534888772248287, · · · , c(50) = 0.02873436866901307 .
The final sum is again over this same wave-function set, but this time associated with the second factor τ. Because the off-diagonal terms of these four sums are symmetrical, it is more efficient computationally to account for that. Table 1 shows the comparison of the present analytical rate coefficient α RA for attaching a positron to antihydrogen to form H + to our prior rate coefficient found via numerical integration [16], at temperatures ranging from 1 to 400 K. There we found that, for a fully two-positron 200-term wave function [18] composed of explicitly correlated exponentials and for T

Comparison with Numerical Integration
The r 1 integral, including coefficients from the second line of (28), for each term in the m-sum follows as [63] We ignored the restriction against R 1 being an element of the reals under the assumption that this integral could likely be considered as a distribution whose integral would smooth out any singularities arising from this restriction. That assumption paid off. Let us consider the gamma functions of negative integer arguments that append each of the 2 F 3 functions [64], whose ratios have the following values for m an integer: Integrating the one remaining term (under the restriction (j) gives us our final result, The rate coefficient is then a quadruple sum. First we have a sum over the seven terms in I 11 + I 21 of (26) and (27) plus their seven permutations α ↔ β, all taking on appropriate values of σ, h ps as used in eq. (28), and being multiplied by the appropriate numerators in (26) and (27) such as 48γ(β + λ). Second we have a sum over the corresponding 14 terms of the second factor that takes on the values of τ , j, pt, multiplied by the corresponding numerators. Third we have a sum over the n terms in the wave function having differing values for the parameters α, β, and γ, and that term's coefficient c whose value is found by diagonalizing the Hamiltonian matrix in order to minimize the ground state energy and then normalized. We display the first and last elements of the fifty-term version we used for both analytical and numerical calculations in case readers wish to check it against their own derivations: α(1) = 0.6878357596671101, · · · , α(50) = 0.37080904876118337 β(1) = 1.2354854281591452, · · · , β(50) = 1.1073078257847453 γ(1) = 0.012984468708341133, · · · , γ(50) = 0.28320160279252 c(1) = −2.534888772248287, · · · , c(50) = 0.02873436866901307 .
The final sum is again over this same wave-function set, but this time associated with the second factor τ . Because the off-diagonal terms of these four sums are symmetrical, it is more efficient computationally to account for that. Table 1 shows the comparison of the present analytical rate coefficient α RA for attaching a positron to antihydrogen to form H + to our prior rate coefficient found via numerical integration [16], at temperatures ranging from 1 to 400 K. There we found, for a fully two-positron 200term wave function [18] composed of explicitly correlated exponentials, that for T 6K the rate coefficient is essentially linear and may be fit by α RA = 0.001050 × 10 −15 cm 3 s −1 T K −1 . One sees the potential for this linear behavior in the analytic result (33) in the factor outside of the sum, but only if the m = 0 term is the dominant contributor to the sum at this temperature. For a fifty-term wave function, the first 6 K, the rate coefficient is essentially linear and may be fit by α RA = 0.001050 × 10 −15 cm 3 s −1 T K −1 . One sees the potential for this linear behavior in the analytic result in Equation (34) in the factor outside of the sum, but only if the m = 0 term is the dominant contributor to the sum at this temperature. For a fifty-term wave function, the first three terms in the analytical sum are 0.00106039 − 1.2144278 × 10 −6 + 1.203002 × 10 −9 × 10 −15 cm 3 s −1 T K −1 , which sum to a value 0.00105917 × 10 −15 cm 3 s −1 T K −1 and the numerical integral for this same fifty-term wave function gives 0.00102625 × 10 −15 cm 3 s −1 T K −1 , a 3% difference.

Comparison with numerical integration
As the temperature increases, this series result that involves powers of the temperature will eventually fail. For a Z = 1 ion, we find that there is no convergence of the m series at T=700 K. For higher Z, higher temperatures are likely possible.

Discussion and Concluding Remarks
We found the analytic form for the rate coefficient α RA for attaching a positron to antihydrogen to form H + as a series convergent for temperatures of 400 K and below, which may be used to estimate formation rates. As our result is for attachment from the 1s antihydrogen state, it depends on the trapped antihydrogen having been held for long enough to ensure that it will have decayed to the ground state from the likely excited state in which it is formed [7,8].
If the positron plasma is held in a Penning trap of the type used to form antihydrogen, such as those reviewed in [69] at a density of n e = 10 16 m −3 in a magnetic field of 1 T at a sub-mm plasma radius, then the positron speeds due to rotation of the plasma can be neglected. The temperature range currently used to form and trap antihydrogen is in the range of 10s of K or lower. At 100 K, the rate coefficient is α RA = 10 × 10 −17 cm 3 s −1 , and if there is unit overlap between the positron plasma and the antihydrogen, then the reaction rate, the product n e α RA , at this temperature would be 10 × 10 −7 s −1 per antihydrogen atom. As the temperature falls to 10 K, the reaction rate falls in an essentially linear manner to be 1 × 10 −7 s −1 per antihydrogen atom. This might just be observable, given the long antihydrogen storage times achieved by ALPHA [8], if all ALPHA's antiprotons could be converted into trapped Hs, while still allowing the anti-atoms to interact with warm positron clouds. Cold p numbers, and hopefully those of trapped H as well, will likely be enhanced by around a factor of 10 2 within the next three years as CERN's AD facility is enhanced by the addition of a further storage ring, ELENA (see, e.g., [70]) that will deliver antiprotons to experiments at an energy nearer 100 keV than the current 5 MeV.