Let $\Omega \subset \mathbb R^d$ be a bounded domain, $d>2$. Let $f \in C^\infty(\Omega)$. If $u \in L^2$ is a distributional solution of $\Delta u = f$ in $\Omega$ then $u \in C^\infty(\Omega)$ and it is a classical solution. This can be easily proved: we can explicitly specify one of the solutions of $\Delta u = f$ of class $С^\infty$, it is $u' = K*f$, where $K$ is the Newton kernel. Next we can easily show that any solution of $\Delta u = 0$ is $C^\infty$ from what we can derive that $u = u' + (u-u') \in C^\infty$. More precisely this procedure is described in these lecture notes, p.8, thm. 2.
The situation is more complicated in the case when $f \in H^k$, where $H^k$ is the Sobolev space, $k \geq 0$. In this case I expect to show that if $u$ is a distributional solution of $\Delta u = f$ then $u \in H^{k+2}$ and $\Delta u = f$ a.e. As above we represent $u = u'+(u-u')$ where $u'$ is some particular distributional solution of $\Delta u = f$. The problem is to show that we can specify such particular solution in the class $H^{k+2}$. Formally we can take convolution of distributions $K*f$, but I'm not sure that $K*f \in H^{k+2}$ since $K$ is not in $L^2$ for $d \geq 4$. Is there a way to generalize the above approach to prove regularity in this Sobolev case?
P.S. I know about other ways to prove regularity, e.g. [Bers, Schechter], but I'm interested if the approach for $C^\infty$ functions is easily extendable.