# Does the $p^+p^-e^-$ system have bound states?

This is batted a bit off the wall, so bear with me. The hydrogen molecular ion, $\mathrm{H}_2^+$, is the simplest three-body system in molecular physics, and in the Born-Oppenheimer approximation its electronic Schrödinger equation is separable (though it's unclear whether it's actually exactly solvable).

Now, an old paper of A.H. Wilson, (The Ionised Hydrogen Molecule, Proc. Roy. Soc. Lond. Ser. A, Math. Phys. 118 no. 780, pp. 635-647 (1928)) shows that the electronic Schrödinger equation $$\left[-\frac12\nabla^2-\frac{Z_1}{\|\mathbf r-\mathbf R_1\|}-\frac{Z_2}{\|\mathbf r-\mathbf R_2\|}\right]\psi(\mathbf r)=E\psi(\mathbf r),$$ when separated in the elliptic coordinates for internuclear distance $2c$ \begin{align} \xi&=\frac{1}{2c}\left(\|\mathbf r-\mathbf R_1\|+\|\mathbf r-\mathbf R_2\|\right) \quad\text{and}\\ \eta&=\frac{1}{2c}\left(\|\mathbf r-\mathbf R_1\|-\|\mathbf r-\mathbf R_2\|\right) \end{align} as $\psi(\xi,\eta,\phi) = X(\xi) Y(\eta)e^{im\phi}$, for $\phi$ the angle about the internuclear axis, reads \begin{align} \frac{\mathrm d}{\mathrm d\xi}\left[(1-\xi^2) \frac{\mathrm d X}{\mathrm d\xi}\right] +\left[\lambda^2\xi^2-\kappa \xi -\frac{m^2}{1-\xi^2}+\mu\right]X&=0 \\ \frac{\mathrm d}{\mathrm d\eta}\left[(1-\eta^2) \frac{\mathrm d Y}{\mathrm d\eta}\right] +\left[\lambda^2\eta^2-\kappa' \eta -\frac{m^2}{1-\eta^2}+\mu\right]Y&=0, \end{align} where $\lambda^2=-2c^2E$ is the eigenvalue, $\mu$ is a separation constant which must be found (and which annoyingly appears in both equations), and \begin{align} \kappa & = 2c(Z_1+Z_2),\\ \kappa' & = 2c(Z_1-Z_2). \end{align}

OK, so that's a lot of set-up, but this is what I want to ask: I notice that the Schrödinger equation above is weirdly symmetric with respect to the combinations $\boldsymbol{Z_1+Z_2}$ and $\boldsymbol{Z_1-Z_2}$ of the two nuclear charges, so there are probably some weird symmetries of the electronic wavefunctions if you flip the charge of one of the protons - i.e. switch it for an antiproton, and pretend that the two are not going to annihilate within some meaningful timespan.

With that in mind, then: has this been explored in the literature? Does the resulting $p^+p^-e^-$ system, with one proton, one antiproton, and one electron, have bound states? (Intuitively, you'd expect the electron to bunch around the proton with most of the probability on the far side from the antiproton, thereby partially shielding the $p^+$-$p^-$ attraction, but hopefully this would be insufficient to completely kill that attraction.) If so, what exactly are the relationships between its electronic eigenstates and those of the hydrogen molecular ion?

• Hat-tip to DavePhD for pointing out that paper in a now-deleted answer. Mar 27 '17 at 18:34
• Well, if the protonium were stable, it would still be neutral, so there would be no overall Coulomb potential for electronic bound states to form. Would a (neutron +electron) have a bound state? Mar 27 '17 at 18:57
• @annav That's not really true - the internal details of the potential matter. The electron can sit (mostly) around the Coulomb well of the potential and in the process, hopefully, not turn away the antiproton. Analogous systems are possible - for example, the hydrogen anion, $\mathrm H^-$, has one (though in this case only one) bound state. One can't just conclude "oh, there's just too much negative charge, this'll fly apart". Mar 27 '17 at 19:03
• The electromagnetic protonium potential well solution will have higher than grounds state open levels, so you are arguing that an electron could be captured on that level, the way it is captured in H-. The mass is crucial in the Bohr radius and thus how tight the hydrogen type energy levels are phys.spbu.ru/content/File/Library/studentlectures/schlippe/… . The protonium has practically 3 orders of magnitude difference from the hydrogen Bohr radius, so mass is important in the charge distribution. Mar 28 '17 at 4:51
• This paper suggests a bound state with binding energy of $0.08793\,\mathrm{eV}$. Apr 3 '17 at 9:35

I have worked previously on high accuracy three-particle atomic and molecular calculations using a series solution method first developed by C. L. Pekeris back in 1958. Myself and others developed a high accuracy code that is able to calculate ground state and excited S states of any three-particle system for any combinations of charges and masses (note, The Born-Oppenheimer approximation is never assumed, all particles are treated on an equal footing) and the $r_{12}$ component is explicitly included in the non-relativistic Schroedinger equation and solved accordingly (i.e. finite mass effects and electron-electron correlation effects are included). These calculations are very close to exact in the non-relativistic limit as you can get without actually solving the N-particle Schroedinger equation exactly.

These calculations are optimised using various Non-Linear variational Parameters (NLP's) to find the minimum along the energy surface. From solving the generalised eigenvalue problem we are left with high accuracy wavefunctions for these systems which can then be used to calculate physical properties.

I ran the system $p^+p^-e^-$ through the code to see what would happen. After optimisation along the energy surface the calculation resulted in a 1S ground state energy of E = -459.038 168 472 $E_h$ (Hartrees) using a 1078$\times$1078 basis set size, which is converging to the two-body protonium energy.

Upon looking at expectation values calculated using the output eigenvector, it reveals that this system is not bound (within this theory). The expectation value of $r_{12}$, ($\left\langle r_{12}\right\rangle$) is 1522608.953 $a_0$ showing that the proton and electron have an enormous spatial separation (showing the electron has detached) and the expectation value of the Dirac delta function for the separation between proton and antiproton, ($\left\langle \delta(r_{1})\right\rangle$) is 246312486.976998. i.e. there is an incredibly high probability of the proton and antiproton coalescing.

These calculations are without application of relativistic corrections from the Breit-Pauli theory or QED corrections, but still offers insight as to the bound state of this system. You may find this paper here an interesting read which talks about the critical mass required to bind three-particles. The reference trail within may be of interest.

Edit my other comment included in this answer:

The ground state energy of protonium is the energy of the lowest continuum threshold of the $p^+p^−e^−$ system. If you calculate the expectation value of the dipole moment resulting from these equal but opposite charges via $p^+,p^-$ is tiny. The dipole polarisability of the $p^+p^−$ atom, that determines the magnitude of the attractive potential between the electron and ground state protonium asymptotically is very small. It is even smaller than the proton-muon system, so there is no chance of a bound state. Furthermore if we analyse the ternary diagram of bound states (shaded region shows bound states) between the three particles $p^+p^-e^-$ it sits near the vertex in the unbound region (taken from this paper):

• I'm not sure that breaking out of Born-Oppenheimer helps much, but if analyzed carefully enough then it can't hurt, either. I'm not sure I understand your results (how can you reach a negative-energy minimum but not have a bound state? Negative with respect to what? Is your electron running into the edge of the box defined by your basis? If so, is your calculation converged?) but your energy is remarkably close to that mentioned here. Jun 29 '17 at 14:23
• That said, if your calculation suggests the $p^+p^-$ system 'coalesces', I'd be inclined to distrust it, since if you take the electron away it will just form a bound protonium, with mean separation $a_0\times m_e/2m_p>0$. Moreover, it is perfectly possible (in principle) to have a loosely-bound electron orbiting a neutral protonium (in the same way that you can loosely bind a second electron to hydrogen), so a large proton-electron distance is not that surprising. Jun 29 '17 at 14:31
• This calculation uses a Coulomb potential between all particles and thus one of the limiting factors will be that there is no explicit inclusion of the strong force between proton and anti-proton which is the dominant force in the protonium system. I have not seen that paper before, and it is interesting that this quick calculation matches the value quite well.
– Yeti
Jun 29 '17 at 14:38
• The problem is posed in electrostatic terms so the solution shouldn't involve the strong force (I'm interested in the problem as posed, not in real protons and antiprotons ;-)). The point was that there is an electrostatic bound state between the proton and the antiproton (identical to hydrogen but with different masses and therefore different scales on all quantities), and your basis needs to include functions that cover the relevant length scales (much tighter than hydrogen) to be accurate. Jun 29 '17 at 14:53
• Moreover it seems that this is what you've found and the ${\sim}459E_{\rm H}$ energy is the protonium ground state, at $\displaystyle \frac{m_p/2}{m_e}$ times the hydrogen ground-state energy, with the $1/2$ coming from the reduced mass of protonium. (Also, apologies, I misquoted above - I meant to link here.) Jun 29 '17 at 14:55

This is certainly not a complete answer, but atleast addresses the existence of a bound state.

My claim is that there always exists a bound state. If the ground state energy were calculated as a function of internuclear distance, then taking the anti proton to infinity would take us back to the hydrogen atom, guaranteeing the existence of a bound state. Now, if at some point, the solution qualitatively changes, giving no bound state solutions, there definitely exists a length scale where the ground state energy vanishes (This follows from the continuity of the Ground state energy as a function of internuclear distance). We would need to establish a preferred length scale where such a "transition" occurs. While the system does have a character length scale $a_0$, there is no reason why electronic length scales have anything to do with this transition length scale. Therefore, the absence of an internuclear length scale makes the transition from bound ground state to free ground state unlikely. This is something we would like to establish.

The Strategy: We establish an upper bound on the ground state energy using the variational principle. We then show that this upper bound is always negative, therefor establishing a negitive energy bound state solution for all internuclear distance.

The variational principle says that the upper bound on the ground state is given by: $$\langle \psi |H|\psi \rangle \geq E_g$$ Where $|\psi \rangle$ is an arbitrarily chosen state vector in Hilbert space, $H$ is the Hamiltonian of the system, and $E_g$ is the ground state energy.The larger the projection of $|\psi \rangle$ onto the the Eigen space of ground state, the stronger the bound.

Explicit calculation: We choose the function $$\psi(r)=\frac{1}{\sqrt{\pi}}e^{-r}$$ as our arbitrarily chosen wave function, since this lies entirely in the ground state Eigen space in the infinite internuclear distance limit, hence promising a somewhat strong bound. We use the translational symmetry of the Hamiltonian to place the proton at origin. We then place the antiproton at a distance $c$ from the origin along the $z$ axis. The Hamiltonian may then be rewritten in spherical coordinates as: $$H=-\frac{1}{2} \nabla ^2-\frac{1}{r}+\frac{1}{\sqrt{r^2+c^2-2rc\cos \theta}}$$ Its easy to see that $$\frac{1}{2}\nabla^2\psi(r)=\left( \frac{e^{-r}}{2}-\frac{e^{-r}}{r} \right)\frac{1}{\sqrt{\pi}}$$ Hence the Inner product is given by: $$\langle \psi |H|\psi \rangle=\frac{1}{\pi} \int_{r=0}^\infty \int_{\theta =0}^{\pi} \int_{\phi =0}^{2 \pi}\left( \frac{e^{-2r}}{2}-\frac{2e^{-2r}}{r}+\frac{e^{-2r}}{\sqrt{r^2+c^2-2rc\cos \theta}} \right) r^2 \sin \theta d\theta d\phi dr$$ This integral is not hard to do and Mathematica says the value of this integral is: $$\frac{2-3c-2(1+c)e^{-2c}}{2c}$$ The numerator is a monotonically decreasing function, that vanishes for $c=0$, hence establishing a negative upper bound on the ground state energy, for all $c>0$

EDIT: I have messed up the algebra here, and the integral really looks like: $$\langle \psi |H|\psi \rangle=\frac{1}{\pi} \int_{r=0}^\infty \int_{\theta =0}^{\pi} \int_{\phi =0}^{2 \pi}\left( -\frac{e^{-2r}}{2}+\frac{e^{-2r}}{\sqrt{r^2+c^2-2rc\cos \theta}} \right) r^2 \sin \theta d\theta d\phi dr$$ This integral is really: $$-\frac{-2+c+2(1+c)e^{-2c}}{2c}$$ Which is only monotonic and negative in $c \geq 1.87 a_0$ However, for small $c$ you can still do better, by "guessing" the ground state wave function intuitively. A simple guess, that minimises the probability for $\theta \in (0,\frac{\pi}{2})$ and maximises the probability $\theta \in (\frac{\pi}{2},\pi)$ is: $$\psi(r,\theta)=N*(1-\cos \theta)r e^{-r}$$ This when integrated gives something like: $$-\frac{3 - 5 c + 4 c^2 + c^3 + (-1 + c) (1 + 2 (-1 + c) c) (3 + 2 (-1 + c) c) e^{2c}}{c^3}$$ Which is negative in the range $0.75 a_0 \leq c\leq 2a_0$

Therefore, for $c>0.75a_0$ we have bound state solutions. What happens below this is an interesting question to ask. I still believe that there will be bound state solutions for all $c>0$, However, a proof looks very elusive. The problem is potentially richer than it originally appears!