My answer Wiener measure and Bochner Minlos can help defining a free massless real scalar field as a random element of $S_0'(\mathbb{R}^2)$. Namely, take the Schwartz space of rapidly decaying test functions $S(\mathbb{R}^2)$ and consider the subspace $S_0(\mathbb{R}^2)=\{f| \widehat{f}(0)=0\}$ of "charge-neutral" test functions. The bilinear form $$ B(f,g)=\int_{\mathbb{R}^2} \frac{d^2\xi}{(2\pi)^2}\ \frac{\overline{\widehat{f}(\xi)}\widehat{g}(\xi)}{|\xi|^2} $$ is continuous and positive on $S_0(\mathbb{R}^2)$ and therefore by the Bochner-Minlos Theorem there exists a unique centered Gaussian probability measure $\mu$ on the topological dual $S_0'(\mathbb{R}^2)$ for which $$ \mathbb{E}(\phi(f)\phi(g))=B(f,g) $$ where $\phi$ is the corresponding random element in $S_0'(\mathbb{R}^2)$.
I did not do the computation (which needs a lot of care) but I suspect that $B(f,g)$ should be a multiple of $$ \int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)(-\log|x-y|)g(y)\ . $$ The issue here is that you have both a UV and an IR problem to deal with. The slow decay of the propagator for large $\xi$ makes it so that $X$ or $\phi$ is a random distribution rather than a random function and, in particular, punctual evaluations $X(x)$ do not make sense. Also, the divergence at $\xi=0$, or zero-mode, makes it so that there is an ambiguity of shifting the field by a constant so the field $X$ does not make sense but its "increments" do. Working with $\partial X$ is another way to circumvent this issue.
Edit: I just did the computation and indeed $$ B(f,g)=\frac{1}{2\pi}\int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)(-\log|x-y|)g(y)\ . $$ This can be done by replacing $1/|\xi|^2$ by $1/|\xi|^{2-\epsilon}$ and taking the $\epsilon\rightarrow 0$ limit by dominated convergence. Then, use the representation $$ \frac{1}{|\xi|^{2-\epsilon}}=\frac{1}{\Gamma\left(\frac{2-\epsilon}{2}\right)} \int_{0}^{\infty} d\alpha\ \alpha^{-\frac{\epsilon}{2}} e^{-\alpha |\xi|^2}\ , $$ write the Fourier transforms as integrals and integrate over $\xi$. One then ends up with the computation of $$ \lim_{\epsilon\rightarrow 0^+} \int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)\left(\frac{1}{\epsilon} |x-y|^{-\frac{\epsilon}{2}}\right)g(y)\ . $$ Since $\int f(x) d^2 x=\int g(x) d^2 x=0$, the last expression is the same as $$ \lim_{\epsilon\rightarrow 0^+} \int_{\mathbb{R}^2\times\mathbb{R}^2} d^2 x\ d^2y\ f(x)\left\{\frac{1}{\epsilon} \left(|x-y|^{-\frac{\epsilon}{2}}-1\right)\right\}g(y)\ . $$ The resulting derivative at $\epsilon =0$ produces the wanted logarithm of $|x-y|$.
Edit 2: A substantial elaboration on the answer I gave above has appeared recently. See the review "Log-correlated Gaussian fields: an overview" by Duplantier, Rhodes, Sheffield and Vargas. This is part of a wider program regarding the study of stationary, isotropic and self-similar Gaussian fields:
- "Fractional Gaussian fields: a survey" by Lodhia, Sheffield, Sun and Watson.
- Chapter 2 of "Construction and Analysis of a Hierarchical Massless Quantum Field Theory", Ph.D. Thesis by Ajay Chandra.
- "Gaussian and their Subordinated Self-similar Random Generalized Fields" by Roland Dobrushin and its follow-up " Multiple Wiener-Itô Integrals" by Péter Major.
