I'm considering a particle in a spatial domain $\Omega \subset \mathbb{R}^2$ which secretes some chemical substance with production rate $q_0$, and is subject to a force proportional to the gradient of its own emitted chemical field, $\nabla c$. I'm assuming that the chemical diffuses in the environment with rate $D$ and decays with constant rate $\alpha$. The concentration of the chemical substance then evolves according to
$$ c_t = q_0 \rho(\mathbf{r},t) - \alpha c + D\nabla^2 c \tag{1} $$
where $\rho(\mathbf{r},t) = \delta(\mathbf{r} - {\mathbf{r}_a}(t))$ represents the chemical point source. Similarly, ${\mathbf{r}_a}$ is the trajectory of the particle, determined by
$$ m\ddot{{\mathbf{r}}}_a = \gamma\nabla c + \mathbf{R}(t), \tag{2} $$
where $\gamma$ is a measure of the particle's response to the concentration gradient and $R$ is a force due to thermal noise. I would like to solve (1) and (2) numerically using a finite difference scheme. The problem, however, is in how to treat the Dirac Delta appearing in (1).
Let $C_{ij}^n \approx c(x_i,y_j,t_n)$ denote the numerical approximation to $c$ at the grid point $(x_i,y_j,t_n)$. Then, if I use an unconditionally-stable time-differencing procedure - let's say the Crank-Nicolson scheme - then one arrives at the following system of difference equations,
$$ \frac{C_{ij}^{n+1} - C_{ij}^{n}}{\delta t} = q_0\frac{1}{2}( \rho(\mathbf{r}_a^{n}) + \rho(\mathbf{r}_a^{n+1}) ) - \alpha\frac{1}{2}(C_{ij}^{n} + C_{ij}^{n+1}) + \frac{1}{2}D(\nabla_h^2C_{ij}^{n} + \nabla_h^2C_{ij}^{n+1} ) \tag{3} $$
where $\nabla_h^2$ is the discretized Laplacian operator. The problem is computing $\rho(\mathbf{r}_a^{n+1})$, as this requires me to have knowledge of the particle's position at the next time step. Is there any way around this? Can I omit this term somehow and still maintain unconditional stability in time?
As a side note, if anyone knows of papers addressing the numerical treatment of Dirac Deltas in time-dependent problems, please provide links.