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

#### **Jack C. Straton**

Department of Physics, Portland State University, Portland, OR 97207-0751, USA; straton@pdx.edu

Received: 16 March 2020; Accepted: 10 April 2020; Published: 20 April 2020

**Abstract:** To overcome the numerical difficulties inherent in the Maxwell–Boltzmann integral of the velocity-weighted cross section that gives the radiative attachment rate coefficient *αRA* for producing the negative hydrogen ion H− or its antimatter equivalent, the positive antihydrogen ion <sup>H</sup>+, we found the analytic form for this integral. This procedure is useful for temperatures below 700 K, the region for which the production of <sup>H</sup><sup>+</sup> 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.

**Keywords:** antihydrogen; radiative attachment; photodetachment; antihydrogen ion; analytical; hydrogen ion

#### **1. Introduction**

The Antiproton Decelerator (AD) facility at CERN [1] has provided the foundation for a variety of experiments (e.g., [2–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–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–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 <sup>H</sup><sup>+</sup> 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–14]. The authors of [15,16] calculated the cross section and rate coefficient for the radiative attachment of a second positron to create the <sup>H</sup><sup>+</sup> ion,

$$
\blacksquare \blacksquare (1s) + \epsilon^{+} \to \blacksquare \blacksquare^{+} \left(1s^{2} \: 1^{\text{c}}\right) + \hbar \omega \; . \tag{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<sup>−</sup> or its antimatter equivalent, the positive antihydrogen ion <sup>H</sup>+, 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 *L*th derivatives of the Laplace transform of the spherical Bessel function *j*1(*kr*),

$$\mathcal{L}\_L(p;k) = \begin{cases} \frac{1}{k^2} \left[ k - p \tan^{-1} \left( \frac{k}{p} \right) \right] & L = 0 \\\\ \frac{1}{k^2} \left[ -\frac{kp}{p^2 + k^2} + \tan^{-1} \left( \frac{k}{p} \right) \right] & L = 1 \\\\ \frac{i(L-2)!}{2k^2} \left[ \frac{p + ikL}{(p + ik)^L} - \frac{p - ikL}{(p - ik)^L} \right] & \\\\ \frac{-(L-2)!}{k \left( k^2 + p^2 \right)^L} \sum\_{j=0}^L (2) \binom{j}{j} \frac{j(L+1)}{j+1} \left( -1 \right)^{j/2} k^j p^{L-j} & L = 2, 3, 4, \dots \end{cases} \tag{2}$$

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<sup>+</sup>*he*−*<sup>σ</sup>rj*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<sup>+</sup>*he*−*<sup>σ</sup>rj*1(*kr*) *R*3/2+*<sup>j</sup> e*−*<sup>τ</sup>rj*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.

#### **2. 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–42] and more recently [43–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

$$
\hbar\omega + \mathcal{H}^- \left( \mathbf{1} s^{2} \, ^1 \mathcal{S}\_o^{\varepsilon} \right) \to \mathcal{H}^- \left( \mathbf{1} s \, \mathbf{k} \mathbf{p} \, ^1 P\_o^{\rho} \right). \tag{3}
$$

Following Ghoshal and and Ho [50], we put the (length gauge1) cross section for photodetachment (or photoionization), *σPI*, laid out by Chandrasekhar [24,25] in atomic units, and obtain

$$\begin{array}{rcl}\sigma\_{PI} &=& \frac{2p\omega a a\_0^2}{3} \left| \left< \psi\_f \left| \not\mathbf{k} \cdot (\mathbf{r}\_1 + \mathbf{r}\_2) \right| \psi\_i \right> \right|^2 \\ &=& 6.81156 \times 10^{-20} \text{ cm}^2 k \left( k^2 + 2I \right) \left| \left< \psi\_f \left| z\_1 + z\_2 \right| \psi\_i \right> \right|^2 \end{array} \tag{4}$$

where *α* is the fine structure constant, *a*<sup>0</sup> 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

<sup>1</sup> 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.

$$2\omega = 2\left(2\pi\nu\right) = 2\frac{\hbar c}{\lambda} = 2\left(\frac{\hbar^2 k^2}{2} + I\right) \equiv \left(p^2 + \gamma\_0^2\right). \tag{5}$$

In the matrix element for photodetachment

$$
\mu\_{Pl} = \int \psi\_f^\*(z\_1 + z\_2) \psi\_l d\tau \tag{6}
$$

in Equation (4), *ψ<sup>f</sup>* 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,

$$\psi\_f = \frac{1}{\sqrt{2\pi}} (e^{ikz\_1 - r\_2} + e^{ikz\_2 - r\_1}) \tag{7}$$

Accurate two-positron initial-state wave functions may take the form

$$
\psi\_H(r\_1, r\_2, r\_{12}) = \frac{1}{\sqrt{2}} \left( 1 - \mathcal{P}\_{12} \right) e^{-ar\_1 - \beta r\_2 - \gamma r\_{12}} \sum\_{l, m, n} c\_{l \text{min}} s^l t^{2m} u^n \tag{8}
$$

where *P*ˆ <sup>12</sup> is the permutation operator for the two identical positrons *α* ↔ *β*, with Hylleraas coordinates [53] given by *s* = *r*<sup>1</sup> + *r*2, *t* = *r*<sup>1</sup> − *r*2, and *u* = *r*<sup>12</sup> ≡ |**r**<sup>1</sup> − **r**2|. One may also express this as sums of powers of *r*<sup>1</sup> and *r*<sup>2</sup> 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,

$$
\psi\_{TS}(r\_1, r\_2, r\_{12}) = \frac{1}{\sqrt{2}} \left(1 - \mathcal{P}\_{12}\right) \sum\_{k} c\_k e^{-a\_k r\_1 - \beta\_k r\_2 - \gamma\_k r\_{12}} \tag{9}
$$

where the parameters in the exponentials are generated in a quasi-random fashion,

$$\begin{array}{rcl}\mathfrak{a}\_{k} &=& \eta \left( (A\_{2} - A\_{1}) \frac{1}{2} \langle k(k+1) \rangle \sqrt{2} + A\_{1} \right) \\\mathfrak{b}\_{k} &=& \eta \left( (B\_{2} - B\_{1}) \frac{1}{2} \langle k(k+1) \rangle \sqrt{3} + B\_{1} \right) \\\mathfrak{m}\_{k} &=& \eta \left( (G\_{2} - G\_{1}) \frac{1}{2} \langle k(k+1) \rangle \sqrt{5} + G\_{1} \right) \end{array} \tag{10}$$

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*<sup>1</sup> = 0.2380, *A*<sup>2</sup> = 1.3240, *B*<sup>1</sup> = 0.9800, *B*<sup>2</sup> = 1.3290, *G*<sup>1</sup> = −0.0720, *G*<sup>2</sup> = 0.288, and *<sup>η</sup>* = <sup>1</sup> − 2.458 × <sup>10</sup>−7. The quasi-random assignment of the 50 values for each of *<sup>α</sup>k*, *<sup>β</sup>k*, and *<sup>γ</sup><sup>k</sup>* 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 *ck* 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 *γ* = *γ<sup>k</sup>* in the sum in Equation (9), the matrix element for ionizing either positron under the influence of the length dipole operator (*z*<sup>1</sup> + *z*2) is the sum of four terms:

$$\mu\_{PI} = 2\left(I\_{11} + I\_{21}\right)\frac{\lambda^{3/2}}{\sqrt{2\pi}} \quad = \int\int d^3 \mathbf{x}\_1 d^3 \mathbf{x}\_2 e^{-\alpha r\_1 - \beta r\_2 - \gamma r\_{12}} \left(z\_1 + z\_2\right) \frac{\lambda^{3/2}}{\sqrt{\pi}} \frac{1}{\sqrt{2}} \quad , \tag{11}$$

$$\quad \times \quad \left[e^{-\lambda r\_2} e^{\bar{\beta}\bar{\pi}\_1} + e^{-\lambda r\_1} e^{\bar{\beta}\bar{\pi}\_2}\right]$$

where the first factor of two comes from *I*<sup>11</sup> = *I*<sup>22</sup> and so on, whose subscripts *j* refer to *zj* in the dipole operator and in the plane wave, respectively, and we have factored out the coefficient *<sup>λ</sup>*3/2 <sup>√</sup>2*<sup>π</sup>* 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 *rj*.

The cross section for radiatively attaching a second positron to H (1*s*) to create the (1*s*2 1*S<sup>e</sup>* ) state of the <sup>H</sup><sup>+</sup> 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 (2*p*2 3*P<sup>e</sup>* ) metastable H<sup>−</sup> state from H (2*s*, 2*p*). For the (1*s*2 1*S<sup>e</sup>* ) case, we have [16],

$$
\sigma\_{RA}(k) = \frac{\mathcal{g}\_1 p\_\omega^2}{\mathcal{g}\_2 p\_\varepsilon^2} \sigma\_{PI} = \frac{6a^2 \left(k^2 + \gamma\_0^2\right)^2}{12 \cdot 2^2 k^2} \sigma\_{PI} \tag{12}
$$

where *g*1/*g*<sup>2</sup> = 6/12 is the statistical weight ratio. Here, the photon momentum relative to the ion is given by *p<sup>ω</sup>* = *h*¯ *ω*/*c* = (*k*<sup>2</sup> + *γ*<sup>2</sup> <sup>0</sup>)/2*c*, and *pe* is the positron momentum *k*. Note that *c* in atomic units is the inverse of the fine structure constant *α*.

To estimate formation rates of <sup>H</sup>+, it is helpful to calculate the positron attachment to <sup>H</sup> 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,

$$\begin{split} a\_{RA}(T) = \langle v\sigma\_{RA} \rangle &= 4\pi \int\_0^\infty dv \, v\sigma\_{RA} \left( k\left( v \right) \right) \left( \frac{m}{2\pi k\_B T} \right)^{3/2} v^2 \exp\left[ -mv^2 / \left( 2k\_B T \right) \right] \\ &= \sqrt{\frac{8k\_B T}{m\_\varepsilon}} \frac{1}{\sqrt{\pi}} \int\_0^\infty dx \, x \frac{g\_{Pl}}{g\_{RA}} \frac{P\_\omega^2}{p\_\varepsilon^2} \sigma\_{PI} \left( \sqrt{2k\_B T x} \right) \exp\left[ -x \right] \ . \end{split} \tag{13}$$

whose overall coefficient is

$$
\sqrt{\frac{8k\_B T}{m\_c}} = \sqrt{\frac{8 \times 8.61733262 \times 10^{-5} \text{ eV } K^{-1} T}{0.51099895000 \times 10^6 \text{ eV}}} 299792458 \times 10^2 \text{ cm/s} \tag{14}
$$

$$
= 1.10113894 \times 10^6 \text{ cm/s} \sqrt{K^{-1} T}
$$

where the temperature *T* is given in Kelvins *K*.

#### **3. Evaluating Integrals**

The algebra in each case is greatly reduced by making the replacement

$$e^{-\gamma r\_{12}}e^{-r\_2(\notin+\lambda)} = \left(-\frac{\partial}{\partial \gamma}\right)\left(-\frac{\partial}{\partial \lambda}\right)\frac{e^{-r\_2(\notin+\lambda)}}{r\_2}\frac{e^{-\gamma r\_{12}}}{r\_{12}}\tag{15}$$

in each term, with the notation

$$I\_{j1} = \left(\frac{\partial}{\partial \gamma}\right) \left(\frac{\partial}{\partial \lambda}\right) \mathcal{R}\_{0j1} \tag{16}$$

where

$$\begin{array}{rcl} R\_{011} &=& \int\int d^3x\_1 d^3x\_2 e^{-ar\_1} \frac{e^{-r\_2(\beta+\lambda)}}{r\_2} \frac{e^{-\gamma r\_{12}}}{r\_{12}} \left(z\_1 e^{ikz\_1}\right) \\\ R\_{021} &=& \int\int d^3x\_1 d^3x\_2 e^{-ar\_1} \frac{e^{-r\_2(\beta+\lambda)}}{r\_2} \frac{e^{-\gamma r\_{12}}}{r\_{12}} \left(z\_2 e^{ikz\_1}\right) \end{array} \tag{17}$$

*Atoms* **2020**, *8*, 13

For the cross-terms *z*2*eikz*<sup>1</sup> = *r*<sup>2</sup> cos *θ*2*eikr*<sup>1</sup> cos *<sup>θ</sup>*<sup>1</sup> , we need the conversion given in Ley-Koo and Bunge [58]

$$\mathbf{Y}\_{k,\mathfrak{m}}\left(\theta\_2,\phi\_2\right) = \sum\_{M=-k}^{k} \left[\mathfrak{D}\_{\mathfrak{m},M}^{(k)}\left(0,\theta\_1,\phi\_1\right)\right]^\* \mathbf{Y}\_{k,M}\left(\theta\_{12},\phi\_{12}\right) \tag{18}$$

where from Edmonds [59] we have

$$\begin{aligned} \left(\mathbf{Y}\_{1,0}\left(\theta\_2,\phi\_2\right)\right)^{\*} &=& \sum\_{M=-1}^{1} \left[\mathfrak{D}\_{0,M}^{(1)}\left(0,\theta\_1,\phi\_1\right)\right]^{\*} \mathbf{Y}\_{1,M}\left(\theta\_{12},\phi\_1\mathbf{1}\_2\right) \\ &=& \sum\_{M=-1}^{1} \left[\sqrt{\frac{4\pi}{3}}\mathbf{Y}\_{1,M}\left(\theta\_{1\prime},\phi\_1\right)\right]^{\*} \mathbf{Y}\_{1,M}\left(\theta\_{12\prime},\phi\_{12}\right) \end{aligned} \tag{19}$$

so that in this case

$$\begin{aligned} \cos\theta\_2 &= 2\sqrt{\frac{\pi}{3}}\mathcal{Y}\_{1,0}\left(\theta\_2,\phi\_2\right) &= \begin{array}{c} \frac{4\pi}{3}\Sigma^1\_{M=-1}\mathcal{Y}^\*\_{1,M}\left(\theta\_1,\phi\_1\right)\mathcal{Y}\_{1,M}\left(\theta\_{12},\phi\_{12}\right) \\ &= \mathcal{P}\_1\left(\cos\theta\_2\right) \\ \end{array} \tag{20} &= \begin{array}{c} \cos\theta\_1\cos\theta\_{12} + \sin\theta\_1\sin\theta\_{12}\cos\left(\phi\_1+\phi\_{12}\right) \end{array} \tag{21} \end{aligned} \tag{20} $$

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),

$$\int d\Omega\_1 \left( r\_1 P\_1 \left( \cos \theta\_1 \right) \right) \left( \sum\_{l=0}^{\infty} (2l+1) i^l j\_l(kr\_1) P\_l \left( \cos \theta\_1 \right) \right) = \left( r\_1 \frac{2}{(2+1)} 2\pi \right) \left( (2+1) i^1 j\_1(kr\_1) \right) \quad. \tag{21}$$

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 around *r*ˆ2 in a frame of reference in which **r**<sup>1</sup> is taken as the polar axis—by *d*Ω<sup>12</sup> = sin *θ*12*dθ*12*dφ*12. One can immediately integrate over *dφ*12. They change variables to cos *θ*<sup>12</sup> = *r*2 <sup>1</sup> + *<sup>r</sup>*<sup>2</sup> <sup>2</sup> − *<sup>r</sup>*<sup>2</sup> 12 / (2*r*1*r*2) giving sin *θ*12*dθ*<sup>12</sup> = (−2*r*12) *dr*12/ (2*r*1*r*2), but one may also change variables to cos *θ*<sup>12</sup> = *u*<sup>12</sup> giving sin *θ*12*dθ*<sup>12</sup> = *du*<sup>12</sup> and simply do that integral using integrals we have not found in the literature:

 <sup>1</sup> −1 *du*<sup>12</sup> *e* −*γ r*2 1−2*r*1*r*2*u*12+*r*<sup>2</sup> 2 *r*2 <sup>1</sup> − <sup>2</sup>*r*1*r*2*u*<sup>12</sup> + *<sup>r</sup>*<sup>2</sup> 2 <sup>=</sup> *<sup>e</sup>*−*γ*|*r*1−*r*2<sup>|</sup> <sup>−</sup> *<sup>e</sup>*−*γ*|*r*1+*r*2<sup>|</sup> *γr*1*r*<sup>2</sup> <sup>1</sup> −1 *du*<sup>12</sup> *e* −*γ r*2 1−2*r*1*r*2*u*12+*r*<sup>2</sup> 2 *r*2 <sup>1</sup> − <sup>2</sup>*r*1*r*2*u*<sup>12</sup> + *<sup>r</sup>*<sup>2</sup> 2 *<sup>u</sup>*<sup>12</sup> <sup>=</sup> <sup>1</sup> *γ*3*r*<sup>2</sup> 1*r*2 2 *e* <sup>−</sup>*γ*|*r*1−*r*2<sup>|</sup> + *e* −*γ*(*r*1+*r*2) *<sup>r</sup>*1*r*2*γ*<sup>2</sup> <sup>−</sup> *<sup>e</sup>* <sup>−</sup>*γ*|*r*1−*r*2<sup>|</sup> + *e* −*γr*1−*γr*<sup>2</sup> + *γ e* <sup>−</sup>*γ*(*r*1+*r*2) (*r*<sup>1</sup> <sup>+</sup> *<sup>r</sup>*2) <sup>−</sup> *<sup>e</sup>* <sup>−</sup>*γ*|*r*1−*r*2<sup>|</sup> <sup>|</sup>*r*<sup>1</sup> <sup>−</sup> *<sup>r</sup>*2<sup>|</sup> [*r*<sup>1</sup> <sup>&</sup>gt; 0, *<sup>r</sup>*<sup>2</sup> <sup>&</sup>gt; <sup>0</sup>] , (22) 1 *γr*1*r*<sup>2</sup> *∂ ∂a* <sup>1</sup> −1 *du*<sup>12</sup> *e* −*γ r*2 1−*a*2*r*1*r*2*u*12+*r*<sup>2</sup> 2 *a*=1 *<sup>u</sup>*<sup>12</sup> <sup>=</sup> <sup>1</sup> *γr*1*r*<sup>2</sup> *∂ ∂a* 1 *a*2*γ*4*r*<sup>2</sup> 1*r*2 2 *e* −*γ r*2 1−2*ar*2*r*1+*r*<sup>2</sup> 2 − *r* 2 <sup>1</sup> − 3*ar*2*r*<sup>1</sup> + *r* 2 2 *γ*2 + *<sup>a</sup>γ*2*r*1*r*<sup>2</sup> <sup>−</sup> <sup>3</sup> *r*2 <sup>1</sup> − <sup>2</sup>*ar*2*r*<sup>1</sup> + *<sup>r</sup>*<sup>2</sup> <sup>2</sup>*γ* − 3 + *e* −*γ r*2 1+2*ar*2*r*1+*r*<sup>2</sup> 2 × *r* 2 <sup>1</sup> + 3*ar*2*r*<sup>1</sup> + *r* 2 2 *γ*<sup>2</sup> + *ar*1*r*2*γ*<sup>2</sup> + 3 *r*2 <sup>1</sup> + <sup>2</sup>*ar*2*r*<sup>1</sup> + *<sup>r</sup>*<sup>2</sup> <sup>2</sup>*γ* + 3 *a*=1

The last integrals to do for the cross section are

$$\begin{array}{rcl} R\_{011} &=& 8\pi^2 \frac{2}{(\beta+\lambda)^2 - \gamma^2} \int\_0^\infty dr\_1 r\_1^2 j\_1 \left(kr\_1\right) \left(e^{-(a+\beta+\lambda)r\_1} - e^{-(a+\gamma)r\_1}\right) \\ R\_{021} &=& 8\pi^2 \frac{4}{((\beta+\lambda)^2 - \gamma^2)^2} \int\_0^\infty dr\_1 \left(j\_1 \left(kr\_1\right) \left(e^{-(a+\beta+\lambda)r\_1} - e^{-(a+\gamma)r\_1}\right)\right) \\ &+& 8\pi^2 \frac{4}{((\beta+\lambda)^2 - \gamma^2)^2} \int\_0^\infty dr\_1 \left(r\_1 j\_1 \left(kr\_1\right) \left((\beta+\lambda)e^{-(a+\beta+\lambda)r\_1} - \gamma e^{-(a+\gamma)r\_1}\right)\right) \\ &+& 8\pi^2 \frac{2}{(\beta+\lambda)^2 - \gamma^2} \int\_0^\infty dr\_1 \left(r\_1^2 j\_1 \left(kr\_1\right) e^{-(a+\beta+\lambda)r\_1}\right) \end{array} \right)$$

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<sup>−</sup> or its antimatter equivalent, the positive antihydrogen ion <sup>H</sup><sup>+</sup> [16].

