Let $u\not\equiv0$ be a non-negative solution of the following equation: $$ -\Delta u+u=|u|^{p-2}u,\quad u\in H^1(\mathbb{R}^N), $$ where $N\geq3$ and $2<p<2^*=2N/(N-2)$. My goal is to prove that $u\in C^2(\mathbb{R}^N)$. The following is my attempt.
First, since $-\Delta u=au$, where $a=|u|^{p-2}-1\in L_{\text{loc}}^{N/2}(\mathbb{R}^N)$, then by Brezis-Kato theorem, we know $u \in L_{\text{loc}}^p(\mathbb{R}^N)$ for all $1 \leq p< \infty$ and then $u \in W_{\text{loc}}^{2,p}(\mathbb{R}^N)$ by standard $W^{2,p}$ estimate. Then we know $u\in C^1(\mathbb{R}^N)$. But how to deduce the $C^2$ regularity? I looked up many books, but I didn't find a similar theorem.
Could anyone please show my a related reference or give some hints of this problem? I'll appreciated it a lot!