The tilte says it all. Given $f\in L^1(R^d)$ (let me restrict to dimension $d\geq 3$ for convenience), what is the optimal $L^p$ regularity for solutions to $$ -\Delta u=f\hspace{3cm}(1)? $$ I'm of course aware that solutions to (1) are only defined up to harmonic functions, so I'm implicitly speaking here of the "integral" solutions defined by means of the Newton potential $$ u(x)=\int\limits_{R^d}\frac{1}{|x-y|^{d-2}}f(y)\,\mathrm{d}y\hspace{3cm}(2) $$ up to normalizing constants.
When I discussed that with one of my colleagues he claimed that the map $f\mapsto u$ is continuous from $L^1$ to $L^p$ for all $p\in [1,\frac d {d-1})$, and that the statement also holds for Radon measures $\mu$ instead of $f\in L^1$ and replacing $f(y)\,\mathrm{d}y$ by $d\mu(y)$ in (2). Unfortunately he couldn't point me to a precise reference, and a quick websearch returns quasilions of results for compactly supported $\mu$ but nothing really relevant to me for $f\in L^1$ (a priori supported in the whole space).
I apologize if the question is trivial, but I'm not familiar with potential theory. As I only need a "black box" result for a specific problem I would greatly appreciate if you could point me directly to a precise statement (e.g. theorem x.x.x page y). Please feel free to close or migrate to SE if you deem it appropriate.
Thank you in advance.
Edit: I know of course about the Hardy-Littlewood-Sobolev inequality, but the latter only gives $L^p$ information on $u$ if $f\in L^q$ with $q>1$, and weak Lebesgue $L^{p,r}$ if $q=1$ so this is off-topic here since I'm only interested in the "classical" Lebesgue $L^p$ regularity of $u$.