Essentially, what is asked is the continuation of my previous MO answer
Reformulation - Construction of thermodynamic limit for GFF
and the solution of the exercise I mentioned at the end of that answer.
There, I explained the construction of Gaussian Borel measures $\mu_m$ on the space $s'(\mathbb{Z}^d)$ of temperate multisequences indexed by the unit lattice in $d$ dimensions. The measure $\mu_m$ is specified by its characteristic function $$ p\longmapsto\exp\left(-\frac{1}{2}\sum_{x,y\in\mathbb{Z}^d}p(x)G_m(x,y)p(y)\right) $$ for $p=(p(x))_{x\in\mathbb{Z}^d}$ in $s(\mathbb{Z}^d)$, the space of multisequences with fast decay. The discrete Green's function $G_m(x,y)$ is defined on $\mathbb{Z}^d\times\mathbb{Z}^d$ by $$ G_m(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}(1-\cos \xi_j)}\ . $$ Here we will assume $m\ge 0$ for $d\ge 3$, and $m>0$ if $d$ is $1$ or $2$.
For any integer $N\ge 1$, define the discrete sampling map $\theta_N:\mathscr{S}(\mathbb{R}^d)\rightarrow s(\mathbb{Z}^d)$ which sends a Schwartz function $f$ to the multisequence $$ \left(f\left(\frac{x}{N}\right)\right)_{x\in\mathbb{Z}^d}\ . $$ This map is well defined and linear continuous. Indeed, $$ \langle Nx\rangle^2=1+\sum_{j=1}^{d} (Nx_j)^2\le N^2\langle x\rangle^2 $$ because $N\ge 1$. So $$ \sup_{x\in\mathbb{Z}^d} \langle x\rangle^k \left|f\left(\frac{x}{N}\right)\right| \le \sup_{z\in\mathbb{R}^d}\langle Nz\rangle^k|f(z)|\ \le N^k\ ||f||_{0,k} $$$$ ||\theta_N(f)||_k:= \sup_{x\in\mathbb{Z}^d} \langle x\rangle^k \left|f\left(\frac{x}{N}\right)\right| \le \sup_{z\in\mathbb{R}^d}\langle Nz\rangle^k|f(z)|\ \le N^k\ ||f||_{0,k} $$ where we used the standard seminorms $$ ||f||_{\alpha,k}=\sup_{z\in\mathbb{R}^d}\langle z\rangle^k|\partial^{\alpha}f(z)|\ . $$
Now$$ ||f||_{\alpha,k}=\sup_{z\in\mathbb{R}^d}\langle z\rangle^k|\partial^{\alpha}f(z)| $$ for Schwartz functions. Now consider the transpose map $\Theta_N=\theta_N^{\rm T}$ from $s'(\mathbb{Z}^d)$ to $\mathscr{S}'(\mathbb{R}^d)$. It is defined by $$ \langle \Theta_N(\psi),f\rangle=\langle\psi,\theta_N(f)\rangle=\sum_{x\in\mathbb{Z}^d}\psi(x)f\left(\frac{x}{N}\right) $$ for all discrete temperate fields $\psi$ and continuum test functions $f$. Essentially, $$ \Theta_N(\psi)=\sum_{x\in\mathbb{Z}^d}\psi(x)\ \delta_{\frac{x}{N}} $$ where $\delta_z$ denotes the $d$-dimensional Dirac Delta Function located at the point $z$. Now $\Theta_N$ is continuous for the strong topologies. Indeed if $A$ is a bounded subset of Schwartz space $$ ||\Theta_N(\psi)||_A=\sup_{f\in A}|\langle \Theta_N(\psi),f\rangle|= \sup_{p\in \theta_N(A)}|\langle \psi,p\rangle| $$ and $\theta_N(A)$ is bounded in $s(\mathbb{Z}^d)$ (because a continuous linear map sends bounded sets to bounded sets). Suppose we are given sequences $m_N$ and $\alpha_N$ dependent on the UV cutoff $N$. Define the Borel measure $$ \nu_N=(\alpha_N\Theta_N)_{\ast}\mu_{m_N} $$ on $\mathscr{S}'(\mathbb{R}^d)$. Its characteristic function is $$ W_N(f)=\int_{\mathscr{S}'(\mathbb{R}^d)}d\nu_N(\phi)\ e^{i\langle\phi,f\rangle} =\int_{s'(\mathbb{Z}^d)}d\mu_{m_N}(\psi)\ e^{i\langle\psi,\alpha_N\theta_N(f)\rangle} $$ by the abstract change of variable theorem. We then get $W_N(f)=\exp\left(-\frac{1}{2}Q_N(f)\right)$ where $$ Q_N(f)=\frac{\alpha_N^2}{(2\pi)^d}\sum_{x,y\in\mathbb{Z}^d} f\left(\frac{x}{N}\right)f\left(\frac{y}{N}\right) \int_{[0,2\pi]^d}d^d\xi\ \frac{e^{i\xi\cdot(x-y)}}{m^2+2\sum_{j=1}^{d}(1-\cos \xi_j)} $$ $$ =\frac{N^{2-d}\alpha_N^2}{(2\pi)^d}\sum_{x,y\in\mathbb{Z}^d} f\left(\frac{x}{N}\right)f\left(\frac{y}{N}\right) \int_{[-N\pi,N\pi]^d}d^d\zeta\ \frac{e^{i\zeta\cdot(\frac{x}{N}-\frac{y}{N})}}{N^2 m_N^2+2N^2\sum_{j=1}^{d}\left(1-\cos \left(\frac{\zeta_j}{N}\right)\right)} $$ after changing $[0,2\pi]^d$ to $[-\pi,\pi]^d$ by periodicity, then changing variables to $\zeta=N\xi$, and finally some algebraic rearrangement.
Pointwise in $\zeta\in\mathbb{R}^d$, we have $$ \lim\limits_{N\rightarrow\infty} 2N^2\sum_{j=1}^{d}\left(1-\cos \left(\frac{\zeta_j}{N}\right)\right) =\zeta^2 $$ and this is why I put an $N^2$ in the denominator. Finally, we can pick the right choice for the sequences $m_N$ and $\alpha_N$. For a fixed $m\ge 0$ (or strictly positive if $d=1,2$) we let $m_N=\frac{m}{N}$. Now we pick $\alpha_N$ so that the prefactor $N^{2-d}\alpha_N^2$ becomes the volume element $N^{-2d}$ for a Riemann sum approximation of a double integral on $\mathbb{R}^d\times\mathbb{R}^d$. Namely, we pick $\alpha_N=N^{-\frac{d}{2}-1}$. Equivalently, going back to $\alpha_N\Theta_N(\psi)$, that means choosing $$ \alpha_N\sum_{x\in\mathbb{Z}^d}\psi(x)\ \delta_{\frac{x}{N}}=\left(\frac{1}{N}\right)^{d-[\phi]} \sum_{x\in\mathbb{Z}^d}\psi(x)\ \delta_{\frac{x}{N}} $$ where $[\phi]=\frac{d-2}{2}$ is the (canonical) scaling dimensionscaling dimension of the free field. I wrote the last equation in a way to explicitly display the lattice spacing $\frac{1}{N}$.
Now an excellent exercise, for graduate students in analysis, is to show that $$ \lim\limits_{N\rightarrow \infty}Q_N(f)=\frac{1}{(2\pi)^d}\int_{\mathbb{R}^d} d^d\zeta\ \frac{|\widehat{f}(\zeta)|^2}{\zeta^2+m^2} $$ where the Fourier transform is normalized as $\widehat{f}(\zeta)=\int_{\mathbb{R}^d}d^dx\ e^{-i\zeta\cdot x} f(x)$. Finally, Fernique's version of the Lévy Continuity Theorem for $\mathscr{S}'(\mathbb{Z}^d)$, shows that the Borel measures $\nu_N$ converge weakly to the one obtained directly in the continuum using the Bochner-Minlos Theorem.