On "duality between Euclidean time and finite temperature (e.g. QFT and thermal-statistical mechanics)", I will do three parts. The first part describes the duality between the Feynman path integral and thermal density matrix equations. The second part is a review on how density matrix cooling works and the third part is how one might get something useful out of density matrix cooling.
Analytic Continuation of the Feynman Path Integral as a Duality Transform
A. Zee says in the book of his you quoted, in his Chapter V.2 "Euclid, Boltzmann, Hawking and Field Theory at Finite Temperature", page 263 in 2003 first edition:
Surely you would hit it big with mystical types if you were to tell
them that temperature is equivalent to cyclic imaginary time. At the
arithmetic level, this connection comes merely from the fact that the
central objects in quantum physics $e^{-iHt}$ and in thermal physics
$e^{-\beta H}$ are related by analytic continuation. Some physicists,
myself included, feel that there may be something profound here that
we have not quite understood.
The concept of duality is that some physics problems can be converted to a different form where they become much easier (or harder). The dualities discussed in the literature at present are generally between two QFTs or between a QFT and a string theory. But the analytic continuation Zee refers to is between QFT and quantum statistical mechanics. So to use it as a duality, we have to put statistical mechanics into a quantum form.
The $e^{-\beta H}$ Zee mentions defines the thermal density matrix $\rho(T)$ by:
\begin{equation}
\rho(T) \propto e^{-H/T}
\end{equation}
where the proportionality constant is $Z = \textrm{tr}(e^{-\beta H})$ and I've chosen units so that Boltzmann's constant is unity so that $\beta = 1/(k_BT) = 1/T$. Since a density matrix is a type of quantum mechanics, this makes the analytic continuation Zee discusses one between the Feynmann path integral, which is a formulation of QFT and the density matrix, which is a formulation of QM.
In using a duality one hopes to convert a difficult quantum problem to an easy one. The difficulty is in showing that the two problems are related. How can these two be related? QFT uses state vectors $|k\rangle$ to represent states; density matrices can be made from these by $\rho = |k\rangle\langle k|$, so that's a start. Quantum field theory is a type of state vector quantum mechanics where the vector components are occupation numbers for each possible energy. So these formulations of quantum mechanics are not entirely unrelated and we can at least imagine that they are related by a duality transformation.
So what kinds of elementary particle problems can we apply density matrices to? The simplest place to look is in theories that do not involve space or time dependency. These kinds of problems seem like they could get you nowhere but you can sometimes get information about the particle types without dealing with spacetime explicitly. For example, you do not need to solve the Dirac equation to realize that its 4x4 matrices imply spin-1/2 electrons and positrons.
The immediate result of the analytic continuation of the Feynmann path integral is to convert its computation from small slices of time to small slices of $\beta = 1/T$. This means that one converts small progress (propagation) in time into small progress in inverse temperature. In other words, what we did in QFT to make things advance in time becomes, in thermal density matrices, a method of reducing the temperature of density matrices.
Cooling Density Matrices
Going to single point spacetime so that the integrals over space reduce to matrix multiplication, the statistical mechanical equivalent of the Feynmann path integral becomes:
\begin{equation}
\rho(T/N) \propto (\rho(T))^N.
\end{equation}
Again the missing / implied proportionality constants are $Z(T/N)$ on the left and $Z(T)^N$ on the right. This defines a method for halving the temperature of density matrix: Square the matrix and then fix its normalization by dividing by the trace. As you repeat the operation the density matrix gets colder and colder until it approaches a limit where $\rho^2=\rho$ and you can make no further progress. This is the cold temperature limit; the result are the pure density matrices.
For spin-1/2, the mixed (finite temperature) density matrices make up the interior of a sphere. The pure states are on the surface of the sphere; given a real unit vector $\vec{u} = (u_x,u_y,u_z)$ the pure state for spin-1/2 in the direction $\vec{u}$ is given by
\begin{equation}
\rho_u = 0.5(1 + u_x\sigma_x + u_y\sigma_y + u_z\sigma_z) = (1+\vec{u}\cdot \vec{\sigma})/2.
\end{equation}
A mixed density matrix can be obtained by making $|\vec{u}|^2 < 1$. So a finite temperature density matrix is given by
\begin{equation}
\rho_u(T) =(1+0.2\vec{u}\cdot \vec{\sigma})/2.
\end{equation}
Squaring $\rho_u(T)$ we use $(\vec{u}\cdot\vec{\sigma})^2=1$ to get
\begin{equation}
(\rho_u(T))^2 =(1+2(0.2)\vec{u}\cdot \vec{\sigma}+(0.2)^2)/4
=(0.52 + 0.2\vec{u}\cdot \vec{\sigma} )/2
\end{equation}
which has trace $0.52$, so to normalize it we divide by the trace and get
\begin{equation}
\rho_u(T/2) =(1+0.3846\;\vec{u}\cdot \vec{\sigma})/2
\end{equation}
and we see that halving the temperature increased the $\vec{u}\cdot \vec{\sigma}$
coefficient from $0.2$ to $0.3846$, that is, it almost doubled the spin part. Continuing to do this keeps increasing this coefficient until it approaches 1 and we have the pure state $\rho_u$.
Density Matrix Symmetry from a Finite Group Algebra
The duality described here is useless until we make additional assumptions. Since the "easy" side of the calculation is density matrix cooling, our assumptions will have to apply to density matrices.
If there are no restrictions on the density matrices, so they can be any NxN matrix that has unit trace and is Hermitian, then their symmetry will be the full SU(N). For example, the pure 2x2 density matrices are made from SU(2) state vectors and so inherit SU(2) symmetry. More generally, general unrestricted density matrices from NxN matrices will transform under the fundamental representation of SU(N) symmetry.
One of the things that density matrices allow that cannot be done with state vectors is an easy representation of statistical mixtures. That is, mixtures that cannot exist in quantum superposition. For example (at the cold temperature limit), it is not possible to create a quantum superposition between particles with different electric charges. They are in different "superselection sectors".
The Hamiltonian has to commute with the charge operator and this divides up the Hilbert space into sectors where it's impossible to get from one to another. But if one has a beam of particles that is half electrons and half neutrinos one can model this as a mixed density matrix. If we use SU(2) to represent the spin-1/2 of the electron and neutrino, the superselection sector restriction for state vectors becomes a restriction on the mixed density matrix that the cross product terms between the neutrino and electron have to be zero. In other words the density matrix has to be block diagonal:
\begin{equation}
\rho = \left(\begin{array}{cccc}
\nu_{11}&\nu_{12}&&\\
\nu_{21}&\nu_{22}&&\\
&&e_{11}&e_{12}\\
&&e_{21}&e_{22}
\end{array}\right)
\end{equation}
where the matrix entries that have to be zero are left blank.
Now it turns out that block diagonal algebras can be created from finite groups as is discussed at length in Hamermesh's book "Group Theory and its Application to Physical Problems". I've got the paperback 1989 Dover edition list price \$22.95. He's considering a finite group $G$ with elements $R$:
$<<<$
Starting from a group $G$, we can construct a system which is called an algebra. The quantities in the algebra are
\begin{equation}
\Sigma_R\;a_R\;R,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(3-162)
\end{equation}
where the coefficients $a_R$ are any complex numbers. By the sum of two quantitites we understand
\begin{equation}
\Sigma_Ra_RR + \Sigma_Rb_RR = \Sigma_R(a_R+b_R)R,
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(3-163)
\end{equation}
and the product implies
\begin{equation}
\Sigma_Ra_RR \cdot \Sigma_Rb_RR = \Sigma_{R,S}a_Rb_S\;RS
=\Sigma_R(\Sigma_S a_{RS}b_{S^{-1}} )R.
\;\;\;\;\;\;(3-164)
\end{equation}
$>>>$
The finite group algebra over the complexes is used in Hamermesh and similar physics group theory texts as a way of defining the irreducible representations of the finite group symmetry. This is because this algebra happens to have every possible irreducible representation included, and each included just one time. Very convenient when you're trying to find all the irreducible representations of a finite group symmetry!
When you diagonalize the group algebra as much as possible, you end up with a block diagonal representation. Each block corresponds to a distinct irreducible representation of the group. And what's more, the complex degrees of freedom in that block diagonal representation correspond exactly to the complex degrees of freedom in the group algebra.
For example, the finite group $S_3$ are the 6 permutations on 3 objects. So the finite group is 6 dimensional. Now an NxN diagonal block will use up $N^2$ complex degrees of freedom and since $S_3$ is not Abelian it will have to have at least one block bigger than 1x1. There's only one way to write 6 as a sum of squares under that restriction, $6 = 1^2 + 1^2 +2^2$ so there are three irreducible representations of $S_3$; two that are 1x1 and one that is 2x2. Thus the finite group algebra gives a convenient way of listing the possible irreducible representations of a finite group.
But from the point of view of density matrices (instead of state vectors) we can think of the algebra as where the mixed particle density matrices live. If we began with a unmixed initial density matrix, say with nonzero entries only in the 2x2 diagonal block, it's clear that the cold temperature limit will be in that block and so will correspond to a particle with internal SU(2) symmetry. And if the initial condition was nonzero only in one of the 1x1 blocks that would be the cold temperature limit. For cases that are mixed, the winner will be the block that dominates the others so the cold temperature limit has three particles defined, one with SU(2) internal symmetry and two with no internal symmetry. We would call them an SU(2) doublet and two SU(2) singlets just as the quarks are two SU(3) color triplets and the leptons are two SU(3) color singlets.
So a possible application of the QFT / thermal statistics analytic continuation duality is a way of defining particle symmetry and representations by a choice of a finite group.
First of all, note that one cannot associate a temperature to a single quantum state (cf "vacuum state of the theory is defined as having zero energy and zero temperature"), and having a zero energy vacuum state is just a convention (as it is cut-off dependent, and thus renormalized).
Furthermore, the OP is confused. Standard (i.e. zero-temperature) QFTs compute observables in the interacting ground-state of the theory, call it $|\Omega\rangle$, which is different from the non-interacting ground-state $|0\rangle$. After that, all the discussion about normal ordering and in and out state in the OP's question is "just" related to technicalities to construct $|\Omega\rangle$ from $|0\rangle$.
Thermal field theory assumes that the observables are obtained from a system the state of which is a thermal distribution, i.e. described by a density matrix $\rho= \exp(-\beta H)/Z$ with $\beta$ the inverse temperature. In the eigenbasis $|n\rangle$ of the interacting hamiltonian (with $|n=0\rangle =|\Omega\rangle$ the same ground-state as above), a single time observable is given by
$$\langle A\rangle = \sum_n \frac{e^{-\beta E_n}}{Z}\langle n|A|n\rangle .$$
In the limit $\beta\to\infty$, one recovers the zero temperature QFT result $\langle A\rangle=\langle \Omega|A|\Omega\rangle$.
Of course, one need to devise a perturbation scheme to relate the interacting eigenstate to the non-interacting eigenstate, more or less copied to that of the zero-temperature QFTs.
Best Answer
Here is a proof following Ojima, "Lorentz Invariance vs. Temperature in QFT", Letters in Mathematical Physics (1986) Vol. 11, Issue 1 (1986) 73-80. The first two pages of the paper are available for free here, but the website wants money for more of the paper. (Click the orange "Look Inside" button if the paper doesn't open automatically.) Fortunately, the proof is on the second page.
Define $$w(A) = tr\bigl(e^{-\beta H} A\bigr)/tr\bigl(e^{-\beta H}\bigr).$$ The KMS condition can be written $$w(\phi(x)\phi(y)) = w(\phi(y)\phi({\tilde x}))$$ where $\tilde x$ is $x$ with the time component shifted by $i\beta$.
Now consider the Fourier transform $$\langle\phi_k\phi_{-k}\rangle = \int d^4x d^4y\ e^{i k \cdot (x-y)} w(\phi(x)\phi(y)).$$ By the KMS condition this is $$\langle\phi_k\phi_{-k}\rangle = \int d^4x d^4y\ e^{i k \cdot (x-y)} w(\phi( y)\phi({\tilde x})).$$ Shifting the time in the $x$ integral, this becomes $$e^{\beta k_0}\int d^4{\tilde x} d^4 y\ e^{i k \cdot ({\tilde x}- y)} w(\phi(y)\phi({\tilde x})) = e^{\beta k_0}\langle\phi_{-k}\phi_{k}\rangle,$$ so we have $$\langle\phi_k\phi_{-k}\rangle = e^{\beta k_0}\langle\phi_{-k}\phi_{k}\rangle.$$ Starting with the right hand side of this equation, if Lorentz invariance holds, $\langle\phi_{-k}\phi_{k}\rangle$ is a scalar under Lorentz transformations for a scalar field $\phi$. However, $e^{\beta k_0}$ is manifestly not a Lorentz scalar since it depends non-covariantly on $k_0$. This implies that the left hand side of the above equation is not a Lorentz scalar, contradicting the Lorentz invariance of $\langle\phi_k\phi_{-k}\rangle$.
This proves that Lorentz covariance cannot hold for finite $\beta$.
Another way to see that Lorentz covariance is broken at finite temperature is to Wick rotate to Euclidean spacetime. The KMS condition then implies periodicity of the Green's function in the time direction. Lorentz transformations in real spacetime are mapped to rotations in Euclidean spacetime by the Wick rotation. Since periodic boundary conditions are imposed in one direction and not the other three directions, rotation symmetry is broken in Euclidean spacetime and therefore Lorentz symmetry is broken in real spacetime.