Solution (source code)

= Solution

Choose an orthonormal basis $(e_n)$ of $H_0^1(D)$ and write the GFF formally as $h=\sum_n\xi_ne_n$, where the $\xi_n$ are independent standard normal variables. The <Green-kernel expansion in the Dirichlet space> gives
$$
\sum_{n\geq1}\left(\int_De_n(x)\rho(dx)\right)^2
=\iint_{D\times D}G_D(x,y)\rho(dx)\rho(dy)<\infty.
$$
Consequently the series
$$
(h,\rho):=\sum_{n\geq1}\xi_n\int_De_n(x)\rho(dx)
$$
converges in $L^2$. Its partial sums are centered Gaussian and their variances converge to the displayed <Green energy>. The $L^2$ limit is therefore Gaussian with mean zero and variance
$$
\iint_{D\times D}G_D(x,y)\rho(dx)\rho(dy).
$$
This constructs the <Finite-Green-energy measure pairing with a Gaussian free field> independently of the chosen orthonormal basis.