Four-body treatment of the antihydrogen-positronium system: binding, structure, resonant states and collisions

We have developed a coupled-rearrangement-channel method allowing the rigorous non-adiabatic treatment of the multi-channel scattering problem for four particles. We present the study of the binding, resonant and collisional properties of the H̄−Ps\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\bar {H}-Ps$\end{document} system with the total angular momentum J = 0+ (singlet positronic configuration). The binding energy, the life-times of the resonant states and the collisional cross sections are calculated and discussed. We present the preliminary cross sections for the elastic and inelastic H̄−Ps\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\bar {H}-Ps$\end{document} scattering, notably for the excitation of Ps and for the rearrangement reaction producing the H̄+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\bar {H}^{+}$\end{document} ions.


Introduction
The experiment aiming to study the Gravitational Behavior of Antihydrogen at Rest (GBAR) involves the step whereby the antihydrogen atoms collide with positronium atoms to form the positronium hydride molecular ionsH + suitable for sympathetic cooling [1][2][3][4]. In view of that we have undertaken the investigation of theH P s system with a method capable to describe the multichannel aspects of theH − P s scattering, including the elastic and inelastic processes, particularly the rearrangement reaction leading to theH + production. TheH P s system has been studied before, both in works concentrating on its ground state [5][6][7][8] and in works treating the scattering [9][10][11][12]. In the present work we use the description that is suitable to study both these aspects. We study the system in terms of several sets of Jacobi coordinates that are natural for its sub-cluster structure; this facilitates analysis of the binding energy, formation and decay of resonances and the treatment of multichannel scattering.
To describe the structure of the theH − P s system we apply the variational approach based on the Gaussian Expansion Method (GEM) [13,14]. The variational manifold explicitly contains contributions from various rearrangement channels expressed in terms of the appropriate Jacobi coordinates. At the same time, the use of Jacobi coordinates allows an efficient and rigorous treatment of the scattering cross sections.
The binding energy is calculated and discussed in terms of the contributions from various rearrangement channels. Channel analysis helps to converge the binding energy and throws light on the structure and collisional properties of the system. The structure of the system is presented in terms of the correlation functions that portray the spatial distribution of the particles. This helps to understand the structure and dynamics of the system, in particular the change of theH + upon binding of an electron, and the coexistence of the atomic and molecular features ofH P s.
The resonant states and their life-times are calculated variationally using GEM in conjunction with the Complex Coordinate Method (CCM). The energies and widths of the resonances are analysed with respect to the channel composition of the resonant wave functions.
Based on the variational description of the four-body system, we apply the Coupled Rearrangement Channels Method (CRCM) to calculate the cross sections forH − P s collisions. The outer part of the totalH P s wave function is made to satisfy the appropriate scattering boundary conditions for the collisional fragments. This is facilitated by the use of Jacobi coordinates in the description of the multi-channel structure of the wave function, that expressly contains various rearrangement channels. The scattering matrix S and the cross sections are obtained from the coupled, non-local integro-differential equations that explicitly couple the collisional channels of interest.

