Reference request: optimal $L^p$ regularity for solutions to $-\Delta u=f$ with $f\in L^1(R^d)$ 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$.
 A: The following result is well known. It is Theorem 5.1 in [1]. It is proved by the method of duality solutions sue to Stampacchia. 

Theorem. Let $\mu$ be a signed Borel measure with the finite total variation in a bounded open set $\Omega\subset\mathbb{R}^{n}$. Then
  the Dirichlet problem
         $$
         \left\{\begin{array}{ccc}
                 \Delta u&=& \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \mu \\
                        u&\in& W^{1,1}_{0}(\Omega)
                \end{array}
         \right.
         $$ has a solution with the property that $u\in W_{0}^{1,p}(\Omega)$, $\Vert \nabla u\Vert_{L^{p}(\Omega)} \leq
 C\Vert\mu\Vert(\Omega)$ for all $1\leq p<n/(n-1)$. If the boundary of
  $\Omega$ is sufficently smooth, then the solution is unique.

Clearly, the result applies to $f\in L^1$, because such a function defines a measure of finite total variation $d\mu = f(x)\, dx$.
As was pointed put by  YangMills in his comment, you can find
a short proof of the result in the case when $d\mu=f(x)\, dx$, $f\in L^1$ in Lemma 14 in https://arxiv.org/abs/0809.2172.
In the case of solutions in $\mathbb{R}^n$ you obtain regularity $u\in W^{1,p}_{\rm loc}$ for all $1<p<n/(n-1)$.
[1] W. Littman, G. Stampacchia, H. F. Weinberger,
Regular points for elliptic equations with discontinuous coefficients. 
Ann. Scuola Norm. Sup. Pisa (3) 17 1963 43–77. 
