Consider a variable $M(x,t)$ driven by space-time Gaussian white noise: $$ \partial_t M(x,t) = -k M + \xi(x,t).$$ The noise has mean $\langle \xi(x,t) \rangle = 0$ and correlation function $\langle \xi(x,t)\xi(y,s) \rangle = 2D\delta(x-y)\delta(t-s)$.
Supposing the initial condition is $M(x,0) = 0$, one can integrate for the solution given a particular realization of the noise: $$ M(x,t) = \int_0^t dt_1 e^{-k(t-t_1)}\xi(x,t_1).$$
My goal is to calculate the second moment $\langle M(x,t)^2 \rangle.$ It's (apparently) easy to square the noise and use the correlation function to calculate the second moment: $$\langle M(x,t)^2 \rangle = \int_0^t dt_1 \int_0^t dt_2 e^{-2kt + k t_1 + k t_2} \langle \xi(x,t_1) \xi(x,t_2)\rangle$$ i.e. $$ \langle M(x,t)^2 \rangle = \int_0^t dt_1 \int_0^t dt_2 e^{-2kt + k t_1 + k t_2} 2 D \delta(x-x)\delta(t_1-t_2).$$ Yet here, in contrast to the case of when the Gaussian noise only depends on time, one encounters a divergence due to the Dirac delta of $0$: $$ \langle M(x,t)^2 \rangle = \frac{D}{k} \delta(0) \big( 1- e^{-2kt}\big) \rightarrow \infty.$$
How can such a simple model with symmetric fluctuations and damping generate infinite fluctuations? I found a very similar calculation in the book "Fractal concepts of surface growth" by Barbasi et al (1995). They follow the exact same steps as me without presenting any divergence (in eq. 4.12). They seem to simply set $\delta(0)=1$. What am I missing here?