Suppose we have $N$ spherical surfaces with radius distributed i.i.d. according to the PDF $$f_{d}(r) = \frac{3r^2}{R_{\text{max}}^3 - R_{\text{min}}^3}$$ on $[R_{\text{min}},R_{\text{max}}]$. On each surface, choose a point $d_i$ uniformly at random.
Choose a particular $d_k$. I am trying to evaluate the expectation $$\mathbb{E}[\left |\sum_{i=1}^N \exp \left ( 2\pi j(d_{i}-d_{k}) \right ) \right |^2]\tag{1}$$ in the case $R_{min} \leq d_{i} \leq R_{max}$. There, (1) simplifies to $$ \mathbb{E}\left [\left(\sum_{i=1}^N \cos \left (\frac{2\pi(d_i-d_k)}{\lambda} \right )\right)^2 \right ] + \mathbb{E}\left [\left(\sum_{i=1}^N \sin \left (\frac{2\pi(d_i-d_k)}{\lambda} \right )\right)^2 \right]^2\tag{2} $$ Note that each expectation is very symmetric: it is the square of a sum over all $N$ points chosen from the distribution, a set from which $d_k$ was chosen. Can we apply facts about equidistant points on a circle to finish?
I don't think splitting up into sines and cosines is usable this early. Instead, note that the magnitude of an absolute value is given by multiplying by the complex conjugate: $$|z|^2=zz^*$$ Now expand the square: \begin{align*} \left|\sum_i{\exp{\!\left(\frac{2\pi j}{\lambda}(d_i-d_k)\right)}}\right|^2&=\sum_i{\exp{\!\left(\frac{2\pi j}{\lambda}(d_i-d_k)+\frac{2\pi j}{\lambda}(d_k-d_i)\right)}}+\\ &\phantom{{}={}}\quad\quad\sum_{i\neq l}{\exp{\!\left(\frac{2\pi j}{\lambda}(d_i-d_k)+\frac{2\pi j}{\lambda}(d_k-d_l)\right)}} \\ &=\sum_i{1}+\sum_{i\neq l}{\exp{\!\left(\frac{2\pi j}{\lambda}(d_i-d_l)\right)}} \end{align*} Interestingly, $d_k$ has vanished: the result is independent of which $d_i$ you choose.
At this point, one can refactorize, but then independent random variables get all intermixed. So we shan't.
By linearity of expectation, it suffices to find the expectation of each summand. The first sum is just $N$; the second has identically-distributed terms. So we seek \begin{align*} \mathbb{E}\left[\left|\sum_i{\exp{\!\left(\frac{2\pi j}{\lambda}(d_i-d_k)\right)}}\right|^2\right]&=N+2(N-1)\mathbb{E}\left[\exp{\!\left(\frac{2\pi j}{\lambda}(d_1-d_2)\right)}\right] \\ &=N+2(N-1)\mathbb{E}\left[e^{\frac{2\pi jd_1}{\lambda}}\right]\mathbb{E}\left[e^{-\frac{2\pi jd_2}{\lambda}}\right] \end{align*} by independence.
To evaluate these expectations, we use the law of the unconscious statistician (note that the radius of $d_1$ and its angle are independent): $$\mathbb{E}\left[e^{\frac{2\pi jd_1}{\lambda}}\right]=\int_{R_{\text{min}}}^{R_{\text{max}}}{\int_0^{2\pi}{\exp{\!\left(\frac{2\pi j}{\lambda}\cdot re^{\theta j}\right)}\left(\frac{d\theta}{2\pi}\right)}\left(f_d(r)\,dr\right)}$$
I will interpret the inner integral as a contour integral. Let $C$ denote the counterclockwise circular contour of radius $\frac{2\pi r}{\lambda}$ starting at $1$. Then \begin{align*} \int_0^{2\pi}{\exp{\!\left(\frac{2\pi j}{\lambda}\cdot re^{\theta j}\right)}\left(\frac{d\theta}{2\pi}\right)}&=\int_0^{2\pi}{\exp{\!\left(\frac{2\pi r}{\lambda}e^{\theta j}\cdot j\right)}\cdot\frac{2\pi r}{\lambda}e^{-\theta j}\,d\left(\frac{2\pi r}{\lambda}e^{\theta j}\right)\cdot\frac{1}{2\pi}} \\ &=\frac{1}{2\pi}\oint_C{e^{jz}\,\frac{dz}{z}} \\ &=1 \end{align*} by Cauchy's integral formula.
Thus the double integral is just $$\mathbb{E}\left[e^{\frac{2\pi jd_1}{\lambda}}\right]=\int_{R_{\text{min}}}^{R_{\text{max}}}{1\cdot f_d(r)\,dr}=1$$ The same argument shows $$\mathbb{E}\left[e^{-\frac{2\pi jd_2}{\lambda}}\right]=1$$
Substituting these values back in, $$\mathbb{E}\left[\left|\sum_i{\exp{\!\left(\frac{2\pi j}{\lambda}(d_i-d_k)\right)}}\right|^2\right]=N+2(N-1)\mathbb{E}\left[e^{\frac{2\pi jd_1}{\lambda}}\right]\mathbb{E}\left[e^{-\frac{2\pi jd_2}{\lambda}}\right]=3N-2$$
Worth noting: this solution implicitly assumes $\{d_i\}_i$ are chosen from 2D spheres (i.e., circles). But the reason the contour integral comes out so nicely is just that $e^z$ is harmonic, so the average over a sphere is the same as the value at the center. That is, we should expect the same result in higher dimensions.