integration - Multiple integral over a disc

I would need some help on this integration problem:

$$I=\int_0^{2\pi}\int_0^{R}\int_0^{2\pi}\int_0^{R}\exp(-a\ r_{12}) \ r_1 \ r_2 \,\mathrm{d}r_1\,\mathrm{d}\phi_1\,\mathrm{d}r_2\,\mathrm{d}\phi_2$$

Here, $r_1, r_2, \phi_1$ and $\phi_2$ are the polar coordinates for the integral on a disc with radius R. And $r_{12}$ is the distance between the coordinates, which can be calculated by

$$r_{12}=\sqrt{r_1^2+r_2^2-2r_1r_2\cos(\phi_1-\phi_2)}.$$

Below a plot of the problem for $a=1$, $R=1$ and fixed coordinates of the first point at $(0,0)$:

Answer

This integral has a closed form in terms of special functions. First, start with the integral in the form
$$ I(\alpha, R) = \int_{D(0,R)}da\int_{D(0,R)}db e^{-\alpha|a-b|} $$
with parameters $\alpha$, $R$, and differentiate with respect to $R$:
$$
\begin{aligned}
\frac{dI}{dR}(\alpha, R) &= \lim_{\epsilon\to0}

\frac{2}{\epsilon}\int_{D(0,R)}da\int_{D(0,R+\epsilon)\setminus D(0,R)}db\, e^{-\alpha|a-b|}
\\&=
4\pi R\int_{D(0,R)}da\, e^{-\alpha|a - R e_1|}
\\&= 4\pi R^3\int_{D(0,1)}da\, e^{-\alpha R|a-e_1|}.
\end{aligned}$$
Here $e_1=(1,0)$ is a unit vector on the unit circle, and the integral simplified because $\int da\,e^{-\alpha|a-b|}$ is independent of where $b$ is on the circle $\partial D(0,R)$ of radius $R$. The factor of $2$ came from the two occurrences of $D(0,R)$ that are symmetric, and $2\pi R$ is the area of the strip $D(0,R+\epsilon)\setminus D(0,R)$ divided by $\epsilon$.

To evaluate this integral, write $a$ in polar coordinates, but with the origin at the point $(1,0)$, so that the unit disk has the form
$$\pi/2 \leq \phi \leq 3\pi/2, \qquad 0 \leq r \leq -2\cos\phi = \rho(\phi). $$
Then the integral for $\frac{dI}{dR}$ becomes

$$ \begin{aligned}
\frac{dI}{dR} &= 4\pi R^3\int_{\pi/2}^{3\pi/2}d\phi\int_0^{-2\cos\phi}r\,dr\, e^{-\alpha R r}
\\&= 4\pi R^3\int_{\pi/2}^{3\pi/2}d\phi\frac{1-(1+\alpha R\rho)e^{-\alpha R\rho}}{\alpha^2 R^2}.
\end{aligned} $$
The integral over $\phi$ can be done with the help of computer algebra:
$$ \frac{dI}{dR} =
\frac{4\pi^2 R}{\alpha^2}\big(1-I_0(2\alpha R)+2\alpha R I_1(2\alpha R)-2\alpha R \mathbf{L}_{-1}(2\alpha R)+\mathbf{L}_0(2\alpha R)\big), $$
where $I$ are the modified Bessel functions, and $\mathbf{L}$ are the modified Struve functions.

Then the integral $I(\alpha,R) = \int \frac{dI}{dR}$ can also be done by a computer in terms of Struve and hypergeometric functions, which yields

$$ I(\alpha,R) = \frac{2\pi R}{\alpha^3}\big(
-2+\alpha \pi R-\alpha \pi R \,{}_0\mathbf{F}_1(2,\alpha^2R^2)+2 \alpha^3 \pi R^3 \,{}_0\mathbf{F}_1(3,\alpha^2 R^2)-2 \alpha \pi R \mathbf{L}_{-2}(2 a R)+\pi \mathbf{L}_1(2\alpha R)\big) $$
where ${}_0\mathbf{F}_1(a;z) = {}_0F_1(a;z)/\Gamma(a)$ is the regularized hypergeometric function.

This is the Mathematica expression for the above formula:

1/a^3 2 \[Pi] R (-2 + a \[Pi] R - 
  a \[Pi] R Hypergeometric0F1Regularized[2, a^2 R^2] + 
  2 a^3 \[Pi] R^3 Hypergeometric0F1Regularized[3, a^2 R^2] - 
  2 a \[Pi] R StruveL[-2, 2 a R] + \[Pi] StruveL[1, 2 a R])

EDIT Differentiating the integral. Let $D=D(0,R)$, $C=D(0,R+\epsilon)\setminus D(0,R)$. Then
$$
\left(\int_{D}+\int_C\,da\right)\left(\int_D+\int_C\,db\right) - \int_Dda\int_C db
= \int_Cda\int_D db + \int_Dda\int_Cdb + \int_C\int_C\,da\,db.$$
The term $\int_C\int_C$ is on the order of $O(\epsilon^2)$ because the area of $C$ is $2\pi R\epsilon$, so it may be ignored.

Because the function $e^{-\alpha|a-b|}$ is symmetric in $a,b$,
$$\int_Cda\int_D db + \int_Dda\int_Cdb = 2\int_Dda\int_Cdb. $$

Also, the integral $\int_D e^{-\alpha|a-b|}\,da$ depends only on the distance between $b$ and the origin, so picking a special value of $b=Re_1$ gives
$$ 2\int_D da\int_C db\,e^{-\alpha|a-b|} = 2\mathop{\mathrm{area}}(C)\int_Dda\,e^{-\alpha|a-Re_1|}. $$