I know that
if $f\colon \mathbb{R}^3 \to \mathbb{R}$ is integrable on any measurable, according to the usual $n$-dimensional Lebesgue measure $\mu_y$, and bounded subset of $\mathbb{R}^3$, and if $g \in C^k(\mathbb{R}^3)$ is compactly supported, then the function$$h:x\mapsto \int_{\mathbb{R}^3} f(y-x)g(y)\,d\mu_y$$belongs to $C^k(\mathbb{R}^3)$, and its partial derivatives of order $\leqslant k$ are given by
$$D^{\alpha} h(x) = \int_{\mathbb{R}^3} f(y-x)D^{\alpha} g(y)\,d\mu_y$$
whose proof, whose author I thank again, is found here. An analogous result holds if $g\in C^k(\mathbb{R}^3)$ is bounded and $f\in L^1(\mathbb{R}^3)$, as this other proof, whose author I also thank again, shows.
I would say that, if $$f(z):z\mapsto \frac{1}{\|z\|}$$ and the partial derivatives $D^\alpha g$, with $0\le\sum_{i=1}^n\alpha_i\le k$, of $g$ decay more rapidly than $\|y\|^{-(2+\varepsilon)}$, $\varepsilon>0$, i.e. if $$\exists C>0\,\exists\varepsilon>0:\forall y\in\mathbb{R}^3\quad\|x\|^{2+\varepsilon}|D^\alpha g(y)|< C$$(a condition weaker than belonging to the Schwartz space, but satisfied by Schwartz functions), then, for all $x\in\mathbb{R}^3$, the function defined by $$y\mapsto f(y-x)D^\alpha g(y) d\mu_y$$is Lebesgue summable on $\mathbb{R}^3$. Can we extend the lemma above to this $f$, i.e. is it true and, if it is, how can it be proved (if possible by using elementary real analysis tools like basic theorems about Lebesgue integrations) that $$D_x^\alpha\int_{\mathbb{R}^3} \frac{g(y)}{\|y-x\|} d\mu_y=\int_{\mathbb{R}^3}\frac{D_y^\alpha g(y)}{\|y-x\|} d\mu_y?$$I am not sure whether that holds and how can it be proved because the $\chi_{L-x_0}$ used in the first linked proof cannot be used this time. I heartily thank any answerer.
What you are looking at here is (up to a constant factor) the so called Newtonian potential of $g$, which is known to be of class $C^1$ if $g$ is bounded and measurable and of class $C^2$ if $f$ is bounded and Hölder continuous for some exponent $\alpha \le 1$. These results are discussed in most books about (elliptic) partial differential equations, since the Newtonian potential is the fundamental solution to the Laplace equation.
A proof of the two above mentioned regularity results is found, e.g., in Gilbarg and Trudinger, Elliptic Partial Differential Equations of second order, chapter 4. The basic idea is to first derive results about the growth of the derivatives of the Newton potential and then use a suitable mollifier to uniformly approximate $D_{x_i}\int f(y)\Gamma(x-y)$ where $\Gamma $ denote the Newtonian potential (which, in general for $n\ge 3$ is given by $$\Gamma(x,y):=C(n) |x-y|^{n-2}$$ (in two dimesions it's the $C(2) \log(|x-y|)$)
The derivative turns out to be $$D_i\int f(y)\Gamma(x-y)dy = \int f(y) D_i \Gamma(x-y) dy$$ (for a proof I'd refer to the book of Gilbarg and Trudinger, or to this link)
(Formulae for higher order derivatives are then easily derived by induction)