While reading the article Finite element approximations for Schrödinger equations with applications to electronic structure computations, I don't understand part of the proof of regularity. For the Schrödinger eigenvalue problem,

\begin{cases}
 -\Delta u+Vu=\lambda u, \ \text{in}\ \Omega \\
 \Vert u\Vert_{0,\Omega}=1
\end{cases}
where $\Omega$ is a bounded domain in $\mathbb{R}^3$, $V =\frac{1}{|x|}$.

Define
$$a(u,v)=\int_\Omega \nabla u\nabla v+Vuv,\ u,v\in H_0^1(\Omega)$$
A number $\lambda$ is called an eigenvalue of the form $a(\cdot, \cdot)$ relative to the form $(\cdot, \cdot)$ if there is a
nonzero vector $u\in H_0^1(\Omega)$, called an associated eigenfunction, satisfying
$a(u,v) = \lambda(u,v)$ for all $v \in H_0^1(\Omega)$.

Now, if $(\lambda, u) \in \mathbb{R}\times H_0^1(\Omega)$ and $Vu \in L^2(\Omega)$, can we deduce from the regularity of the Poisson equation that $u=(-\Delta )^{-1} (Vu+\lambda u) \in H^2(\Omega)$?

<cite authors="Gong, Xingao; Shen, Lihua; Zhang, Dier; Zhou, Aihui">_Gong, Xingao; Shen, Lihua; Zhang, Dier; Zhou, Aihui_, Finite element approximations for Schrödinger equations with applications to electronic structure computations, J. Comput. Math. 26, No. 3, 310-323 (2008). [ZBL1174.65047](https://zbmath.org/?q=an:1174.65047).</cite>