Can this expression be simplified? (double integral over a sphere in $\mathbb{R}^3$. Suppose $f:\mathbb{R}^3\to\mathbb{R}$ is an integrable function (I'm open to further restricting $f$, e.g. $f$ can be assumed to be smooth).
For some constants $0<r\leq R$, is there a way to simplify the expression
$$
\int_{\lvert y \rvert = R}\int_{\lvert x\rvert =r} f(x+y) \,\mathrm{d}S(x)\mathrm{d}S(y)
$$
to an expression not involving surface integrals?
 A: 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$$
A: 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$.
A: If you work in Fourier space you can write $f(x+y)=\int e^{ik(x+y)}\hat{f}_k dk$ and you can separate the integral to $\int dk \hat{f}_k \int e^{ikx} dS(x) \int e^{iky}dS(y) $.
