Is there any physical meaning for such a correlation function? Consider a thermal scalar field theory, we have the partition functional
$$Z={\rm tr}(e^{-\beta H}).$$
We can build this theory as an Euclidean quantum field theory
$$Z=\int\mathcal{D}\Phi\,e^{-S_E[\Phi]},$$
where we have the imaginary time $\tau\in[0,\beta]$. This imaginary time is in some sense an auxiliary concept. 
I wonder, does the following correlation function
$$\langle \Phi(\tau_1,\vec{x})\,\Phi(\tau_2,\vec{x})\rangle$$
have any physical meaning? A correlation at different imaginary times? 
 A: A Euclidean correlation function may be interpreted as a Lorentzian expectation value by "cutting" the path integral and continuing the time coordinate. 
Let me review how this procedure relates Euclidean correlation functions on a closed manifold $M\times S^1_\beta$ to thermal expectation values in a Lorentzian quantum field theory on $M \times \mathbb{R}_t$. Take $M = S^1_l$ for simplicity. Then the claim is that Euclidean correlation functions on $T^2=S^1_l\times S^1_\beta$ compute expectation values of local operators on $S^1_l\times \mathbb{R}_t$ in a thermal state.
Begin with the Euclidean path integral on $T^2$ with no operator insertions. Imagine cutting the torus on a circular cross-section at fixed Euclidean time, say $\tau = 0$. We obtain a cylinder of length $\beta$ bounded by two circles at $\tau = 0$ and $\beta$. This path integral prepares the thermal density matrix $e^{-\beta H}$, which takes a state in $\mathcal{H}(S^1_l)$ and evolves it forward in Euclidean time by $\beta$. The thermal partition function of the QFT on $S^1_l \times \mathbb{R}_t$, $Z(\beta) = \mathrm{tr}_{\mathcal{H}(S^1_l)} e^{-\beta H}$, is then computed by gluing the two ends of the path integral together (due to the trace), which simply produces the $T^2$ Euclidean partition function:
$$ \mathrm{tr}_{\mathcal{H}(S^1_l)} e^{-\beta H} = \sum_\phi \langle \phi | e^{-\beta H} | \phi \rangle = \int_{T^2} \mathcal{D} \Phi\, e^{-S[\Phi]}.$$
Now insert some local operators $O_i(\theta_i, 0)$ at positions $\theta_i$ on the circle at $t_i = \tau_i = 0$ where we cut the torus. Then by the same reasoning as above, the Euclidean path integral on $T^2$ with these operator insertions computes the thermal expectation value of the operators in the QFT on $S^1_l\times \mathrm{R}_t$:
$$\mathrm{tr}_{\mathcal{H}(S^1_l)}\left( e^{-\beta H} O_1(\theta_1,0)\cdots O_n(\theta_n,0) \right) = \int_{T^2} \mathcal{D} \Phi\, e^{-S[\Phi]}O_1(\theta_1,0)\cdots O_n(\theta_n,0).$$
Finally, we can compute thermal expectation values of operators $O_i(\theta_i,t_i)$ at different Lorentzian times $t_i$ by continuing the Euclidean correlator on $T^2$ to $\tau_i = i t_i$:
$$\langle  O_1(\theta_1,t_1)\cdots O_n(\theta_n,t_n)\rangle_\beta = \int_{T^2} \mathcal{D} \Phi\, e^{-S[\Phi]}O_1(\theta_1,it_1)\cdots O_n(\theta_n,it_n).$$
In other words, we compute the Euclidean correlation function $\langle  O_1(\theta_1,\tau_1)\cdots O_n(\theta_n,\tau_n)\rangle$ on $T^2$ and interpret the answer as a function of complex variables $\tau_i$. Then to compute the Lorentzian expectation value we evaluate the answer on $\tau_i = i t_i$. 
We have considered here Euclidean correlation functions on $T^2$ since your question was concerned with thermal expectation values. But this procedure is completely general. If we know how to compute Euclidean correlation functions on a closed manifold $\mathcal{M}$, then we can cut $\mathcal{M}$ on a slice $M$ at fixed Euclidean time $\tau = 0$. Then by continuing the Euclidean correlation function to Lorentzian time, we can interpret it as an expectation value in a state of the Lorentzian theory on $M\times \mathbb{R}_t$. 
Formally, the expectation value $\langle \Psi_L |\mathcal{T} O_1(x_1,t_1)\cdots O_n(x_n,t_n) | \Psi_R \rangle$ corresponds to a path integral contour in the complex $t$-plane. First the state $|\Psi_R\rangle$ at $t=\tau = 0$ is prepared by a Euclidean path integral in the imaginary time direction. Then the contour turns and runs up along the real (Lorentzian) time direction where the operators are inserted, then turns back and returns to $t=0$, where the Euclidean path integral that prepares the state $\langle \Psi_L |$ is glued on.
