I'm trying to prove the statement in the title in as simple a way as possible. It is Theorem 3.2.9 in Helein's book "Harmonic maps, conservation laws, and moving frames", although it is not proved there. The statement is as follows.
Suppose $\phi\in\mathbb{R}^m$ is a weak solution of $\Delta \phi = f \in \mathcal{H}^1$, where $\mathcal{H}^1$ is the standard Hardy space on $\mathbb{R}^m$. Then $ \Big\lVert\frac{\partial^2\phi}{\partial x^\alpha \partial x^\beta}\Big\rVert_{L^1(\mathbb{R}^m)} \le C\lVert f \rVert_{\mathcal{H}^1(\mathbb{R}^m)}. $
My idea is to use convolution with the kernel of the Laplacian, and then differentiate, estimate in $L^1$ and somehow interpolate between the $\mathcal H^1$ and $BMO$ norms. Then since the kernel of the Laplacian is in $BMO$, I am finished. However there are two problems with my proof: I don't know how to prove that one can interpolate a convolution between $\mathcal H^1$ and $BMO$ (atomic decomposition?) and I don't know how to prove that the kernel of the Laplacian is in $BMO$.
Does anyone have either a better way to prove this theorem, or a way to fix up my proof? Thanks!