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')$.

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')$.

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

Given this form, it appears to be a multiplicative colored noise process. So writing down a Milstein scheme may be straightforward.

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

I've implemented the scheme, but there's something incorrect regarding the it$\hat{o}$ term. I'm wondering two things:

Is there a better way for me to derive the $\Delta t$ correction term in the Milstein?

Is this a case where switching to the Stratonovich form may be advantageous because I no longer have to deal with the correction term?