I am broadly working on a dynamic process where we want to see how a field $\rho(r)$ changes in space in time with thermal noise. The system is biased around a thermodynamic saddle point dictated by $H[\rho]$, and $\nabla^2$ is the typical Laplacian operator that introduces diffusion.

$\frac{\partial \rho(r)}{\partial t} = \nabla \cdot \left[ \rho(r)\nabla \left( \frac{\partial H[\rho(r)]}{\partial \rho}  \right)\right]+\eta(t,r)$

In this case, $\eta(t,r)$ is not white Gaussian noise.

$\left<\eta(t,r)\right> = 0$

$\left<\eta(t,r) \eta(t',r')\right> = -2 \delta(t,t')\nabla\cdot \left[ \rho(r)\nabla \delta(r,r')\right]$

I've left out factors of $\Delta V$ related to discretization. Typically, we adopt an it$\hat{\text{o}}$ interpretation when we numerically integrate the system given above. We typically end up with the Euler–Maruyama scheme:

$\rho_{n+1}(r) = \rho_{n}(r)+\Delta t \nabla^2 \left( \frac{\partial H[\rho(r)]}{\partial \rho} \right)+R_n(r)$ 

I've only discretized in time. Where $R_n$ is just a random number computed using the construction above. The main struggle with this stochastic partial differential equation is that $\rho$ describes a physical quantity ie. density, and cannot be negative. From my understanding, it is a strong order problem. Thus I would like to make an approximation to the noise:

$\left<\eta(t,r) \eta(t',r')\right> = -2 \rho(r)\delta(t,t')\nabla^2 \delta(r,r')$.