Free Boson Correlator $ \langle X(z)X(w) \rangle =- \ln |z - w| $ In physics papers, the massless free boson has a definition involving an action:
$$ S(X) = \frac{1}{8\pi} \int d\sigma^2\,   \partial X \overline{\partial X}$$
The random functions $X(z)$ are sampled according to the Gaussian distribution: $e^{-S(X)}$.  In this case we can define the correlator to be the vacuum expectation value:
$$ \mathbb{E}[X(z,\overline{z})X(w, \overline{w})] \equiv 
\langle \varnothing | X(z,\overline{z})X(w, \overline{w}) |\varnothing \rangle =
- \tfrac{1}{2}\ln |z - w|$$
See Ginsparg, applied Conformal Field Theory. I am trying to fill in the details of this derivation... using probability theory, complex analysis or anything else.

Since $X$ is a bosonic field I will just model it as a function... judging from the norm it should live in some Hardy space $X \in H^2$ - the complexification of $L^2$ - where we should fill in some Riemann surface with flat structure. I think Ginsparg uses $H^2(\mathbb{R}\times S^1)$.
The equations of motion ought to be the Cauchy Riemann equations.  The fields in this theory are holomorphic functions. $X(z) = \sum a_n z^n$.  Then I would like to expand.  $\partial X(z) = \sum n a_n z^n$ and get
$$ \langle \partial X(z)\partial X(w) \rangle = \sum n\mathbb{E}[a_n^2]\,(z \overline{w})^n  = \frac{1}{|1-z\overline{w}|^2}$$
There are holes everywhere in this argument - a real physicists's derivation.  How can I make this more rigorous?  Especially the intuition about the Gaussian random function in $L^2$ ? 
Thanks.
 A: In massless 2d scalar field theory, the expression $\langle \phi(x)\phi(y) \rangle = -\ln|x-y|$ is best regarded as a mnemonic.  (A first hint of trouble:  $|x-y|$ has dimensions of length, so you can't just stick it into a logarithm.)  Likewise, the expression $e^{-S}$ for the measure.  The Euclidean measure is usually supported on functions for which $S$ isn't defined.
The idea behind this shorthand is that any Euclidean correlation function in the massless 2d scalar model can be written as a limit of a correlation function in the massive model.  The massive 2d theory is Gaussian and completely characterized by the 2-point function $\langle \phi(x)\phi(y)\rangle$, which blows up like $-\ln(m|x-y|)$ as $m \to 0$.  This means that some correlation functions will not remain finite in this limit, which reduces the size of the Hilbert space and the collection of field operators.  
Concretely, this works as follows.. 
Recall the Wightman/GNS construction of the Hilbert space and field operators:  The collection of Euclidean correlation functions specifies a linear functional $E$ on the algebra $\mathcal{A}_+$ of Euclidean observables which are supported at positive time.  Reflection about $t=0$ gives a map $\theta: \mathcal{A}_+ \to \mathcal{A}_-$, which defines a (degenerate) sesquilinear form $s: \mathcal{A}_+ \times \mathcal{A}_+ \to [0,\infty]$.  Call an element $a$ of $\mathcal{A}_+$ $L^2$-finite if $E(\theta(a)a) < \infty$, by analogy with the case where $E$ is the expectation value of a normed functional measure.  If $s$ is reflection positive, we get an inner product space by restricting attention to the space $\mathcal{V}$ of $L^2$-finite vectors and quotienting out by the null ideal of $s$.   Completing gives a Hilbert space $\mathcal{H}$.  The elements of $\mathcal{V}$ act on a space of smooth vectors within $\mathcal{H}$ via right multiplication.
When we take the massless limit, we reduce the size of $\mathcal{V}$, as the evaluation observables $ev_x(\phi) = \phi(x)$ cease to $L^2$-finite.  (Note that $E(\theta(ev_x) ev_y) = \langle \phi(x)\phi(y)$.)  This is why people say $\phi$ isn't a quantum field in the massless case.  However, many local observables do survive this limit, such as the famous vertex operators.  (And various idiocies occur, such as zero becoming part of the Hamiltonian's continuous spectrum, thus 'no vacuum vector'.)
It's harmless abuse of notation to write expressions in terms of $\phi$ (or its Fourier components, which you've written $a_n$), as long as you remember that the computation has to commute with the construction above.
A: 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, pointwise 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.

