I was solving this problem
Find the solution $u(x,y,t)$ of the problem $$ \begin{cases} u_t=D\nabla^2u \\[5pt] u(x,y=0,t) = 0\\[5pt] u(x,y,0) = \phi(x,y) \end{cases} $$ where $$ \phi(x,y) = \begin{cases} c &(x-x_0)^2+(y-y_0)^2\le R^2 \\[5pt] 0 & (x-x_0)^2+(y-y_0)^2\gt R^2 \end{cases}$$ then find an approximate solution in the region $$(x-x_0)^2+(y-y_0)^2\gg R^2$$ $\underline{\textbf{Hint:}}$ Note that the Dirac delta function is given by the limit when $R\rightarrow \infty$ of the function $$\frac{1}{\pi R^2}\chi_R (x,y)$$ where $$\chi_R(x,y) = \begin{cases} 1& x^2+y^2\le R \\ 0&x^2+y^2\gt R\end{cases}$$
and I got stuck evaluating the following integral: taking the convolution between the heat kernel and the initial condition $$u(x,y,t) = \frac{1}{4\pi Dt}\iint_{\mathbb R^2} e^{-\frac{|\underline{x}-\underline{x}'|^2}{4Dt}}u(\underline{x}',0)\,\mathrm d^2\underline x' = \frac{c}{4\pi Dt}\iint\limits_{\Omega} e^{-\frac{|\underline{x}-\underline{x}'|^2}{4Dt}}\,\mathrm d^2\underline x'$$ where $\Omega$ is the region limited by the disk centered in $(x_0, y_0)$ of radius $R$.
I've tried going to polar coordinates, which seems to me the right choice, but I haven't got very far: $$\begin{cases} x' = x_0+r\cos\theta \\ y'=y_0+r\sin\theta\end{cases}$$ which gives me $$u(x,y,t) = \frac{c}{4\pi Dt}e^{(x^2+y^2)+(x_0^2+y_0^2)-2(x_0+y_0)}\int\limits_0^R\int\limits_0^{2\pi}r^2\sin(\theta) e^{-2r[\cos\theta(1-x_0)+\sin\theta(1-y_0)]+r^2}\,\mathrm dr\,\mathrm d\theta$$ here's where I'm stuck.
How do I proceed? Are my calculations correct?
Sorry for the late correction, I forgot the Jacobian when changing the variables.
WLOG, suppose $x_0=0,\ y_0>0,\ D=\frac12,\ c\pi R^2=1$.
By Green's function of the heat equation $u(x,y\le0,t)=0$.
For $y>0$, the $u$ in the upper half plane is equivalent to the upper half plane solution of whole plane solution of from the original source and its negative image $-\Omega$ reflected across the line $y=0$, since the latter satisfies the PDE and its zero boundary condition on $y=0$. \begin{align} u(x,y>0,t) &= \frac c{\sqrt{2\pi t}}\Big(\int_{\Omega}dudv\, e^{-\frac{(x-u)^2+(y-v)^2}{2t}}-\int_{-\Omega}dudv\, e^{-\frac{(x-u)^2+(y-v)^2}{2t}}\Big) \\ &\underset{x^2+y^2\to\infty}\sim\ \frac1{\sqrt{2\pi t}}\Big( e^{-\frac{x^2+(y-y_0)^2}{2t}}-e^{-\frac{x^2+(y+y_0)^2}{2t}}\Big) \\ &=\sqrt{\frac2{\pi t}}\,e^{-\frac{x^2+y^2+y_0^2}{2t}}\sinh\Big(\frac{yy_0}t\Big). \end{align} We have treated the source as one point, as at large distance the exact shape $\Omega$ of the source impacts little on the final approximation.
The disk is cut into two parts by $y=0$. We still treat the upper and lower half of the plane separately. This time the total "energy" of positive source is a fraction proportional to its area above $y=0$. $u$ no longer vanishes in the lower half plane but is evaluated similarly to that in the upper half plane.