#### **4. Doing the** *r*<sup>1</sup> **Integrals Last**

The conventional path to crafting an analytical rate coefficient *αRA* for producing the positive antihydrogen ion <sup>H</sup>+, 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

$$\int\_0^\infty \frac{e^{-\mathbf{x}} \left(\gamma p^2 + 2k\_B T \mathbf{x}\right)^3}{4\sqrt{2}(k\_B T)^{5/2} \mathbf{x}^{3/2}} \tan^{-1}\left(\frac{\sqrt{2k\_B T \mathbf{x}}}{p\_f}\right) \tan^{-1}\left(\frac{\sqrt{2k\_B T \mathbf{x}}}{p\_i}\right) d\mathbf{x} \tag{24}$$

and

$$\int\_0^\infty \frac{e^{-x} \left(\gamma\_0^2 + 2k\_B T x\right)^3 \tan^{-1}\left(\frac{\sqrt{2k\_B T x}}{p\_i}\right)}{4k\_B^2 T^2 x \left(p\_f^2 + 2k\_B T x\right)} dx \quad , \tag{25}$$

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 √1 *kr*<sup>1</sup> *π* <sup>2</sup> *J* <sup>3</sup> 2 (*kr*1) = *j*<sup>1</sup> (*kr*1) ([60], p. 673, Equation (C.2)):

