Setup
Let $\Omega$ be a bounded domain with smooth boundary and $\zeta \in C^\infty_c(\Omega)$ be distributed according to a Gaussian measure such that $\zeta \sim \mathcal{N}(0,\mathcal{C}_\zeta)$ where $\mathcal{C}_\zeta$ is the covariance operator. Using the Karhunen-Loeve expansion, we can write $$ \zeta = \sum_k \lambda_k^{1/2}\eta_k\phi_k $$ where $\{\lambda_k,\phi_k\}$ are the eigenvalues/eigenvectors of $\mathcal{C}_\zeta$ (self-adjoint, trace-class operator) and $\eta_k \sim \mathcal{N}(0,1)$ is independent and identically distributed.
Now let $A_\zeta := I-\Delta + \zeta\nabla$ be a linear p.d. operator with domain $D_A = C^2_c(\Omega)$, parameterized by the random field $\zeta$. For example, if I am thinking right, the following holds for the expectation $$\mathbb{E}(A_\zeta) = I-\Delta. $$
Question
What I am trying to figure out, is a formula for the Gaussian measure that describes the range of $A_\zeta$, $\text{Ran}(A_\zeta)$. That is, if $$A_\zeta x = y \quad,\quad y \in C_c(\Omega)$$ defines a linear PDE (assume zero-Dirichlet boundary conditions), then $x = A_\zeta^{-1}y$ and there must be a Gaussian measure $\mu = \mathcal{N}(m,\mathcal{C})$, with mean $m$ and covariance operator $\mathcal{C}$ such that $$ x \sim \mu $$ Thus, in the end, we should be able to express $x$ via the K-L expansion as before and write $$ x = m + \sum_k \lambda_k^{1/2}\eta_k\phi_k $$ where $\{\lambda_k,\phi_k\}$ are now the eigenvalues/eigenvectors of $\mathcal{C}$.
Note
For a numerical approximation of the above, one could sample $\zeta$ using the K-L expansion, and compute the corresponding $x$ by solving the linear PDE in order to estimate the measure $\mu$. However, since this problem is linear I think I should be able to find a formula for the covariance operator $\mathcal{C}$.
Any ideas how I could formulate that?
Fixing $y \in C_c(\Omega)$, this is as far as I can get right now...
mean
$$ m = \mathbb{E}[x] = \mathbb{E}[A_\zeta^{-1}]y = (I-\Delta)^{-1}y \quad \text{since} \quad \mathbb{E}[\zeta] \equiv 0 \quad \text{(zero function)} $$
covariance op.
$$ \mathcal{C} = \mathbb{E}[(x-m)(x-m)^*] = \mathbb{E}[xx^*] - mm^* = \mathbb{E}[A_\zeta^{-1}xx^*(A_\zeta^{-1})^*] - mm^* $$
Does that seem correct?
I guess we would like to make the covariance operator $\mathcal{C}_\zeta := \mathbb{E}[\zeta\zeta^*]$ appear, but this inversion does not permit that right now.