Reformulation - Construction of thermodynamic limit for GFF I've posted a question about the thermodynamic limit for Gaussian Free Fields (GFF) a couple days ago and I haven't got any answers yet but I kept thinking about it and I thought it would be better to reformulate my question and exclude the previous one, since now I can pose it in a more concrete way. The problem is basically give mathematical meaning to the infinite volume Gaussian measure associated to the Hamiltonian of the discrete GFF. In what follows, I will formulate the problem and then state the question.
A (lattice) field is a function $\phi: \Lambda_{L} \to \mathbb{R}$, where $\Lambda_{L} := \mathbb{Z}^{d}/L\mathbb{Z}^{d}$. Thus, the space of all fields is simply $\mathbb{R}^{\Lambda_{L}}$. The discrete Laplacian is the linear operator $\Delta_{L}:\mathbb{R}^{\Lambda_{L}}\to \mathbb{R}^{\Lambda_{L}}$ defined by:
\begin{eqnarray}
(\Delta_{L}\phi)(x) := \sum_{k=1}^{d}[-2\phi(x)+\phi(x+e_{k})+\phi(x-e_{k})] \tag{1}\label{1}
\end{eqnarray}
If $\langle \cdot, \cdot \rangle_{L}$ denotes the usual inner product on $\mathbb{R}^{\Lambda_{L}}$, we can prove that $\langle \phi, (-\Delta_{L}+m^{2})\phi\rangle_{L} > 0$ if $\langle \phi,\phi\rangle_{L}> 0$ and $m \neq 0$. Thus, $-\Delta_{L}+m^{2}$ defines a positive-definite linear operator on $\mathbb{R}^{\Lambda_{L}}$. We can extend these ideas to $\mathbb{R}^{\mathbb{Z}^{d}}$ as follows. A field $\phi: \mathbb{Z}^{d}\to \mathbb{R}$ is called $L$-periodic if $\phi(x+Ly) = \phi(x)$ for every $y \in \mathbb{Z}^{d}$. Let $\mathcal{F}_{per}$ be the set of all $L$-periodic fields, so that $\mathcal{F}_{per} \subset \mathbb{R}^{\mathbb{Z}^{d}}$. Now, using the same expression in (\ref{1}) we can define 'infinite volume' Laplacians $\Delta_{per}$ and $\Delta$ on $\mathcal{F}_{per}$ and $l^{2}(\mathbb{Z}^{d}):=\{\phi:\mathbb{R}^{d}\to \mathbb{R}:\hspace{0.1cm} \sum_{x \in \mathbb{Z}^{d}}|\psi(x)|^{2}<\infty\}$, respectivelly. In addition, if $\phi \in \mathcal{F}_{per}$, its restriction $\phi|_{\Lambda_{L}}$ can be viewed as an element of $\Lambda_{L}$, and the action of $\Lambda_{per}$ to $\phi|_{\Lambda_{L}}$ is equivalent to the action of $\Delta_{L}$ to $\phi|_{\Lambda_{L}}$. 
The Hamiltonian for the GFF in the lattice $\Lambda_{L}$ is given by:
\begin{eqnarray}
H_{\Lambda_{L}}(\phi) = \frac{1}{2}\langle \phi, (-\Delta_{L}+m^{2})\phi\rangle_{\Lambda} \tag{2}\label{2}
\end{eqnarray} 
The first step is to define finite volume measures on $\mathbb{R}^{\mathbb{Z}^{d}}$. For each finite $\Lambda \subset \mathbb{Z}^{d}$, let $C_{\Lambda} =(C_{xy})_{x,y \in \Lambda}$ be the matrix with entries $C_{xy} := (-\Delta_{per}+m^{2})_{xy}$, where $(-\Delta_{per}+m^{2})_{xy}$ is the Kernel of $-\Delta_{per}+m^{2}$ on $\mathbb{R}^{\mathbb{Z}^{d}}$. Because the Kernel of $-\Delta_{per}+m^{2}$ is the same as $-\Delta_{L}+m^{2}$, each $C_{\Lambda}$ is a positive-definite matrix and, thus, define a Gaussian measure $\mu_{\Lambda}$ on $\mathbb{R}^{\Lambda}$. Because this family of Gaussian measures $\mu_{\Lambda}$ is consistent, we can use Kolmogorov's Extension Theorem to obtain a Gaussian measure $\mu$ on $\mathbb{R}^{\mathbb{Z}^{d}}$ (with the product $\sigma$-algebra). Moreover, we can also obtain a family $\{f_{\alpha}\}_{\alpha \in \mathbb{Z}^{d}}$ of random variables such that $\mu_{\Lambda}$ is the joint probability distribution of $\{f_{\alpha}\}_{\alpha \in \Lambda}$. As it turns out, it is possible to prove that these random variables are given by $f_{\alpha}(\phi) = \phi(\alpha)$, $\alpha \in \mathbb{Z}^{d}$. In summary, if $A$ is a Borel set in $\mathbb{R}^{\Lambda}$, we must have:
\begin{eqnarray}
\mu_{\Lambda_{L}}(A) = \frac{1}{Z}\int_{A}e^{-\frac{1}{2}\langle \phi, (-\Delta_{L}+m^{2})\phi\rangle_{L}}d\nu_{L}(\phi) = \mu(A\times \mathbb{R}^{\mathbb{Z}^{d}\setminus \Lambda}) \tag{3}\label{3}
\end{eqnarray} 
with $\nu_{L}$ being the Lebesgue measure on $\mathbb{R}^{\Lambda}$. The Gaussian measure $\mu$ is our a priori measure on $\mathbb{R}^{\mathbb{Z}^{d}}$ and, by (\ref{3}), it can be interpreted as a finite volume over $\mathbb{R}^{\Lambda}$. 
Now, let $G(x,y)$ the Green function of $-\Delta+m^{2}$ in $l^{2}(\mathbb{Z}^{d})$. If $s_{m}:=\{\psi \in \mathbb{R}^{\mathbb{N}}:\hspace{0.1cm} \sum_{n=1}^{\infty}n^{2m}|\psi_{n}|^{2}\equiv ||\psi||_{m}^{2}\infty\}$, define $s:=\bigcup_{m\in \mathbb{Z}}$ and $s':=\bigcap_{m\in \mathbb{Z}}s_{m}$. Let u $K=(K_{xy})_{x,y \in \mathbb{Z}^{d}}$ be an 'infinite matrix' given by $K_{xy}:= G(x,y)$. If we order $\mathbb{Z}^{d}$, we can consider $K$ to be an 'infinite matrix' $K = (K_{ij})_{i,j \in \mathbb{N}}$. Now, define the following map:
\begin{eqnarray}
s \times s \ni (\psi,\varphi) \mapsto (\psi,K\varphi):= \sum_{i,j=1}^{\infty}\psi_{i}K_{ij}\varphi_{j}\tag{4}\label{4}
\end{eqnarray}
Let $W(\phi):= e^{-\frac{1}{2}(\phi,K\phi)}$. It is possible to prove that $W$ is a positive-definite function on $s$, so that, by Minlos Theorem, there exists a Gaussian measure $\tilde{\mu}_{K}$ on $s'$ such that $W$ is the Fourier transform of $\tilde{\mu}_{K}$. 
[Question] I would like to establish a connection between $\mu$ and $\tilde{\mu}_{K}$ (where, here, $\mu$ is the restriction of $\mu$ to $s'\subset \mathbb{R}^{\mathbb{Z}^{d}}$ with its natural $\sigma$-algebra). It seems to me that $\tilde{\mu}_{K}$ is the infinite volume measure of $\mu$, in the sense that when we take $L\to \infty$ one should obtain $\tilde{\mu}_{K}$. In other words, $\tilde{\mu}_{K}$ is the infinite volume Gibbs measure obtained by taking the thermodynamic limit of the measures $\mu_{L}$. But, if I'm not mistaken, to prove that $\tilde{\mu}_{K}$ is the infinite volume Gibbs measure, I should prove that:
\begin{eqnarray}
\lim_{L\to \infty}\int f(\phi)d\mu_{L}(\phi) = \int f(\phi)d\tilde{\mu}_{K}(\phi) \tag{5}\label{5}
\end{eqnarray}
i.e. I should prove that $\mu$ converges weakly to $\tilde{\mu}_{K}$. And I don't know how to prove it.
Note: The above setup is a result of some of my thoughts about the problem during the last few days. I've been using a lot of different references and each one work the problem in a different way, with different notations and different objectives, so I'm trying to put it all together in one big picture. It is possible that my conclusions are not entirely correct or I can be going in the wrong direction, idk. But any help would be appreciated. 
Note 2: I think it could be easier to prove a more particular limit such as $\lim_{L\to \infty}\int\phi(x)\phi(y)d\mu_{L}(\phi) = \int \phi(x)\phi(y)d\tilde{\mu}_{K}(\phi)$ and this would be enough to establish the existence of infinite volume correlation functions, which is one of the most important quantities in statistical mechanics. However, I don't think I could conclude that $\tilde{\mu}_{K}$ is the associate infinite volume Gibbs measure for the system just from this limit. Don't I need to prove it for a general $f$ as above?
 A: For $x\in\mathbb{Z}^d$ I will denote by $\bar{x}$ the corresponding equivalence class in the discrete finite torus $\Lambda_{L}=\mathbb{Z}^d/L\mathbb{Z}^d$.
I will view a field $\phi\in\mathbb{R}^{\Lambda_L}$ as a column vector with components $\phi(\bar{x})$ indexed by $\bar{x}\in\Lambda_L$.
The discrete Laplacian $\Delta_L$ then acts by
$$
(\Delta_L\phi)(\bar{x})=\sum_{j=1}^{d}\left[
-2\phi(\bar{x})+\phi(\overline{x+e_j})+\phi(\overline{x-e_j})
\right]\ .
$$
Now take the column vectors
$$
u_k(\bar{x})=\frac{1}{L^{\frac{d}{2}}} e^{\frac{2i\pi k\cdot x}{L}}
$$
for $k=(k_1,\ldots,k_d)\in\{0,1,\ldots,L-1\}^d$.
They give an orthonormal basis in $\mathbb{C}^{\Lambda_L}$ which diagonalizes the Laplacian matrix.
Let $C_L=(-\Delta_L+m^2{\rm I})^{-1}$ and denote its matrix elements by $C_L(\bar{x},\bar{y})$.
We then have, for all $x,y\in\mathbb{Z}^d$,
$$
C_L(\bar{x},\bar{y})=
\frac{1}{L^d}\sum_{k\in\{0,1,\ldots,L-1\}^d}
\frac{e^{\frac{2i\pi k\cdot(x-y)}{L}}}{m^2+
2\sum_{j=1}^{d}\left[1-\cos\left(\frac{2\pi k_j}{L}\right)\right]}
=:G_L(x,y)
$$
where we used the formula to define $G_L$ on $\mathbb{Z}^d\times\mathbb{Z}^d$.
Because we assumed $m>0$, we have the trivial uniform bound in $L$ saying
$$
|G_L(x,y)|\le \frac{1}{m^2}\ .
$$
Now let $\nu_L$ denote the centered Gaussian probability measure on $\mathbb{R}^{\Lambda_L}$ with covariance matrix $C_L$. We also define an injective continuous linear map
$$
\tau_L:\mathbb{R}^{\Lambda_L}\longrightarrow s'(\mathbb{Z}^d)
$$
which sends $\phi\in\mathbb{R}^{\Lambda_L}$ to $\psi\in s'(\mathbb{Z}^d)$ defined by
$\psi(x)=\phi(\bar{x})$ for all $x\in\mathbb{Z}^d$.
Of course $\mathbb{R}^{\Lambda_L}$ has its usual finite-dimensional space topology, whereas
$s'(\mathbb{Z}^d)$ is given the strong topology and the resulting Borel $\sigma$-algebra.
As I explained in my answer to the previous MO question we can use such a map to push forward probability measures. We thus go ahead and define $\mu_L=(\tau_L)_{\ast}\nu_L$ which is a Borel probability measure on $s'(\mathbb{Z}^d)$.
Now we switch gears and consider the Green's function $G_{\infty}(x,y)$ for $-\Delta+m^2$ on $\mathbb{Z}^d$. More explicitly,
$$
G_{\infty}(x,y)=\frac{1}{(2\pi)^d}
\int_{[0,2\pi]^d}d^d\xi\ 
\frac{e^{i\xi\cdot(x-y)}}{m^2+
2\sum_{j=1}^{d}\left(1-\cos\xi_j\right)}\ .
$$
The function
$$
\begin{array}{crcl}
W_{\infty}: & s(\mathbb{Z}^d) & \longrightarrow & \mathbb{C} \\
 & f & \longmapsto & \exp\left(
 -\frac{1}{2}\sum_{x,y\in\mathbb{Z}^d} f(x)\ G_{\infty}(x,y)\ f(y)
 \right)
\end{array}
$$
satisfies all the hypotheses of the Bochner-Minlos Theorem for $s'(\mathbb{Z}^d)$. Therefore it is the characteristic function of a Gaussian Borel probability measure $\mu_{\infty}$ on $s'(\mathbb{Z}^d)$.
Finally after all these preliminaries we can state the main result which the OP asked for.
Theorem: When $L\rightarrow\infty$, the measure $\mu_L$ converges weakly to $\mu_{\infty}$.
The proof uses the Lévy Continuity Theorem for $s'(\mathbb{Z}^d)$ which is due to Xavier Fernique. One only has to prove that for all discrete test function $f\in s(\mathbb{Z}^d)$,
$$
\lim\limits_{L\rightarrow \infty} W_L(f)\ =\ W_{\infty}(f)
$$
where $W_L$ is the characteristic function of the measure $\mu_L$.
By definition, we have
$$
W_L(f)=\int_{s'(\mathbb{Z}^d)} \exp\left[
i\sum_{x\in\mathbb{Z}^d}f(x)\psi(x)
\right]\ d[(\tau_L)_{\ast}\nu_L](\psi)\ .
$$
By the abstract change of variable theorem,
$$
W_L(f)=\int_{\mathbb{R}^{\Lambda_L}} \exp\left[
i\sum_{x\in\mathbb{Z}^d}f(x)\phi(\bar{x})
\right]\ d\nu_L(\phi)
$$
$$
=\exp\left(
-\frac{1}{2}\sum_{\bar{x},\bar{y}\in\Lambda_L}
\tilde{f}(\bar{x})\ C_L(\bar{x},\bar{y})\ \tilde{f}(\bar{y})
\right)
$$
where we introduced the notation $\tilde{f}(\bar{x})=\sum_{z\in\mathbb{Z}^d}f(x+Lz)$.
Hence
$$
W_L(f)= \exp\left(
 -\frac{1}{2}\sum_{x,y\in\mathbb{Z}^d} f(x)\ G_{L}(x,y)\ f(y)
 \right)\ .
$$
Since the function
$$
\xi\longmapsto \frac{1}{(2\pi)^d}
\frac{e^{i\xi\cdot(x-y)}}{m^2+
2\sum_{j=1}^{d}\left(1-\cos\xi_j\right)}
$$
is continuous on the compact $[0,2\pi]^d$ and therefore uniformly continuous, we have that, for all fixed $x,y\in\mathbb{Z}^d$, the Riemann sums $G_L(x,y)$ converge to the integral $G_{\infty}(x,y)$.
Because of our previous uniform bound on $G_L(x,y)$ and the fast decay of $f$, we can apply the discrete Dominated Convergence Theorem in order to deduce
$$
\lim\limits_{L\rightarrow\infty}
\sum_{x,y\in\mathbb{Z}^d} f(x)\ G_{L}(x,y)\ f(y)\ =\ 
\sum_{x,y\in\mathbb{Z}^d} f(x)\ G_{\infty}(x,y)\ f(y)\ .
$$
As a result $\lim_{L\rightarrow \infty}W_L(f)=W_{\infty}(f)$ and we are done.
Note that we proved weak convergence which as usual means
$$
\lim\limits_{L\rightarrow \infty}
\int_{s'(\mathbb{Z}^d)}F(\psi)\ d\mu_L(\psi)\ =\ 
\int_{s'(\mathbb{Z}^d)}F(\psi)\ d\mu_{\infty}(\psi)
$$
for all bounded continuous functions $F$ on $s'(\mathbb{Z}^d)$.
One also has convergence for correlation functions or moments because of the Isserlis-Wick Theorem relating higher moments to the second moment and the previous argument where we explicitly treated the convergence for second moments.
Finally, note that the extension map $\tau_L$ used here is the periodization map, but there are lots of other choices which work equally well. A good exercise, is to construct the massive free field in the continuum, i.e., in $\mathscr{S}'(\mathbb{R}^d)$, as the weak limit of suitably rescaled lattice fields on $\mathbb{Z}^d$ with a mass adjusted as a function of the (rescaled) lattice spacing.
