I'm interested in Poisson's equation $-\Delta u=f$ set in the whole space $R^d$ (let's say $d\geq 3$ for simplicity) when $f$ has very little integrability, specifically $f\in L^{1+\varepsilon}$ for small $\varepsilon$. I'm also interested in $f\in L^1\log L^1$ or the limit case $f\in L^1$ but not in general measure data, so let me focus here on $L^{1+\varepsilon}$ for simplicity.
Since the solution is defined up to harmonic functions one usually uses the convention that $$ u(x)=\int\limits_{R^d} G(x-y)f(y)dy $$ where $G(x-y)=\frac{1}{|x-y|^{d-2}}$ is the Poisson's Kernel (up to normalizing constants). Unless I'm mistaken my hypothesis $f\in L^{1+\varepsilon}$ implies by the Hardy-Littlewood-Sobolev lemma that $u\in L^{d/(d-2)+\varepsilon_1}$ and $\nabla u\in L^{d/(d-1)+\varepsilon_2}$ for some small $\varepsilon_1,\varepsilon_1>0$ related to $\varepsilon$.
Now assume, for given $f\in L^{1+\varepsilon}$, that there is a distribution $u'$ solving $-\Delta u'=f$ in the sense of distributions such that $u'\in L^{d/(d-2)+\varepsilon_1}$ but also $\nabla u'\in L^2$. Unless I missed something this readily implies $f\in H^{-1}$.
Question: Can I guarantee that $u'=G*f$, that is $u'=u$?
This is of course a uniqueness statement. I am aware that uniqueness is not automatically guaranteed for such $f$'s. I've heard of Stampacchias's adjoint method and Boccardo-Gallouet's approximation method, but I'm not familiar with any of those (it seems that the former guarantees uniqueness but only works for linear equations, while the latter does not ensure uniqueness but works for nonlinear problems). Since I have here $\nabla u'\in L^2$ (hence $f\in H^{-1}$) I'm hoping to avoid this kind of things, is that hopeless?
This question is also related to this other post
Thank you in advance for your time and thoughts
Edit: just added the link
