we consider driving by stochastic sources and correspondingly density matrices obtained from, for example, frequency, spatial or time averages (when considering time dependent, in particular, non-harmonic, stochastic sources) [@EMC15]. Note that we always assume here that $\hat{T}$ itself is not explicitly time dependent, even if the driving is.

Consider a source correlation function $\Gamma_0$ with $$\label{Gx1x2}
\Gamma_0(x_1,x_2) = \braket{\psi_0^-(x_1) \psi^-_0(x_2)^*},$$ where $\psi_0^-(x)$ is the source wave field as defined in (\[psim\]) and $\braket{.}$ denotes an ensemble average over time intervals, frequency or local spatial averaging. Note, that while the source distribution $\ket{\psi_0^-}$ show strong spatial fluctuations, the averaged quantity ${\Gamma_0}$ can be a fairly smooth function of both $x_1$ and $x_2$. Starting from (\[prop\]), we write the system correlation function, including multiple reflections at boundaries, in terms of the source correlation function and the transfer operator $\hat{T}$. That is, $$\label{corr-def}
\hat{\Gamma}  = (1 - \hat{T})^{-1} \hat{\Gamma}_0 (1 - \hat{T}^{\dagger})^{-1}.$$ Formally, one can write this in terms of the source correlation function $\Gamma_0$ being propagated along outgoing waves and undergoing multiple reflections at boundaries. After geometrically expanding the right hand side of Eq. (\[corr-def\]), we write $$\label{G0}
\hat{\Gamma}  = \sum_{n,n' = 0}^{\infty} \hat{T}^{n} \hat{\Gamma}_0 (\hat{T}^{\dagger})^{n'} = \hat{K} + 
\sum_{m=1}^{\infty} \left(\hat{T}^m \hat{K} 