Solution (source code)

= Solution

The <Fourier transform> turns $\nabla^2$ into multiplication by $-q^2$. Each Cartesian component is an <Ornstein-Uhlenbeck process>:
$$
\dot p_{q,i}=-r(q)p_{q,i}+f_{q,i},\qquad\boxed{r(q)=\Gamma(a+\kappa q^2).}
$$
Since $r(q)>0$, the transient from the remote past vanishes, and the integrating-factor solution is
$$
\boxed{p_q(t)=\int_{-\infty}^t f_q(s)e^{-r(q)(t-s)}\,ds.}
$$
Using the stated <Gaussian white noise> covariance and assuming $t_1\geq t_2$,
$$
\begin{aligned}
\langle p_{q,i}(t_1)p_{-q,j}(t_2)\rangle
&=2\Gamma k_BT\,\delta_{ij}\int_{-\infty}^{t_2}e^{-r(q)(t_1+t_2-2s)}\,ds\\
&=\delta_{ij}\frac{\Gamma k_BT}{r(q)}e^{-r(q)(t_1-t_2)}.
\end{aligned}
$$
Interchanging the times gives the stationary two-time <correlation function>
$$
\boxed{\langle p_{q,i}(t_1)p_{-q,j}(t_2)\rangle=\delta_{ij}\frac{k_BT}{a+\kappa q^2}e^{-\Gamma(a+\kappa q^2)|t_1-t_2|}.}
$$
At equal times this is the equilibrium <Gaussian field theory> covariance, and the decay time is $r(q)^{-1}$.