General outline of the method
The scattering cross sections were calculated using the Coupled Rearrangement Channel Method [15,16]. In this method one constructs the total scattering wave function in terms of the inner part (that describes the internal, highly correlated part of the 4-body system) and the outer part (that describes the asymptotic motion of the various possible scattering fragments), see (1).
The inner part is expanded in terms of the square-integrable, 4-body basis functions v which are taken to be the eigenfunctions of the matrix eigenvalue problem for the total Hamiltonian of the system. Their construction is described in the next section. The important aspect of our approach is the use of several Jacobi sets of coordinates in the expansion of one and the same total wave function, and each function v contains contributions from several functional manifolds expressed in different sets of Jacobi coordinates (corresponding to different arrangement channels, see Fig. 1). The manifolds corresponding to various sets of coordinates are different from each other, even if the basis functions that span these manifolds are similar (they are all gaussians). In other words, various sets of Jacobi coordinates describe different parts of the configuration space of the system. Consequently, the problem of linear dependencies is alleviated and we do not suffer from the overcompletness.
The outer part is given in terms of the (a priori unknown) channel functions χ c (R c ) of the relative motion between various scattering fragments. The channel functions are coupled to each other by being matched to the inner part and becoming part of the total wave function . It is the asymptotic behavior of the channel functions that give the information about the Fig. 1 Rearrangement channels forH P s. We see the various arrangements of Jacobi coordinates, and the angular momenta for the wave functions in these coordinates. For the symmetrization purposes, we also include rearrangement channels with the positrons exchanged scattering cross sections. To this end, the total scattering wave function forH P s is written as In the above, all 3 sums are over 4-body basis functions, but these functions are differently constructed. The first sum contains the variational solutions of the 4-body problem, obtained in the basis that contains functions from as many Jacobi arrangement channels as needed for the good description of the binding energy. In the present calculation we have used 7 Jacobi arrangement channels (out of 15 possible for the system) in the single calculation that couples all of them. The remaining summations are over 4-body functions that describe the asymptotic character of those open scattering channels that one wishes to consider and couple. Here, function φ (n) c describes the n-body fragment of the scattering channel c; r c , q c and R c are the 3 Jacobi coordinates in that channel.
The second sum is over channels that asymptotically contain two 2-body fragments. In the present case, these fragments areH and P s and are described by functions φ (2) c (r c ), φ (2) c (q c ), respectively. The relative motion of the centers of mass of these fragments is described by the function χ c in terms of the third Jacobi coordinate R c that connects the mass centers. The three Jacobi coordinates for the channel functions in this group look graphically like the letter H (see Fig. 1), we abbreviate the summation as c ∈ H .
The third sum is over channels that asymptotically contain a three-body fragment (H + or P s + ). The wave functions for these fragments (φ (3) c (r c , q c )) are obtained by separate 3body calculations. The motion of the 4th particle is described by the function χ c in terms of the third Jacobi coordinate R c that connects the lone particle to the CM of the 3-body fragment. The three Jacobi coordinates for the functions in this group look graphically like the letter K (see Fig. 1), we abbreviate the summation as c ∈ K.
As will be seen in the following, expansion coefficients b v and functions χ c in (1) are determined simultaneously via the self-consistent, integro-differential procedure.
We first describe the variational procedure that generates the expansion functions v for the inner part of the scattering wave function.

