This question was previously asked here
https://physics.stackexchange.com/questions/251927/two-point-function-of-a-free-massless-scalar-field-in-euclidean-space-time
though I did not get there an answer.
Let $\phi(x)$ be a quantum free massless scalar field (i.e. operator valued function) on a Euclidean space-time of dimension $D>2$. Thus it satisfies the equation
$$\Delta\phi(x)=0,$$
where $\Delta =\sum_{i=1}^D(\partial_i)^2$ is the Laplacian.
I would like to compute the function $$D(x):=\langle 0|\phi(x)\phi(0)|0\rangle$$
(there is no time ordering!).
The argument below shows that $D(x)$ is a constant function; that does not sound to make sense. Indeed one has
$$\Delta D(x)=\langle 0|\left(\Delta \phi(x)\right)\phi(0)|0\rangle=0.$$
Since any harmonic generalized function is smooth (elliptic regularity), $D(x)$ is a smooth function. Moreover $D(x)$ is $SO(D)$-invariant. Hence on $\mathbb{R}^D\backslash\{0\}$ as a function of $r=|x|$, $D(r)$ satisfies some linear second order differential equation (which can easily be written down). The space of solutions is 2-dimensional. Hence $D(x)$ is a linear combination of the constant function and of $\frac{1}{|x|^{D-2}}$. Since $D(x)$ is smooth, it must be constant.
Remark. To compare with, if instead of Euclidean metric one considers Minkowski metric on $\mathbb{R}^{D}=\mathbb{R}^{1+d}$, then the above argument shows $\Box D(x)=0$ where $\Box =\partial_0^2-\sum_{i=1}^d\partial_i^2$, and $D(x)$ is Lorentz-invariant. In fact the correct answer is $D(x)=\int \frac{d^dp}{(2\pi)^d}\frac{1}{2|p|}e^{-i(|p|x^0-\sum_{i=1}^dp_ix^i)}=\int \frac{d^{d+1}p}{(2\pi)^d}\delta(p^2)\theta(p_0)e^{-i\sum_{\mu=0}^dp_\mu x^\mu}$; it is obtained by an explicit construction of the quantum field $\phi(x)$. Here $\delta(p^2)\theta(p_0)$ can be considered as the only (up to proportionality) Lorentz invariant measure on the light cone; it is finite in a neighborhood of 0 for $D>2$ hence its Fourier transform is well defined in sense of distributions.
Best Answer
Your argument starts with the assumption that there is such a thing as a "Euclidean free scalar field", as an operator valued function or distribution. And this is where it goes wrong.
Obviously, you are thinking of a Minkowski space free scalar quantum field, an operator valued distribution $\phi(x)$, satisfying $\square \phi(x) = 0$, which you would then like to "analytically continue" to complex values of the coordinates $x$ (Wick rotation). After the Wick rotation, you expect to still have an operator valued distribution, but satisfying $\Delta \phi(x) = 0$ instead, which you could then use to take products, expectation values etc. Or rather, as your initial hypothesis indicates, you would like to skip all these preliminary steps and reason directly with the Euclidean operator valued $\phi(x)$.
As I'm sure you realize, analytic continuation is not applicable to arbitrary distributions. So, expecting the distribution $\phi(x)$ to have a well-defined analytic continuation (Wick rotation) in the $x$ coordinates is wishful thinking.
What is actually done in practice is the following. First, one computes the Minkowski time-ordered 2-point function $G(x) = \langle 0 | T \phi(x) \phi(0) | 0 \rangle$. Then, either by direct computation in specific constructions or axiomatically, one affirms that it is a very special distribution given by the boundary value of an analytic function and it is this analytic function that can be interpreted as the analytic continuation (Wick rotation) of the distribution $G(x)$. On the other hand, as you noted, it satisfies $\Delta G(x) = \delta(x)$, and not $\Delta G(x) = 0$.
On the other hand, the distributional properties of $W(x) = \langle 0 | \phi(x) \phi(0) | 0 \rangle$ are different and don't allow such a simple notion of analytic continuation. I think you could represent $W(x)$ as a linear combination of boundary values of analytic functions, but with the limits taken in particular ways. Such a decomposition $W(x) = \sum_i W_i(x)$, where each $W_i(x)$ is the boundary value of an analytic function, may be possible but is then not required to satisfy $\Delta W_i(x) = 0$ for each individual term. The same goes for the operator valued distribution $\phi(x)$ itself.