Distribution of collision times for 2-D diffusion in a area of radius $R$ with concentric disc of radius $r<R$ and absorbing boundaries

Question

I'm trying to derive the time probability distribution function of collision times $p_R(t)$ and $p_r(t)$ for a time-zero uniform random concentration of $N$ diffusing particles with diffusion coefficient $D$ for a circular area of radius $R$ and concentric disc $r<R$. Once a particle hits a surface at $r$ or $R$ it is annihilated/absorbed. $p_\rho(t)dt$ is the probability of a particle hitting and being absorbed by the circle of radius $\rho$ (where $\rho$ is outer $R$ or inner $r$) between time $t$ and $t+dt$. I've written a Monte Carlo engine to simulate but would like a closed-form solution.

Answer 1

You want $\dot\rho = D\nabla^2 \rho$, where $\rho$ is the density; given the radial symmetry, $\nabla^2 \rho =\partial_r^2 \rho + \frac{1}{r}\partial_r \rho$. Decaying eigenfunctions (such that $\dot\rho = -D k^2\rho$) will be solutions to $\partial_r^2 \rho + \frac{1}{r}\partial_r \rho + k^2 \rho = 0$, which are the ordinary Bessel functions of the first and second kinds, $J_0(kr)$ and $Y_0(kr)$. (The Bessel functions of the second kind are singular at the origin, so we can only include them because we're working in $0 < R_{\text{min}} \le r \le R_{\text{max}}$.) You need linear combinations that vanish at $R_{\text{min}}$ and $R_{\text{max}}$: $$a_k J_0(k R_{\text{min}})+b_k Y_0(k R_{\text{min}}) = a_k J_0(k R_{\text{max}})+b_k Y_0(k R_{\text{max}}) = 0.$$ Only for (infinitely many) special values of $k$ will the necessary condition be met for a nontrivial solution: $$J_0(kR_{\text{min}}) Y_0(kR_{\text{max}}) - Y_0(kR_{\text{min}}) J_0(kR_{\text{max}})=0.$$ The functions $a_k J_0(kr) + b_k Y_0(kr)$ at these special values will form an orthogonal basis; in particular, you'll find that $$ \int_{R_{\text{min}}}^{R_{\text{max}}}(a_k J_0(kr) + b_k Y_0(kr))(a_l J_0(lr) + b_l Y_0(lr))\cdot r dr \propto \delta_{k,l}. $$ You can use this to write your initial condition as $\rho(r, t=0)=\sum C_k (a_k J_0(kr)+b_k Y_0(kr))$ for some coefficients $C_k$, and then $\rho(r,t)=\sum C_k(a_k J_0(kr) + b_k Y_0(kr))\exp(-Dk^2t)$ for all later times.