Binding energy and structure ofHPs
To obtain the inner part of the total wave function, we solve the Schrödinger equation forH P s by means of the variational approach, using the Gaussian Expansion Method in Jacobi coordinates [13,14]. The expansion functions v are taken to be the solutions of the eigenvalue problem with respect to the total Hamiltonian H projected onto the subspace P spanned by the Gaussian basis. The variational expansion is given by where ϕ c i (r c μ ) are the basis functions for the Jacobi coordinate r c μ , μ = 1, 2, 3, in the channel c. The Schrödinger equation takes form and is solved as the matrix eigenvalue problem. The basis functions ϕ c i in (2) are either primitive Gaussians or the atomic orbitals formed out of these Gaussians (these two sets are connected via an unitary transformation). The sum over channels c ≡ C; l 1 , l 2 , l 3 runs over the various Jacobi sets of coordinates (we call them rearrangement channels and number them with C, see Fig. 1) and various associated triple configurations of the angular momenta l c 1 , l c 2 , l c 3 compatible with the coupling to the total angular momentum J .
The bracket around the triple products in (2) means that they are coupled to the total angular momenta J, M and appropriately symmetrized. In the present paper we consider the J = 0 + symmetry, i.e. the states with the total (orbital) angular momentum 0 and (natural) parity +1, with the positrons in the singlet (spin 0) configuration. The spin is not treated explicitly, rather we symmetrize the spacial wave functions to be space symmetric with respect to the permutation of the spacial coordinates of the two positrons. The latter symmetrization is achieved by "doubling" the Jacobi rearrangement channels, i.e. by including their companions obtained by permutation P of the two positrons (except for those arrangements that contain the e + 1 − e + 2 coordinate and need not be symmetrized, since ϕ(r c ).
There are alltogether 15 Jacobi rearrangement channels for theH P s system. Our studies of the convergence of binding energy revealed that the most important ones are channels 13,1,10 and 15 illustrated on Fig. 1. The channel numbers C = 13, 1, 10, 15 conform to the previous literature [13]. Thus we use 7 arrangements: arrangements 13, 1,10; their companions obtained by permutation of the two positrons; and arrangement 15 (which contains the e + 1 − e + 2 coordinate and does not need symmetrization). Our calculations show that the most important contribution to the binding energy comes from the molecularH − P s channel (C = 13). Next in importance is the atomic channel with the electron orbiting theH + atomic ion core (C = 1). Next comes the positronium ion channel, P s + +p (C = 10). After that comes the repulsion (or total break up) channel where the CM of two mutually repelling positive particles is connected to the CM of two mutually repelling negative particles (C = 15).
The convergence of the binding energy is shown in Table 1 and illustrated on Fig. 2. The achieved accuracy is on the order of few micro-Hartrees. Last column shows which channels are included. The notation + 13/1/10/15-(220,022,202) means that we added the essential combinations with lλL = 220, 022, 202 for channels 13, 1, 10 and 15. The adopted proton mass is m p = 1836.152 672 47

Calculations of the scattering length and the scattering cross sections for antihydrogen -positronium collisions
The outer part of the total wave function describes the asymptotic behavior of the system (cf. (1)). It is a sum over channel functions that satisfy the correct scattering boundary conditions. The sum runs over open physical channels c, i.e. over various possible fragmentations of the system and the angular momenta of the relative motion of the fragments. We include both the so called H channels (fragmentation into two two-body fragments) and K channels (fragmentation into a three-body fragment and a single particle). Hence φ (2) in (1) are the wave functions of isolated 2 and 3 body fragments in channel c; r c and q c are the internal coordinates of these fragments. The relative motion of the fragments in various scattering channels is given by the functions χ c (R c ), with R c being the Jacobi coordinate for this motion. If these fragments are neutral (e.g. the atomic fragments) the wave functions of relative motion satisfy  Table 1. Our experience shows that it is more important to increase the number of channels than to merely increase the number of the Gaussian functions spanning each channel with u l (k c R c ) being the Ricatti-Hankel functions [17], υ c i , υ c being the velocities of the fragments in the initial and final channel, respectively, and S c,c i being the elements of the scattering S-matrix.
If the fragments are charged (as in partitioning into a 3-body system and a lone particle), the function of the relative motion satisfies where u (C+) l (k c R c ) is the outgoing spherical Coulomb function [18]. Actually, the S-matrix elements in the above equations depend not only on the initial and final angular momenta given by c i and c, but also on the total angular momentum of the system J , its parity, and the intermediate couplings of the angular momenta of the collision fragments. More about this is coming later in the text.
It is important to emphasize that functions χ c (R c ) are not requested to vanish in the inner region (except for R c = 0). In that way they can contribute to the description of the inner region spanned by functions v , and at the same time facilitate a smooth transition between the inner and outer regions.
Solving the Schrödinger equation (3) with expressed as in (1) becomes equivalent to the determination of expansion coefficients b v and the channel functions χ c . The strategy has been discussed in ref. [19] for the case ofH H . The solution of the eigenvalue problem is converted into the solution of the system of non-local, coupled integro-differential equations. We show the general form of these equations in the Appendix, postponing the derivations and description of numerical procedures to the more comprehensive publication. The objective is then to solve these equations as to simultaneously and self-consistently determine the outer part (χ c and thereby S c,c ) and the inner part (b v ). The asymptotic functions χ c are determined by numerical integration using the Compact Finite Difference Method (CFDM) that takes into account the boundary conditions.
The expansion coefficients b v can be expressed in terms of channel functions χ c [19]. This suggests the use of iterative procedure. However, the attempts to solve the problem iteratively did not succeed, we were not able to achieve convergence. Instead, we have solved the problem by casting the numerical integration into algebraic equations which were solved as the matrix problem, determining the values of the functions χ c on the grid, in form of one vector for all functions χ c . Subsequently, the scattering matrix elements S c,c i , and the cross sections are determined, i.e. χ c −→ S c,c i −→ σ c,c i . To construct the cross sections we need a more subtle notation for the channel angular momenta. We replace the symbol c by a pair α, λ where α defines the angular momenta of the collisional fragments l and L (see Table 2 and Fig. 1), and λ is the angular momentum of their relative motion. We also need to take into account that the three angular momenta that characterize the channels c can be coupled in different ways. We adopt the coupling scheme where the momenta of the scattering fragments (l and L) in the channel state α are added first and the result ( ) is added to the momentum of the relative motion of the fragments λ, to give the total angular momentum J .
The cross sections, for scattering from the initial state α i in the partial wave λ i into a final state state α in a partial wave λ f , are given by The above partial cross sections need to be calculated for each λ i = 0, 1, 2..., for the associated values of |λ i − i | ≤ J ≤ (λ i + i ) and for the allowed values of λ f , compatible with J and f . The above partial cross sections can be summed up to give the cross section for scattering from the initial state α i to a final state α Finally, the state to state cross sections can be summed up over the final states to give the total cross section where θ is the Heaviside step function and E α is the collision energy threshold above which the reaction to state α becomes open.
In the present work we consider the scattering ofH (1s) and P s(1s) below theH (n = 2) threshold and for this case the notation in (6) -(8) becomes significantly simpler, because Fig. 3 The energy spectrum of theH P s system of the J = 0 + symmetry with positrons in the singlet (spin 0) configuration. There are two energy scales: the total 4-body energy, and the energy relative to the lowest threshold. The orange dots indicate the resonances. There exists only one bound state, just below the threshold for the dissociation into antihydrogen and positronium atoms in their ground states. The various branches of the continuous spectrum are indicated by blue lines, originating at the thresholds. We see the thresholds for positronium excitation, which accumulate at the 3-body fragmentation threshold i = l i = L i = 0, J = λ i and f = L f . Therefore, for this case the equation (7) becomes In the above equation, J is equivalent to the relative angular momentum in theH (1s) + P s(1s) scattering (J = λ i ) and thus σ λ f ,λ i α,α i (J = λ i ; E) can be regarded as a partial cross section for scattering from the initial state α i = α 1 in J -wave to the final state α in λ f -wave.
In the following we present the cross sections for the s-wave scattering betweenH (1s) and P s(1s), thus setting α i = α 1 (see Table 2) and λ i = 0. This situation conforms entirely to the J = 0 symmetry (the wave functions possess the J = 0 + , singlet positronic configuration symmetry). The excitation of the positronium states P s (n f , L f ) is possible but only in the final partial waves with the angular momenta satisfying λ f = L f . The considered range of collisional energies exceeds the threshold forH + formation. This means that we need to solve the coupled channel problem for all seven opened channels, and determine the 7 × 7 scattering matrix S. The energy spectrum ofH P s showing the scattering channels is shown on Fig. 3. The channel labeling for the included coupled channels is given in Table 2. (Please note that the number of open scattering channels only accidentally coincides (and should not be confused) with the 7 arrangement channels included in the calculation of the inner part of the wave function.) The results for the s-wave elastic scattering from the one channel calculation are presented in Fig. 4. The one channel calculation is exact below the first excitation threshold for P s(n = 2). The construction of the inner part of the scattering wave function is similar to that for the bound state calculation (shown in Table 1). The latter was constructed using 41 648 fourbody basis functions from 16 channels in 7 sets of Jacobi coordinates, including angular momentum configurations with l c μ ≤ 2. In the scattering calculations we use the similar size expansion (44 496 basis functions) to describe the closed part of the wave function (the first summation in (1)). However the number of channels is reduced, so that the closed part can be considered to be a restriction of the one used in the calculation of the binding energy of H P s. On the other hand we have extended the size of the manifolds along the coordinates describing the relative motion of the scattering fragments, especially by additional longrange functions along theH + −e coordinate as to describe the expected oscillation structure of the resonant states.
The basis functions for the 2-body (atomic) and 3-body (ionic) fragments (appearing in the second and third summation in (1)) are obtained by separate 2 and 3 body calculations. These functions are compact but of good quality. Thus for instance the fragment of antihydrogen ionH + is expanded in 736 basis functions, resulting in the binding energy -0.527 442 a.u. (and thus differing from the best available reference value only by 4 micro Hartrees).
Suma sumarum, at the present stage of calculations, our scattering expansions are less elaborate than the one for the bound state. This is why the accuracy of the binding energy of H P s is on the order of 5 ppm (parts per million) whereas the non-unitarity of the scattering matrix for the here reported scattering calculations is on the order of 5% (at worst). 1 Still, the quality of the 4-body inner part is such that already the one channel calculation is able to discern the resonant features both below and above the first excitation threshold, and at the same time give the very competitive scattering length. We mention in passing that the good quality of the inner part of the scattering wave function means that it alone can be used to get the first estimation of the elastic cross section. This is shown on Fig. 4 Table 3 which include the results for the elastic cross section obtained by the simple Projection Method combined with the Effective Range Theory (PM/ERT) [20]. The importance of the good description of the inner part (reflected e.g. in the resulting binding energy) on the scattering length has been noted before, e.g. in ref. [10].

and in
To judge the accuracy of the obtained cross sections at low energies we report the scattering length a, related to the low energy limit of the elastic scattering according to σ 0,0 11 (E → 0) = 4πa 2 . Comparison of our scattering length to other calculations is given in Table 3.
Since the scattering length obtained from a variational calculation is tacitly assumed to be an upper bound for this quantity, our result is seen to be highly accurate. The accuracy of the cross sections at higher energies is estimated from the residual non-unitarity of the scattering matrix S(E).
On Fig. 6 we illustrate the importance of channel coupling and the influence of inelastic scattering on the elastic one. Even though the elastic scattering dominates, the inelastic one has a substantial effect. We see that the pure one-channel elastic scattering (black broken On Fig. 7 we show the cumulative scattering cross sections resulting from the 3-channels calculation, i.e. the sum of σ 0,0 11 , σ 0,0 21 and σ 1,0 31 . The sum of these 3 cross sections does not correspond to the pure partial wave scattering. Rather, we include those partial waves that are compatible with the total angular momentum J = 0 for the entire system, since at the moment we consider collisions ofH (1s) with P s(1s) with λ i = 0 for the initial relative motion (low energy collisions). Hence, the outcoming productsH (1s) and P s(2p) are kept in the λ f = 1 partial wave as to ensure conservation of the total angular momentum.
The three different shades of Fig. 7 show the separate contributions and the cumulative effect of the presented cross sections. Each of them is obtained from the common coupled 3channels calculation. We notice that σ 0,0 21 , σ 1,0 31 start to contribute above the lowest inelastic channel (that of P s(n = 2) excitation) and that the behavior of the cross sections across the thresholds for P s(n = 2) and P s(n = 3) excitations is markedly affected by the nearby resonances.
To analyse the resonant structure of the cross section we have located the resonances through an independent calculation, using the complex coordinate method and solving the complex eigenvalue problem. As seen from the upper panel of Fig. 7 the complex eigenvalues that determine the positions and the life-times of the resonances correspond very well to the resonant features seen in scattering cross sections obtained by the CRC method.
The results of the 7-channels calculation, coupling theH + + e channel and all open channels below, is presented on Fig. 8. In the energy interval shown in Fig. 8  satisfies the unitarity condition with the accuracy of 5% at worst. We notice that, once thē H + production sets in, it becomes the strongest inelastic process with the cross section on the order of 1 a.u. .

Summary
We have presented the first round of the 4-body, Coupled Rearrangement Channels calculations for theH P s system. The simultaneous use of several sets of Jacobi coordinates in the same calculation allows the systematic analysis of the binding, structure and scattering properties of the system. It also alleviates the problem of linear dependencies (at the present level of the description we have not noticed it at all). Using the 7 most essential rearrangement channels we have achieved the binding energy of the ground state E b = 0.788 867 5 a.u. which indicates the accuracy on the order of micro-Hartrees.
We have then used the 4-body solutions of the eigenvalue problem forH P s as the basis for the expansion of the inner part of the total scattering wave function. The asymptotic behavior of that wave function, rendering the scattering matrix elements, was determined through solving the coupled, integro-differential equations. We have solved these equations including all open scattering channels for energies above the threshold for the rearrangement reactionH + + e. We have calculated the s-wave cross sections for the elastic scattering, the inelastic scattering ending with excitation of the positronium, and the rearrangement reaction ending with theH + ion. We have seen that, notwithstanding some resonant features, above the threshold for theH + production this reaction is the strongest inelastic process, with the cross section on the order of 1 a.u. The resonant features, attributed to the presence of charged fragments in theH + + e channel, are seen already in the purely elastic scattering, but, as seen from our coupled-channel procedure, are then substantially modulated by the presence of inelastic processes. The accuracy of our cross-section calculations can be inferred from the residual nonunitarity of the S-matrix. In the presented energy range of the scattering calculations, the unitarity of the S-matrix has been satisfied to within 5% (at worst).
For the low energy scattering, the indication of the accuracy may be obtained from the consideration of the low energy limit of the elastic scattering, which becomes σ 0,0 11 (E = 0) = 231.38 a 2 0 . The corresponding scattering length is a = 4.291 a 0 which is currently the lowest informal upper bound for this quantity.