Can this expression be simplified? (double integral over a sphere in $\mathbb{R}^3$.

Let's start by writing the integral as $$ I = \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} f(x+y) \delta(|x|-r) \delta(|y|-R)\, {\rm d}^3x {\rm d}^3y$$ where $\delta$ is Dirac's delta function. Changing the variables to $(z,y) = (x+y,y)$ we get \begin{align} I &= \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} f(z) \delta(|z-y|-r) \delta(|y|-R)\, {\rm d}^3y {\rm d}^3z \\ &= \int_{\mathbb{R}^3} f(z) \Big(\int_{\mathbb{R}^3}\delta(|z-y|-r) \delta(|y|-R)\, {\rm d}^3y \Big) {\rm d}^3z\end{align} For $0\le r \le R$, it can be calculated that To calculate $ \int_{\mathbb{R}^3}\delta(|z-y|-r) \delta(|y|-R)\, {\rm d}^3y $ that let us use spherical coordinates for $y$ where angle $\theta$ is measured from the direction of $z$, so that $$|y| = \rho $$ $$|z-y| = \sqrt{(|z|-\rho\cos\theta)^2 + (\rho\sin\theta)^2} = \sqrt{\rho^2 - 2\rho|z|\cos\theta + |z|^2}$$ We have then \begin{align} & \int_{\mathbb{R}^3}\delta(|z-y|-r) \delta(|y|-R)\, {\rm d}^3y = \\ &\quad = \int_0^\infty {\rm d}\rho\int_0^{\pi}{\rm d}\theta\int_0^{2\pi}{\rm d}\phi\,\rho^2\sin\theta \delta(\sqrt{\rho^2 - 2\rho|z|\cos\theta + |z|^2}-r) \delta(\rho-R) = \\ &\quad = 2\pi R^2 \int_0^{\pi}{\rm d}\theta \sin\theta \;\delta(\sqrt{R^2 - 2R|z|\cos\theta + |z|^2}-r)=^{u:=\sqrt{R^2 - 2R|z|\cos\theta + |z|^2}}\\ &\quad = 2\pi R^2 \int_{|R-|z||}^{R+|z|} {\rm d}u \frac{u}{R|z|} \delta(u-r)= \\ &\quad = \frac{2\pi R}{|z|} \int_{-\infty}^\infty {\rm d}u \,u \,\theta((R+|z|)-u)\theta(u - \big|R-|z|\big|)\delta(u-r)=\\ &\quad = \frac{2\pi R r}{|z|} \,\theta((R+|z|)-r)\theta(r - \big|R-|z|\big|)=\\ &\quad = \frac{2\pi R r}{|z|} \,\theta((R+|z|)-r)\theta(r - (R-|z|))\theta(r - (|z|-R))=\\ &\quad = \frac{2\pi R r}{|z|} \,\theta(|z|+(R-r))\theta(|z| - (R-r))\theta(R+r - |z|) = ^{\text{since }|z|>0, R-r>0}\\ &\quad = \frac{2\pi R r}{| z|} \theta\big(|z|-(R-r)\big)\theta\big(R+r-|z|\big) \end{align} where $\theta$ is Heaviside's theta function.

In the end you get $$ I = 2\pi R r \int_{|z|\in[R-r, R+r]} \frac{f(z)}{|z|} {\rm d}^3z$$

Here's a re-derivation of Adam Latosiński's formula without using distributions. Ultimately what needs to be done here is just a change of variables, but I'm going to use probabilistic language because I think it makes the whole thing easier to understand.

Observe that $$ E[f] = \frac{1}{4\pi R^2} \int_{|y|=R} \frac{1}{4\pi r^2}\int_{|x|=r} f(x+y)\,dS(x)\,dS(y) $$ is the expected value of $f(x+y)$ if we randomly choose vectors $x$ and $y$ for which $|x|=r$ and $|y|=R$. We want to change variables to $$ \rho=|x+y|\qquad\text{and}\qquad u = \frac{x+y}{\rho}. $$ Note that $u$ is a unit vector which can be chosen uniformly at random from the unit sphere in $\mathbb{R}^3$. Clearly $\rho \in [R-r,R+r]$, but the probability distribution of $\rho$ is not uniform on this interval, so we have to do some work to figure out how to sample $\rho$.

For this purpose, we may assume by symmetry that $y=(R,0,0)$ and $x$ is uniformly chosen on the sphere $|x|=r$. Now, it is a famous observation of Archimedes that each coordinate of the unit sphere in $\mathbb{R}^3$ is uniformly distributed on the interval $[-1,1]$. It follows that the first coordinate $t$ of $x$ is uniformly distributed on the interval $[-r,r]$. But $$ \rho = \sqrt{(R+t)^2 + (r^2 -t^2)} = \sqrt{R^2+r^2 + 2Rt}\tag*{(*)} $$ so we can sample $\rho$ by choosing $t$ uniformly at random from $[-r,r]$ and then using this formula.

We conclude that $$ E[f] = \frac{1}{2r}\int_{-r}^r \frac{1}{4\pi} \int_{|u|=1} f(\rho u)\,dS(u)\,dt $$ where $\rho$ is the formula (*) above. Note that the $t$ here actually represents the component of $x$ in the direction of $y$, i.e. $t = (x\cdot y)/R$. Combining this with our original formula for $E[f]$, we deduce that $$ \int_{|y|=R} \int_{|x|=r} f(x+y)\,dS(x)\,dS(y) = 2\pi rR^2 \int_{-r}^r \int_{|u|=1} f(\rho u)\,dS(u)\,dt $$ where $\rho$ is given by the formula above. Now observe that $$ d\rho = \frac{d}{dt}\Bigl[\sqrt{R^2+r^2 + 2Rt} \Bigr]\,dt = \frac{R}{\sqrt{R^2+r^2+2Rt}}\,dt = \frac{R}{\rho}\,dt. $$ Substituting $dt = (\rho/R)\,d\rho$, the integral becomes $$ 2\pi rR \int_{R-r}^{R+r} \int_{|u|=1} \rho\,f(\rho u)\,dS(u)\,d\rho $$ If we like, we can change this to a single volume integral over $z=\rho u$. Since $dV(z) = \rho^2\,dS(u)\,d\rho$ and $\rho=|z|$, the integral is $$ 2\pi rR \int_{|z|\in[R-r,R+r]} \frac{f(z)}{|z|}\,dV(z) $$ where $z$ represents $x+y$.