Let $(\Omega,\mathcal F, P)$ be a probability space, let $\mathcal G_1$ and $\mathcal G_2$ be two (not necessarily independent) centered real-valued Gaussian processes, i.e., measurable functions $\Omega\rightarrow \mathbb R^{[0,1]}$ (here $\mathbb R^{[0,1]}$ denotes the space of all functions $[0,1]\rightarrow\mathbb R$), and let $F_1$ and $F_2$ be two continuous functions $[0,1]\rightarrow\mathbb R$ of bounded variation. Put $$I_1(\omega) := \int_0^1\mathcal G_1(\omega)(t)\,\mathrm dF_1(t)\qquad\text{and}\qquad I_2(\omega) := \int_0^1\mathcal G_2(\omega)(t)\,\mathrm dF_2(t),$$i.e., $I_1$ and $I_2$ are the Stieltjes integrals of the Gaussian processes $\mathcal G_1$ and $\mathcal G_2$ with respect to the functions $F_1$ and $F_2$, respectively. Note that both $(\mathcal G_1(\cdot)(t))_{t\in[0,1]}$ and $(\mathcal G_2(\cdot)(t))_{t\in[0,1]}$ are a collections of random variables. My goal is now to compute $$\mathsf E[I_1 I_2] := \int I_1(\omega) I_2(\omega)\,\mathrm dP.$$
My attempt so far: By definition of the Stieltjes integral, $$I_k(\omega) = \lim_{n\rightarrow\infty}\sum_{[a_k,b_k]\in\Pi_{k,n}} \mathcal G_k(\omega)(c_k)\big(F_k(b_k) - F_k(a_k)\big),$$ where $(\Pi_{k,n})_{n\in\mathbb N}$ is a sequence of partitions of $[0,1]$ and $c_k$ is always between $a_k$ and $b_k$. Hence \begin{align*}\mathsf E[I_1 I_2] &= \mathsf E\left[\left(\lim_{n_1\rightarrow\infty}\sum_{[a_1,b_1]\in\Pi_{1,n_1}} \mathcal G_1(\omega)(c_1)\big(F_1(b_1) - F_1(a_1)\big)\right)\left(\lim_{n_2\rightarrow\infty}\sum_{[a_2,b_2]\in\Pi_{2,n_2}} \mathcal G_2(\omega)(c_2)\big(F_2(b_2) - F_2(a_2)\big)\right)\right] \\ &= \mathsf E\left[\lim_{n_2\rightarrow\infty}\lim_{n_1\rightarrow\infty}\left(\sum_{[a_1,b_1]\in\Pi_{1,n_1}} \mathcal G_1(\omega)(c_1)\big(F_1(b_1) - F_1(a_1)\big)\right)\left(\sum_{[a_2,b_2]\in\Pi_{2,n_2}} \mathcal G_2(\omega)(c_2)\big(F_2(b_2) - F_2(a_2)\big)\right)\right] \\ &= \mathsf E\left[\lim_{n_2\rightarrow\infty}\lim_{n_1\rightarrow\infty}\sum_{[a_1,b_1]\in\Pi_{1,n_1}}\sum_{[a_2,b_2]\in\Pi_{2,n_2}} \mathcal G_2(\omega)(c_2) \mathcal G_1(\omega)(c_1)\big(F_1(b_1) - F_1(a_1)\big)\big(F_2(b_2) - F_2(a_2)\big)\right]. \end{align*} Here I used the definition for the first equality, the assumption that the integrals exist for the second equation, and then simple algebra for the last equality. These steps should be fine, I guess. In the next step I would like to move out the limits, but I am not sure whether this is possible. If it was possible, I would have to compute $\mathsf E[\mathcal G_1(\cdot)(c_1)\mathcal G_2(\cdot)(c_2)]$. If both processes are independent, this expectation is zero. Otherwise, the expectation would be equal to the corresponding covariance value. Is my reasoning correct? What can I use to move out the limits out of the expectation?
PS: I would also appretiate suggestions to improve the notation for the Stieltjes integral sum.
The expectation you want to calculate is, by measure theoretic Fubini, $$\int_0^1\int_0^1\mathbb E[{\cal G}_1(t)\,{\cal G}_2(s)]\,dF_1(t)\,dF_2(s)=\int_0^1\int_0^1\operatorname{Cov}[{\cal G}_1(t)\,{\cal G}_2(s)]\,dF_1(t)\,dF_2(s)$$ because the Gaussian processes are centered.