$$\begin{array}{rcl} I\_{11} &=& 4\pi f\_{0}^{\prime\prime} \, dr\_{1} 4\pi \left(\frac{1}{\sqrt{k}} \sqrt{\frac{\pi}{2}} f\_{\frac{1}{2}} \left(kr\_{1} \right) r\_{1}^{3/2} \right) \\ & \times \quad \left(\frac{8\gamma \left(\oint + \lambda\right) \varepsilon^{-\left(a+\beta+\lambda\right)r\_{1}}}{\left(\left(\oint + \lambda\right)^{2} - \gamma^{2}\right)^{3}} - \frac{8\gamma \left(\oint + \lambda\right) \varepsilon^{-\left(a+\gamma\right)r\_{1}}}{\left(\left(\oint + \lambda\right)^{2} - \gamma^{2}\right)^{3}} + \frac{2r\_{1}\gamma \varepsilon^{-\left(a+\beta+\lambda\right)r\_{1}}}{\left(\left(\oint + \lambda\right)^{2} - \gamma^{2}\right)^{2}} + \frac{2r\_{1}\left(\left(\oint + \lambda\right) \varepsilon^{-\left(a+\gamma\right)r\_{1}}\right)}{\left(\left(\oint + \lambda\right)^{2} - \gamma^{2}\right)^{2}}\right) \end{array} \tag{26}$$

