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)$?

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.