I mentioned "Wiener" process in quotes because the Wiener term that appears in most SDEs are functions of time. In fact, my SDE is a function of space ($z$). Suppose I have the following first-order stochastic vector differential equation (in space): $$ \frac{d\textbf{u}}{dz} = \left[\textbf{G}+\textbf{H}\xi(z)\right]\textbf{u} $$ where $\textbf{u}$ is a $n$-dimensional column vector, $\textbf{G}$ and $\textbf{H}$ are constant $n\times n$ matrices, and $\xi(z)$ is a (possibly colored) Gaussian noise process with zero mean and autocorrelation function $$ R_{\xi}(\zeta)=\left<\xi(z)\xi(z+\zeta)\right>. $$ If I define an equivalent Wiener process such that $$ W(z)=\int_{0}^{z}\xi(z^{\prime})dz^{\prime}, $$ how do I show that the mean of $W(z)$ squared is derived as $$ \left<W^{2}\right> = \left<\int_{0}^{z}\xi(z_{1})dz_{1}\int_{0}^{z}\xi(z_{2})dz_{2}\right> = 2\int_{0}^{z}(z-\zeta)R_{\xi}(\zeta)d\zeta $$ ?
For reference, the paper that I'm following is here: https://opg.optica.org/josab/fulltext.cfm?uri=josab-15-8-2269&id=35539 and the equations are (A1)-(A5) in the appendix