Context
I recently read a paper (resume) on quantum chemistry and an incredible analytical approximation has been made for computing ground state of Hydrogen to Boron atoms. I used this result to write a paper (equation (43)), expressing quantum mechanics in 3D space and obtained the same results. For Helium atom, we arrive at a radial Newton-Schrodinger type equation and we know that we have almost the solution using Hylleraas result in Hartree-Fock model (roughly an exponential times a polynomial).
Helium differential equation
According to Rioux's result and the equivalent treatment in my paper, we have for Helium wave function the following Schrödinger equation in atomic units (I pass the details on Laplace expansion):
$$ \left[ - \frac{\Delta_r}{2} - \frac{Z}{r} + \frac{1}{2} \left( \frac{1}{r} \int_0^r r'^2 R(r')^2 d r' + \int_{r}^{+ \infty} r' R(r')^2 d r' \right) - \epsilon \right] R(r) = 0 $$
with $Z$ the number of proton (2 for Helium) and $\epsilon < 0$ the eigenvalue to compute. This is equivalent to the system:
$$ \left[ - \frac{\Delta_r}{2} - \frac{Z}{r} + \frac{1}{2} V(r) - \epsilon \right] R(r) = 0 $$ $$ \Delta_r V(r) = - R(r)^2 $$
So this is a radial Schrodinger-Newton problem with a Coulomb potential $- Z/r$. We know that it is non-linear and there is no closed form but it is only in the case $Z = 0$. There is no boundary condition, it is an eigenvalue problem with an infinite set of solutions, exactly like hydrogen atom. The solution is normalized as:
$$ \int_0^{+ \infty} r^2 R(r)^2 d r = 1 $$
Question
Do you have an idea to solve this non-linear differential eigenvalue equation ?
A hint on the possible solution
We know that in Hartree-Fock theory (multi-electronic wave function) the wave function is something like an exponential times a polynomial and we obtained a really good precision of $0.07$ % on the smallest eigenvalue (with Z = 2 which is $\epsilon \simeq - 2,90/2$, because there is two electron, the total energy in this formulation is $2 \epsilon \simeq - 2.90$). So as guesses the function that can suit the first non-linear ODE could be something like:
$$ R(r) =^? \frac{1}{N} \sum_{n = 0}^{+ \infty} a_n r^n e^{-b_n r} $$
or:
$$ R(r) =^? \frac{1}{N} \sum_{n = 0}^{+ \infty} (1 + a_n r)^n e^{-b_n r} $$
with $a_n$ real numbers, $b_n >0 $ real numbers, $N$ normalization constant. Do you have an idea to find the solution of the first equation ? It is extremelly difficult at first sight but the extreme accuracy ($0.07$ % of error) of the exponential function variational fit let think that maybe $R(r)$ is composed of decreasing exponential... and probably we can use a series method to solve it.
Few tracks to observe for polynomials specialists
Laguerre polynomials: In quantum mechanics, for hydrogen atom or harmonic oscillator, the solution is analogous: Laguerre polynomials, the argument is different ($r$ for hydrogen atom, $r^2$ for harmonic oscillator). Anyway the solution very probably admit a Laguerre decomposition (see Laguerre transform).
Orthogonality: the solutions ($L^2$ integrable function) should be an orthonormal set. The Rodrigues formula (which exist for a huge range of polynomials like Laguerre, Hermite and others) using derivative ensure an orthogonality relation and should maybe used to find the solution.
Another track using radial Ladder operator
Ladder operator formulation can be used to compute ground state of Hydrogen atom, Harmonic oscillator (in 1D and 3D). We can maybe use this formulation to compute Helium ground states, probably by a reccursive formula. See this paper to compute the radial operator for the Hydrogen atom.