and

*I*<sup>21</sup> = 4*π* <sup>∞</sup> 0 *dr*14*π* √1 *k <sup>π</sup>* <sup>2</sup> *J* <sup>3</sup> 2 (*kr*1)*r* 3/2 1 48*γ* (*β* + *λ*)*e*−(*α*+*β*+*λ*)*r*<sup>1</sup> *r*2 <sup>1</sup> ((*β* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) <sup>4</sup> <sup>−</sup> <sup>48</sup>*γe*−(*α*+*γ*)*r*<sup>1</sup> (*<sup>β</sup>* <sup>+</sup> *<sup>λ</sup>*) *r*2 <sup>1</sup> ((*β* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) <sup>4</sup> <sup>+</sup> <sup>48</sup>*<sup>γ</sup>* (*<sup>β</sup>* <sup>+</sup> *<sup>λ</sup>*) <sup>2</sup> *e*−(*α*+*β*+*λ*)*r*<sup>1</sup> *r*<sup>1</sup> ((*β* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) 4 <sup>−</sup> <sup>48</sup>*γ*<sup>2</sup> (*<sup>β</sup>* <sup>+</sup> *<sup>λ</sup>*)*e*−(*α*+*γ*)*r*<sup>1</sup> *r*<sup>1</sup> ((*β* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) <sup>4</sup> <sup>+</sup> <sup>16</sup>*<sup>γ</sup>* (*<sup>β</sup>* <sup>+</sup> *<sup>λ</sup>*)*e*−(*α*+*β*+*λ*)*r*<sup>1</sup> ((*β* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) <sup>3</sup> <sup>+</sup> <sup>8</sup>*<sup>γ</sup>* (*<sup>β</sup>* <sup>+</sup> *<sup>λ</sup>*)*e*−(*α*+*γ*)*r*<sup>1</sup> ((*β* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) <sup>3</sup> <sup>+</sup> <sup>2</sup>*r*1*γe*−(*α*+*β*+*λ*)*r*<sup>1</sup> ((*β* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) 2 . (27)

The first term of *I*<sup>11</sup> 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*<sup>11</sup> and the second-to-last of *I*<sup>21</sup> cancel, thus one is left with seven terms comprising three powers of *r*<sup>1</sup> 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

$$\begin{split} B(\mathbf{T},\mathbf{a},\sigma,\mathbf{b},\tau\_{i}\boldsymbol{\beta},\gamma,\lambda,\mathbf{a},\tau\_{0},\tau\_{1},\tau\_{0},\mathbf{p},\mathbf{p},\mathbf{t}) &= \frac{6.81159\varepsilon \times 10^{-26} \left(1.1014 \times 10^{6} \sqrt{T}\right)\_{\frac{\mathbf{r} \times \mathbf{r}}{\mathbf{r} \times \mathbf{r}}}}{\sqrt{\pi} \left(2^{2} \boldsymbol{\varepsilon}^{2}\right)\_{\frac{\mathbf{r}}{\mathbf{r} \times \mathbf{r}}}} \frac{\mathbf{r}\prime\prime}{g\varepsilon\Lambda} \\ &\times 128\pi^{4}\lambda^{3} \int\_{0}^{\infty} d\Omega\_{1} \int\_{0}^{\infty} d\nu\_{1} \int\_{0}^{\infty} d\tau\_{1} \int\_{0}^{\infty} d\tau\_{2} \frac{\left(2k\_{1}\mathrm{Tr}\,\tau\_{1} + \tau\_{0}\right)^{2} J\_{\frac{\mathbf{r}}{\mathbf{r}}} \left(\sqrt{2}\sqrt{4\pi\chi\prime\prime}x\_{\mathrm{Tr}}\right) J\_{\frac{\mathbf{r}}{\mathbf{r}}\prime} \left(\sqrt{2}\sqrt{4\mathrm{Tr}\,\tau\_{1}\mathrm{Tr}\,\tau\_{2}}\right) e^{-\sigma\varepsilon\_{1} - \tau\Omega\_{1} - \varepsilon\varpi}}{k\_{B}T \left((d+\lambda)^{2} - \gamma^{2}\right)^{\mathrm{Re}} \left((f+\lambda)^{2} - \gamma\_{0}^{2}\right)^{\mathrm{Re}}} \end{split} \tag{28}$$

where we replace *β* in the denominators with *d* and *f* to allow for multiplication of sums of terms whose positrons are exchanged by the *P*ˆ <sup>12</sup> 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 *γ<sup>b</sup>* 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*<sup>1</sup> 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]

$$J\_{\frac{3}{2}}\left(kr\_{1}\right) = \frac{1}{\Gamma\left(\frac{5}{2}\right)} \left(\frac{kr\_{1}}{2}\right)^{3/2} \,\_0F\_{1}\left(\cdot; \frac{5}{2}; -\frac{1}{4}\left(kr\_{1}\right)^{2}\right) \tag{29}$$

and combine their product as [64,65]

$$\begin{array}{c} \frac{1}{\Gamma\left(\frac{5}{2}\right)} \left(\frac{r\_{1}\sqrt{2k\_{B}Tx}}{2}\right)^{3/2} \,\_{0}F\_{1}\left(\cdot\frac{5}{2}; -\frac{1}{4}\,\left(r\_{1}\sqrt{2k\_{B}Tx}\right)^{2}\right) \frac{1}{\Gamma\left(\frac{5}{2}\right)} \left(\frac{R\_{1}\sqrt{2k\_{B}Tx}}{2}\right)^{3/2} \,\_{0}F\_{1}\left(\cdot\frac{5}{2}; -\frac{1}{4}\,\left(R\_{1}\sqrt{2k\_{B}Tx}\right)^{2}\right) \\\ =\quad \frac{4\sqrt{2}\left(r\_{1}\sqrt{k\_{B}Tx}\right)^{3/2}\left(R\_{1}\sqrt{k\_{B}Tx}\right)^{3/2} \sum\_{m=0}^{\infty} \frac{2^{-m}\left(-r\_{1}^{2}\right)^{m}(k\_{B}Tx)^{m}}{m!\left(\frac{5}{2}\right)\_{m}} \,\_{2}F\_{1}\left(-m-\frac{3}{2},-m; \frac{5}{2}; \frac{R\_{1}^{2}}{r\_{1}^{2}}\right) \end{array} \tag{30}$$

After expanding

$$\left(2k\_B T \mathbf{x} + \gamma\_0^2\right)^3 = 8k\_B^3 T^3 \mathbf{x}^3 + 12k\_B^2 T^2 \mathbf{x}^2 \gamma 0^2 + 6k\_B T \mathbf{x} \gamma 0^4 + \gamma 0^6 \tag{31}$$

the integral over powers *n* of *x* in this sum is [66] (p. 364 No. 3.381.4)

$$\int\_0^\infty d\mathbf{x} \epsilon^{-a\mathbf{x}} \mathbf{x}^{m+n+\nu} = a^{-m-n-\nu-1} \Gamma(m+n+\nu+1)$$

*Atoms* **2020**, *8*, 13

The *r*<sup>1</sup> integral, including coefficients from the second line of Equation (28), for each term in the *m*-sum follows as [67]

*B*<sup>2</sup> = *π*32 19 <sup>2</sup> <sup>−</sup>*mλ*<sup>3</sup> *<sup>a</sup>*−*m*−*n*<sup>−</sup> <sup>5</sup> <sup>2</sup> (*kBT*) *m*+ 3 2 Γ *m* + *n* + <sup>5</sup> 2 9kB*Tm*! 5 2 *<sup>m</sup>* ((*d* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) (*<sup>f</sup>* <sup>+</sup> *<sup>λ</sup>*)<sup>2</sup> <sup>−</sup> *<sup>γ</sup>*b2 *<sup>R</sup>j*+<sup>3</sup> <sup>1</sup> *e* −*τR*1 <sup>∞</sup> 0 *dr*1*rh*+<sup>3</sup> 1 −*r* 2 1 *<sup>m</sup> e* <sup>−</sup>*σr*<sup>1</sup> <sup>2</sup>*F*<sup>1</sup> <sup>−</sup>*<sup>m</sup>* <sup>−</sup> <sup>3</sup> 2 , −*m*; 5 2 ; *R*2 1 *r*2 1 <sup>=</sup> *<sup>π</sup>*7/2(−1)*<sup>m</sup>*2*m*<sup>+</sup> <sup>25</sup> <sup>2</sup> *λ*3*Rj*+<sup>3</sup> <sup>1</sup> *e* <sup>−</sup>*τR*<sup>1</sup> *<sup>a</sup>*−*m*−*n*<sup>−</sup> <sup>5</sup> <sup>2</sup> *σ*−*h*−2*m*−<sup>4</sup>(*kBT*) *m*+ 1 2 Γ *<sup>m</sup>* <sup>+</sup> *<sup>n</sup>* <sup>+</sup> <sup>5</sup> 2 <sup>×</sup> <sup>1</sup> 9*m*! 5 2 *<sup>m</sup>* Γ −*<sup>m</sup>* − <sup>3</sup> 2 ((*d* + *λ*)<sup>2</sup> − *γ*<sup>2</sup>) (*<sup>f</sup>* <sup>+</sup> *<sup>λ</sup>*)<sup>2</sup> <sup>−</sup> *<sup>γ</sup>*b2 (32) × −2(*m* + 1) Γ(−2*m* − 3) <sup>Γ</sup>(−*m*) <sup>Γ</sup>(*<sup>h</sup>* <sup>+</sup> <sup>2</sup>*<sup>m</sup>* <sup>+</sup> <sup>4</sup>) <sup>2</sup>*F*<sup>3</sup> <sup>−</sup>*<sup>m</sup>* <sup>−</sup> <sup>3</sup> 2 , −*m*; 5 2 , − *h* <sup>2</sup> <sup>−</sup> *<sup>m</sup>* <sup>−</sup> <sup>3</sup> 2 , − *h* <sup>2</sup> <sup>−</sup> *<sup>m</sup>* <sup>−</sup> 1; <sup>1</sup> 4 *σ*2*R*<sup>2</sup> 1 + 3*R*<sup>4</sup> 1 − 1 *R*2 1 − *h* <sup>2</sup> <sup>−</sup>*mσh*+2*m*+<sup>4</sup> <sup>−</sup> (*<sup>h</sup>* <sup>+</sup> <sup>2</sup>)Γ(*<sup>h</sup>* <sup>+</sup> <sup>1</sup>) cos *<sup>π</sup><sup>h</sup>* <sup>2</sup> + *πm h* + 2*m* + 7 Γ(−*h* − 2*m* − 5) <sup>Γ</sup>(−*m*) <sup>2</sup>*F*<sup>3</sup> *h* <sup>2</sup> <sup>+</sup> <sup>1</sup> 2 , *h* <sup>2</sup> <sup>+</sup> 2; <sup>1</sup> 2 , *h* <sup>2</sup> <sup>+</sup> *<sup>m</sup>* <sup>+</sup> 3, *<sup>h</sup>* <sup>2</sup> <sup>+</sup> *<sup>m</sup>* <sup>+</sup> <sup>9</sup> 2 ; 1 4 *σ*2*R*<sup>2</sup> 1 <sup>−</sup> (*<sup>h</sup>* <sup>+</sup> <sup>3</sup>)*σ*Γ(*<sup>h</sup>* <sup>+</sup> <sup>2</sup>) sin *<sup>π</sup><sup>h</sup>* <sup>2</sup> + *πm* − <sup>1</sup> *R*2 1 (*h* + 2*m* + 8) Γ(−*h* − 2*m* − 6) <sup>Γ</sup>(−*m*) <sup>2</sup>*F*<sup>3</sup> *h* <sup>2</sup> <sup>+</sup> 1, *<sup>h</sup>* <sup>2</sup> <sup>+</sup> <sup>5</sup> 2 ; 3 2 , *h* <sup>2</sup> <sup>+</sup> *<sup>m</sup>* <sup>+</sup> <sup>7</sup> 2 , *h* <sup>2</sup> <sup>+</sup> *<sup>m</sup>* <sup>+</sup> 5; <sup>1</sup> 4 *σ*2*R*<sup>2</sup> 1 ⎤ ⎥ ⎥ <sup>⎦</sup> . (*h*) + 1 > 0 ∧ (*h* + 2*m*) + 4 > 0 ∧ (*σ*) > 0 ∧ *R*2 <sup>1</sup> <sup>∈</sup>/ <sup>R</sup> ∨ *R*2 1 ≤ 0 <sup>∧</sup> (*R*<sup>1</sup> <sup>∈</sup>/ <sup>R</sup> <sup>∨</sup> ( (*R*1) <sup>=</sup> <sup>0</sup> <sup>∧</sup> *<sup>R</sup>*<sup>1</sup> = 0))

We ignored the restriction against *R*<sup>1</sup> 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. That assumption paid off. Let us consider the gamma functions of negative integer arguments that append each of the <sup>2</sup>*F*<sup>3</sup> functions ([66], p. 946 No. 8.334.3; [68]),

$$\begin{array}{rcl} \Gamma(-2m-3) &=& -\frac{\pi \csc(\pi(2m+3))}{\Gamma(2m+4)} \\\\ \Gamma(-m) &=& -\frac{\pi \csc(\pi m)}{\Gamma(m+1)} \\\\ \Gamma(-h-2m-5) &=& -\frac{\pi \csc(\pi(h+2m+5))}{\Gamma(h+2m+6)} \\\\ \Gamma(-h-2m-6) &=& -\frac{\pi \csc(\pi(h+2m+6))}{\Gamma(h+2m+7)} \end{array} \tag{33}$$

whose ratios have the following values for *m* an integer:

$$-\frac{\pi\csc(\pi(2m+3))}{\Gamma(2m+4)} \left(-\frac{\pi\csc(\pi m)}{\Gamma(m+1)}\right)^{-1} = \frac{(-1)^{m+1}\Gamma(m+1)}{2\Gamma(2m+4)}$$

$$-\frac{\pi\csc(\pi(h+2m+5))}{\Gamma(h+2m+6)} \left(-\frac{\pi\csc(\pi m)}{\Gamma(m+1)}\right)^{-1} = 0$$

$$-\frac{\pi\csc(\pi(h+2m+6))}{\Gamma(h+2m+7)} \left(-\frac{\pi\csc(\pi m)}{\Gamma(m+1)}\right)^{-1} = 0$$

Integrating the one remaining term (under the restriction (*j*) > −4 ∧ (((*h* + *j*) < −2 ∧ (*σ* − *τ*) ≤ 0 ∧ (*σ* + *τ*) ≥ 0) ∨ ((*h* + *j*) ≥ −2 ∧ (*σ* − *τ*) < 0 ∧ (*σ* + *τ*) > 0)) ∧ ((*τ*) > 0 ∨ ((*τ*) = 0 ∧ (*j* + 2*m*) < −6)) gives us our final result,

*B*(T, a, *σ*, h, *τ*, j, *β*, *γ*, *λ*, *α*, *γb*, ps, pt) = 6.811556×10−<sup>20</sup> 1.10114 × <sup>10</sup><sup>6</sup> √*T* <sup>√</sup>*<sup>π</sup>* (22*c*2) *gPI gRA* (34) × (*kBT*) 1 <sup>2</sup> *<sup>π</sup>*7/2*λ*<sup>3</sup> ((*d* + *λ*)<sup>2</sup> − *γ*2) ps (*<sup>f</sup>* + *<sup>λ</sup>*)<sup>2</sup> − *<sup>γ</sup>*<sup>2</sup> *b* pt × ∞ ∑ *m*=0 2*m*+25/2(*m* + 1)*a*−*m*<sup>−</sup> <sup>5</sup> <sup>2</sup> Γ(*m* + 1)(*kBT*)*<sup>m</sup>* 9*m*! 5 2 *<sup>m</sup>* Γ <sup>−</sup>*<sup>m</sup>* <sup>−</sup> <sup>3</sup> 2 Γ(2*m* + 4) × 8*k*<sup>3</sup> *<sup>B</sup>T*3<sup>Γ</sup> *m* + <sup>5</sup> <sup>2</sup> + 3 *<sup>a</sup>*<sup>3</sup> <sup>+</sup> 12*k*<sup>2</sup> *BT*2*γ*<sup>2</sup> 0Γ *m* + <sup>5</sup> <sup>2</sup> + 2 *a*2 + 6*kBTγ*<sup>4</sup> 0Γ *m* + <sup>5</sup> <sup>2</sup> + 1 *<sup>a</sup>* <sup>+</sup> *<sup>a</sup>*0*γ*<sup>6</sup> 0Γ *m* + 5 <sup>2</sup> <sup>+</sup> <sup>0</sup> <sup>×</sup> *<sup>τ</sup>*−*j*−4Γ(*<sup>j</sup>* <sup>+</sup> <sup>4</sup>)*σ*−*h*−2*m*−4Γ(*<sup>h</sup>* <sup>+</sup> <sup>2</sup>*<sup>m</sup>* <sup>+</sup> <sup>4</sup>) × <sup>4</sup>*F*<sup>3</sup> *j* <sup>2</sup> <sup>+</sup> 2, *<sup>j</sup>* 2 + 5 2 , <sup>−</sup>*<sup>m</sup>* <sup>−</sup> <sup>3</sup> 2 , −*m*; 5 2 , −*h* <sup>2</sup> <sup>−</sup> *<sup>m</sup>* <sup>−</sup> <sup>3</sup> 2 , −*h* <sup>2</sup> <sup>−</sup> *<sup>m</sup>* <sup>−</sup> 1; *<sup>σ</sup>*<sup>2</sup> *τ*2 

The rate coefficient is then a quadruple sum. First, we have a sum over the seven terms in *I*<sup>11</sup> + *I*<sup>21</sup> 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:

$$\begin{array}{llll} a(1) = 0.6878357996671101, \cdots & a(50) = 0.37080904876118337\\ \beta(1) = 1.2354854281591452, \cdots & \beta(50) = 1.1073078257847453\\ \gamma(1) = 0.012984468708341133, \cdots & \gamma(50) = 0.28320160279252\\ \varepsilon(1) = -2.534888772248287, \cdots & \varepsilon(50) = 0.02873436866901307 \end{array} \tag{35}$$

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.

#### **5. Comparison with Numerical Integration**

Table 1 shows the comparison of the present analytical rate coefficient *αRA* for attaching a positron to antihydrogen to form <sup>H</sup><sup>+</sup> 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 - 6 K, the rate coefficient is essentially linear and may be fit by *<sup>α</sup>RA* = 0.001050 × <sup>10</sup>−<sup>15</sup> cm3s−<sup>1</sup> *<sup>T</sup>* <sup>K</sup><sup>−</sup>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 × <sup>10</sup>−<sup>6</sup> + 1.203002 × <sup>10</sup>−<sup>9</sup> × <sup>10</sup>−<sup>15</sup> cm3s−<sup>1</sup> *T K*−*1*, which sum to a value 0.00105917 × <sup>10</sup>−<sup>15</sup> cm3s−<sup>1</sup> *T K*−*<sup>1</sup>* and the numerical integral for this same fifty-term wave function gives 0.00102625 × <sup>10</sup>−<sup>15</sup> cm3s−<sup>1</sup> *T K*−*1*, a 3% difference.

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.

**Table 1.** The rate coefficient *<sup>α</sup>RA* (10−<sup>15</sup> cm3s−1) for attaching a positron to antihydrogen to form <sup>H</sup><sup>+</sup> for a fifty-term wave function.


#### **6. Discussion and Concluding Remarks**

We found the analytic form for the rate coefficient *αRA* for attaching a positron to antihydrogen to form <sup>H</sup><sup>+</sup> 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 1*s* 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 *ne* = 1016 m−<sup>3</sup> 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 *<sup>α</sup>RA* = <sup>10</sup> × <sup>10</sup>−<sup>17</sup> cm3s−1, and if there is unit overlap between the positron plasma and the antihydrogen, then the reaction rate, the product *ne αRA*, at this temperature would be 10 × <sup>10</sup>−<sup>7</sup> <sup>s</sup>−<sup>1</sup> per antihydrogen atom. As the temperature falls to 10 K, the reaction rate falls in an essentially linear manner to be 1 × <sup>10</sup>−<sup>7</sup> <sup>s</sup>−<sup>1</sup> 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<sup>2</sup> 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.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares that there is no conflict of interest.

#### **References and Note**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
