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

Question

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?

Answer 1

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

Answer 2

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!