Let $\Phi:\mathbb{R}^d\to \mathbb{R}$ be the fundamental solution to the Laplace equation, i.e the unique function $\phi$ such that $\Delta \phi = \delta_0$ in the sense of distributions. The solution $u$ to the equation $\Delta u = f$ in the full space $\mathbb{R}^d$ is given by the convolution $\Phi * f$, where $f$ is a function such that the above convolution is well defined (say, $f$ is compactly supported and bounded).
I am tasked with showing that $\|D^2u\|_{L^2(\mathbb{R}^d)} = \sum_{i,j = 1}^{n}\left \| \frac{\partial^2 u }{\partial x_i\partial x_j}\right\|_{L^2(\mathbb{R}^d)} = \|f\|_{L^2(\mathbb{R}^d)}$.
What is the idea here? The most obvious thing to do seems to be integrate by parts,i.e do something like: $$ \|D^2u\|_{L^{2}(\mathbb{R}^d)} = \int_{\mathbb{R}^d} u_{ij}u_{ij} = - \int_{\mathbb{R}^d} u_{iij}u_j = \int_{\mathbb{R}^d} u_{ii}u_{jj} = \int_\mathbb{R^d} (\Delta u)^2 = \|f\|_{L^2(\mathbb{R}^d)} $$ With the left hand side obeying the summation convention. But this doesn't seem right. The above method only works for $C^3$ solutions, and moreover, these solutions need to be compactly supported. We certainly don't have either guarantee. What's the idea here?