Two point function of a free scalar field in Euclidean space-time 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.
 A: 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.
A: I would like to answer my own question. I thank Igor Khavkine for his answer and other people for their comments which helped me to clarify the question and isolate the main difficulty. 
The answer to the above question is that in Euclidean space-time the non-time-ordered (Wightman) two point function $D(\vec x,\tau):=\langle 0|\phi(\vec x,\tau)\phi(\vec 0,0)|0\rangle$ is not well defined (badly diverges) for $\tau<0$. For $\tau>0$ it is well defined and of course coincides in that domain with the time-ordered two point function. (The same is true for any derivatives of $\phi$.)
Below I will first illustrate this claim by a non-rigorous very general argument which works for not necessarily free fields, and then give explicit computation for free scalar fields.
(1) Let $\phi$ be a field, not necessarily free or scalar. The non-time-ordered product is
$$\langle 0|\phi(\vec x,\tau)\phi(\vec 0,0)|0\rangle=\langle 0|\phi(\vec x,0)e^{-\tau H}\phi(\vec 0,0)|0\rangle.$$
The spectrum of the Hamiltonian $H$ is bounded from below and unbounded from above. Generically the vector $\phi(\vec x,0)|0\rangle$
contains states of arbitrarily high energy. Thus the vector $e^{-\tau H}\phi(\vec 0,0)|0\rangle$ becomes of very infinite norm for $\tau<0$, and the two point function should badly diverge.
(2) Let now $\phi$ be a free scalar field of mass $m\geq 0$ in Euclidean space-time of dimension $D=d+1$. We exclude from the discussion the case when simultaneously $D=2$ and $m=0$. While the latter case is most important for CFT technically it requires separate treatment and will be briefly discussed below in (5).   We denote by $(\vec x,\tau)$ coordinates on Euclidean space and $(\vec x,t)$ coordinates on the corresponding Minkowski space. The Wick rotation is $t=-i\tau$. On the Minkowski space-time the field $\phi(\vec x,t)$ is quantized as follows (see e.g. Peskin-Schroeder book "An introduction to QFT", Ch.2):
\begin{eqnarray}\label{1}
\phi(\vec x,t)=\int\frac{d^d\vec p}{(2\pi)^d}\frac{1}{\sqrt{2E_p}}\left(a_pe^{-ip\cdot x}+a_p^\dagger e^{ip\cdot x}\right)\big|_{p_0=E_p},
\end{eqnarray}
where $E_p=\sqrt{|\vec p|^2+m^2}$ and
\begin{eqnarray*}
[a_p,a_q]=[a_p^\dagger,a_q^\dagger]=0,\, [a_p,a_q^\dagger]=(2\pi)^d\delta^{(d)}(\vec p-\vec q),\\
a_p|0\rangle=0.
\end{eqnarray*}
The Euclidean field $\phi(\vec x,\tau)$ is obtained by substituting $t=-i\tau$. Making that and using the above commutation relations we obtain
\begin{eqnarray}\label{E:DE}
D(\vec x,\tau)=\int\frac{d^dp}{(2\pi)^d}\frac{1}{2E_p}e^{-E_p\tau+i\vec p\cdot \vec x}.(*)
\end{eqnarray}
It is obvious that for $Re(\tau)>0$ the integral absolutely converges while for $Re(\tau)<0$ it badly diverges. Notice that for purely imaginary $\tau$ it still makes sense and is the non-time-ordered two point function in the Minkowski space-time.
(3) One can show that the function $D(\vec x,\tau)$ in the formula ($*$) in the half-space $\tau>0$ is $SO(D)$-invariant, i.e. depends only on $|\vec x|^2+\tau^2$ (e.g. for $D=2$ it is shown in Rotation invariance of an integral ).
Let us show that this implies that the scalar free field does not admit a quantization which is equivariant under the whole Poincare group, i.e. all orientation preserving isometries of $\mathbb{R}^D$. This is in sharp contrast with the Minkowski case, and I guess this is what was meant by Igor Khavkine when he said that there is no operator formalism in Euclidean case. Indeed otherwise for any given point in $x\in\mathbb{R}^D\backslash\{0\}$ we always could choose a time direction such that such that the time coordinate would be positive: $x=(\vec x,\tau), \tau>0$. Then $\langle 0|\phi(x)\phi(0)|0\rangle$ would be well convergent, this is a contradiction. However my impression is that if one fixes the time direction then apparently the free scalar field can be quantized in a way equivariant under the smaller group of all translations and rotations preserving this fixed time direction.
(4) The formula ($*$) defines the function $D$ which is $SO(D)$-invariant in the half-space $\tau>0$ and clearly satisfies $(\Delta-m^2)D=0$. By rotation invariance we can extend $D$ uniquely to $\mathbb{R}^D\backslash\{0\}$. This extension satisfies the same equation. One may ask what is the physical meaning of this extended function? Obviously it is the time-ordered two point function which well converges everywhere but $0$.
(5) Let us discuss briefly the case $D=2,m=0$. It is clear that the integral ($*$) has a divergency at $p=0$ in this case (as opposed to all previously discussed cases). Already in the Minkowski 2d space-time there is an (infrared) divergency which makes the above procedure of quantization impossible. It was shown by Wightman ("Introduction to some aspects of quantizes fields", in "Lectures notes, Cargese Summer School, 1964") that one has to abandon some of the usual requirements imposed on quantization; he has shown that it is possible to quantize the free massless scalar field in 2d Minkowski space-time provided one abandons the assumption of positivity of the scalar product in Hilbert space. I did not try to check how to apply the Wick rotation as in (2) to Wightman's procedure of quantization.
In CFT one works with the periodic boundary conditions: $\phi(x+L,\tau)=\phi(x,\tau)$. Isometries (and even conformal transformations) of the cylinder are just translations preserving the time direction. Thus the possibility of operator formalism of quantization does not contradict to what was said by Khavkine due to the discussion in (3). 
