# Heuristics behind Dirac delta function in Master equation for probability?

I'm reading this paper [Phys. Rev. Lett. 106, 160601 (2011)] and it studies simple diffusion where a particle stochastically resets to its initial position $x_0$ at a constant rate $r$. As you can see, equation (1) is the master equation for $p(x,t|x_0)$, the probability that the particle is at $x$ at time $t$, having begun from $x_0$:

$$\frac{\partial p(x,t|x_0)}{\partial t}=D\frac{\partial^2p(x,t|x_0)}{\partial x^2}-rp(x,t|x_0)+r\delta(x-x_0),$$

I understand the origin of the LHS and of the first term of the RHS, they come from the simple diffusion process. But what about the second and third terms of the RHS? I think the second one has to do with a negative flux out each $x$ (due to stochastic resetting, as the paper says), and the third one has to do with the positive flux into $x_0$, but is there a (heuristic) way to understand a little further the use of a Dirac delta function for this last term? For example, what if the resetting occurs to a set of points? How this modify the master equation?

Here's how I would come to some intuition for it. I would think about the rate of "probability flow" into a region by integrating the equation over a region in space. For now, let's suppose that no diffusion occurs at all, since that is more complicated (although directly doable and understandable). Then $$\int_a^b dx\frac{\partial p(x,t|x_0)}{\partial t}=\int_a^bdx\left(-rp(x,t|x_0)+r\delta(x-x_0)\right),$$ which we can write as $$\frac{\partial}{\partial t}\int_a^b dx~p(x,t|x_0) = -r\int_a^bdx~p(x,t|x_0)+r\int_a^bdx~\delta(x-x_0),$$ which we can further write as $$\frac{\partial}{\partial t}P(a\leq x\leq b,t|x_0) = -rP(a\leq x\leq b,t|x_0)+r\int_a^bdx~\delta(x-x_0),$$ where $P(a\leq x\leq b,t|x_0)$ is the probability of there being a particle between $a$ and $b$ at time $t$ given that it started at $x_0$.
Now, if $x_0$ is not contained in the interval $[a,b]$, then this becomes $$\frac{\partial}{\partial t}P(a\leq x\leq b,t|x_0) = -rP(a\leq x\leq b,t|x_0),$$ which just says that the probability of the particle being in the interval decays with time exponentially, which is exactly their resetting process, represented by the exponential decay of being at a particular spot other than $x_0$.
On the other hand, if $x_0$ is contained in the interval, then the equation becomes $$\frac{\partial}{\partial t}P(a\leq x\leq b,t|x_0) = r(1-P(a\leq x\leq b,t|x_0)).$$ There are two components to the probability flow: there's the one we had before, that says that we are "losing probability" because there is an exponential process to jump away from $x$, but we are gaining probability because particles are jumping to $x_0$.
If we consider a centered interval about $x_0$, then $$\frac{\partial}{\partial t}P(x_0-\epsilon\leq x\leq x_0+\epsilon,t|x_0) = r(1-P(x_0-\epsilon\leq x\leq x_0+\epsilon,t|x_0)),$$ and if we narrow in on $x_0$ by taking $\epsilon$ to be small, then this is approximately given by $$\frac{\partial}{\partial t}P(x_0-\epsilon\leq x\leq x_0+\epsilon,t|x_0) \approx r,$$ which just states that the rate at which particles are resetting is exactly $r$.