text:
# References

Relativistic Doppler Effect in Quantum Communication
ASHER PERES and DANIEL R. TERNO

Department of Physics, Technion—Israel Institute of Technology, 32000 Haifa, Israel

Abstract.  When an electromagnetic signal propagates in vacuo, a polarization detector cannot be rigorously perpendicular to the wave vector because of diffraction effects. The vacuum behaves as a noisy channel, even if the detectors are perfect. The “noise” can however be reduced and nearly cancelled by a relative motion of the observer toward the source. The standard definition of a reduced density matrix fails for photon polarization, because the transversality condition behaves like a superselection rule. We can however define an effective reduced density matrix which corresponds to a restricted class of positive operator-valued measures. There are no pure photon qubits, and no exactly orthogonal qubit states.

1. Introduction

The long range propagation of polarized photons is an essential tool of quantum cryptography . Usually, optical fibers are used, and the photons may be absorbed or depolarized due to imperfections. In some cases, such as communication with space stations, the photons must propagate in vacuo . The beam then has a finite diffraction angle of order $`\lambda /a`$, where $`a`$ is the aperture size, and new deleterious effects appear. In particular a polarization detector cannot be rigorously perpendicular to the wave vector, and the transmission is never faithful, even with perfect detectors. Moreover, the “vacuum noise” depends on the relative motion of the observer with respect to the source.

The relativistic effects reported here are essentially different from those for massive particles  because massless particles have only two linearly independent polarization states. The properties that we discuss are kinematical, not dynamical. At the statistical level, it is not even necessary to involve quantum electrodynamics. Most formulas can be derived by elementary classical methods as shown below. It is only when we need to consider individual photons, for cryptographic applications, that quantum theory becomes essential.

This article consists of two parts. First we consider the propagation of a classical electromagnetic wave. The wave vector cannot be constant because of diffraction effects. A polarization detector cannot unambiguously distinguish orthogonal polarizations, even if the detector is perfect. The vacuum behaves as a noisy channel. We then show that this “noise” can be reduced and nearly cancelled by a relative motion of the observer toward the source.

In the second part of this paper, we investigate the transmission of a single photon. The diffraction effects mentioned above lead to superselection rules which make it impossible to define a reduced density matrix for polarization. It is still possible to have “effective” density matrices; however, the latter depend not only on the preparation process, but also on the type of detection that is used by the observer.

2. Classical electromagnetic signals

Assume for simplicity that the electromagnetic signal is monochromatic. In a Fourier decomposition, the Cartesian components of the wave vector $`k_\mu `$ (with $`\mu =0,1,2,3`$) can be written in term of polar angles:

$$k_\mu =(1,\mathrm{sin}\theta \mathrm{cos}\varphi ,\mathrm{sin}\theta \mathrm{sin}\varphi ,\mathrm{cos}\theta ),$$
(1)
where we use units such that $`c=1`$ and $`k_0=1`$. Let us choose the $`z`$ axis so that a well collimated beam has a large amplitude only for small $`\theta `$, and let us rotate the $`x`$ and $`y`$ axes so that a particular k in which we are interested has $`\varphi =0`$ (we shall later return to arbitrary k with $`\varphi 0`$). The Fourier transform of the electric field is perpendicular to

$$𝐤=(\mathrm{sin}\theta ,0,\mathrm{cos}\theta ).$$
(2)
If the emitter (“Alice”) selects a linear polarization angle $`a`$, the Fourier transform of E is proportional to

$$𝐄_𝐤=(\mathrm{cos}a\mathrm{cos}\theta ,\mathrm{sin}a,\mathrm{cos}a\mathrm{sin}\theta ).$$
(3)
The magnetic field Fourier transform is proportional to $`𝐁_𝐤=𝐤\times 𝐄_𝐤`$. The Poynting vector P is parallel to k and gives the energy flux, as usual. We shall henceforth omit the subscript k.

Suppose that the receiver (“Bob”) has an infinite flat detector parallel to the $`xy`$ plane. Then only the component

$$P_zE_x^2+E_y^2=\mathrm{cos}^2a\mathrm{cos}^2\theta +\mathrm{sin}^2a$$
(4)
is absorbed by the detector. Moreover, if Bob selects a polarization angle $`b`$ in the $`xy`$ plane, the flux detected by him will be proportional to

$$(E_x\mathrm{cos}b+E_y\mathrm{sin}b)^2=(\mathrm{cos}a\mathrm{cos}b\mathrm{cos}\theta +\mathrm{sin}a\mathrm{sin}b)^2.$$
(5)
Under ideal conditions ($`\theta =0`$), the maximal signal behaves as $`\mathrm{cos}^2(ab)`$. For small but finite $`\theta `$, the fraction of the signal that is lost, when $`a=b`$, is $`\theta ^2\mathrm{cos}^2a`$. More generally it is $`\theta ^2\mathrm{cos}^2(a\varphi )`$, where $`\varphi `$ is the azimuthal angle of k, which was set to zero in the preceding calculation by a suitable choice of the $`x`$ axis.

Apart from the above loss in intensity, the angle mismatch may introduce errors. Consider again the case $`\varphi =0`$. When Alice emits a wave that is linearly polarized along the $`x`$ or $`y`$ axes, Bob’s intensity pattern varies as $`A_x^2\mathrm{cos}^2b`$ and $`A_y^2\mathrm{sin}^2b`$, respectively. These intensity distributions would be unambiguously recognized as pertaining to the corresponding linear polarizations. However, for a general linear polarization state that Alice may send ($`0a\pi /2`$) the intensity will be distributed according to Eq. (5). As a result, the angle $`a^{}`$ that Bob will ascribe to it would be related to the actual polarization direction by $`\mathrm{tan}a^{}=\mathrm{tan}a/\mathrm{cos}\theta `$.

In a real experiment, the angles $`\theta `$ and $`\varphi `$ are distributed in a continuous way around the $`z`$ axis (exactly how depends on the properties of the laser) and one has to take a suitable average over them. Since the definition of polarization explicitly depends on the direction of k, taking the average over many values of k leads to an impure polarization and therefore may cause not only an attenuation of the beam, but also identification errors.

Let us now consider the effect of a motion of Bob relative to Alice, with a constant velocity $`𝐯=(0,0,v)`$. The Lorentz transformation of $`k_\mu `$ in Eq. (1) yields new components

$$k_0^{}=\gamma (1v\mathrm{cos}\theta )\text{and}k_z^{}=\gamma (\mathrm{cos}\theta v),$$
(6)
where $`\gamma =(1v^2)^{1/2}`$. Considering again a single Fourier component, we have, instead of the unit vector k, a new unit vector

$$𝐤^{}=(\frac{\mathrm{sin}\theta }{\gamma (1v\mathrm{cos}\theta )},\mathrm{\hspace{0.33em}0},\frac{\mathrm{cos}\theta v}{1v\mathrm{cos}\theta }).$$
(7)
In other words, there is a new tilt angle $`\theta ^{}`$ given by

$$\mathrm{sin}\theta ^{}=\mathrm{sin}\theta /\gamma (1v\mathrm{cos}\theta ).$$
(8)
For small $`\theta `$, such that $`\theta ^2|v|`$, we have

$$\theta ^{}=\theta \sqrt{\frac{1+v}{1v}}.$$
(9)
The square root is the familiar relativistic Doppler factor. For large negative $`v`$, the diffraction angle becomes arbitrarily small, and sideway losses can be reduced to zero.

It is noteworthy that the same Doppler factor was obtained by Jarett and Cover  who considered only the relativistic transformations of bit rate and noise intensity, without any specific physical model. This remarkable agreement shows that information theory should properly be considered as a branch of physics.

3. Quantized electromagnetic signals

We now turn to the quantum language in order to discuss applications to secure communication. The ideal scenario is that Alice sends isolated photons (one particle Fock states). In a more realistic setup, the transmission is by means of weak coherent pulses containing on the average less than one photon each.

A basis of the one-photon space is spanned by states of definite momentum and helicity,

$$|𝐤,\mathit{ϵ}_𝐤^\pm |𝐤|\mathit{ϵ}_𝐤^\pm ,$$
(10)
where helicity states $`|\mathit{ϵ}_𝐤^\pm `$ are explicitly defined in Eq. (15) below. The momentum basis is normalized by

$$𝐪|𝐤=(2\pi )^3(2k^0)\delta ^{(3)}(𝐪𝐤).$$
(11)
Polarization states that correspond to different momenta belong to distinct Hilbert spaces and cannot be superposed (an expression such as $`|\mathit{ϵ}_𝐤^\pm +|\mathit{ϵ}_𝐪^\pm `$ is meaningless if $`𝐤𝐪`$). The complete basis (10) does not violate this superselection rule, owing to the othogonality of the momentum basis. Therefore, a generic one-photon state is given by a wave packet

$$|\mathrm{\Psi }=𝑑\mu (𝐤)f(𝐤)|𝐤,𝜶(𝐤),$$
(12)
where the polarization state $`|𝜶(𝐤)`$ corresponds to the 3-vector

$$𝜶(𝐤)=\alpha _+(𝐤)\mathit{ϵ}_𝐤^++\alpha _{}(𝐤)\mathit{ϵ}_𝐤^{},$$
(13)
$`|\alpha _+|^2+|\alpha _{}|^2=1`$, and the explicit form of $`\mathit{ϵ}_𝐤^\pm `$ is given below. The Lorentz–invariant measure is

$$d\mu (k)=\frac{d^3𝐤}{(2\pi )^32k^0},$$
(14)
and normalized states satisfy $`𝑑\mu (k)|f(𝐤)|^2=1`$. Since diffraction angles and frequency spreads are usually small , $`f(𝐤)`$ significantly differs from zero only in the vicinity of a certain momentum that we shall denote by $`𝐤_A`$.

Lorentz transformations of quantum states are most easily computed by referring to some standard momentum, which for photons is $`p^\nu =(1,0,0,1)`$. Accordingly, standard right and left circular polarization vectors are $`\mathit{ϵ}_p^\pm =(1,\pm i,0)/\sqrt{2}`$. If we are interested in linear polarization, all we have to do is to use Eq. (13) with $`\alpha _+=(\alpha _{})^{}`$, so that $`𝜶(𝐤)`$ is real. In general, $`𝜶(𝐤)`$ corresponds to elliptic polarization.

Under a Lorentz transformation $`\mathrm{\Lambda }`$, these states become $`|𝐤_\mathrm{\Lambda },𝜶(𝐤_\mathrm{\Lambda })`$, where $`𝐤_\mathrm{\Lambda }`$ is the spatial part of a four-vector $`k_\mathrm{\Lambda }=\mathrm{\Lambda }k`$, and the new polarization vector can be obtained by an appropriate rotation that is described below. For each $`𝐤`$ a polarization basis is labeled by the helicity vectors,

$$\mathit{ϵ}_𝐤^\pm =R(\widehat{𝐤})\mathit{ϵ}_p^\pm .$$
(15)
and the corresponding quantum states are just $`|𝐤,\mathit{ϵ}_𝐤^\pm `$. As usual, $`\widehat{𝐤}`$ denotes the unit 3-vector in the direction of $`𝐤`$. The matrix that rotates the standard direction $`(0,0,1)`$ to $`\widehat{𝐤}=(\mathrm{sin}\theta \mathrm{cos}\varphi ,\mathrm{sin}\theta \mathrm{sin}\varphi ,\mathrm{cos}\theta )`$ is

$$R(\widehat{𝐤})=\left(\begin{array}{ccccc}\mathrm{cos}\theta \mathrm{cos}\varphi & & \mathrm{sin}\varphi & & \mathrm{cos}\varphi \mathrm{sin}\theta \\ \mathrm{cos}\theta \mathrm{sin}\varphi & & \mathrm{cos}\varphi & & \mathrm{sin}\varphi \mathrm{sin}\theta \\ \mathrm{sin}\theta & & 0& & \mathrm{cos}\theta \end{array}\right),$$
(16)
and likewise for $`\widehat{𝐤}_\mathrm{\Lambda }`$.

Under a general Lorentz transformation, be it a rotation or a boost, helicity is preserved, but the states (and corresponding geometric vectors) acquire helicity-dependent phases,

$$\alpha _+\mathit{ϵ}_𝐤^++\alpha _{}\mathit{ϵ}_𝐤^{}\alpha _+e^{i\xi (\mathrm{\Lambda },\widehat{𝐤})}\mathit{ϵ}_{𝐤_\mathrm{\Lambda }}^++\alpha _{}e^{i\xi (\mathrm{\Lambda },\widehat{𝐤})}\mathit{ϵ}_{𝐤_\mathrm{\Lambda }}^{},$$
(17)
where the explicit expressions for $`\xi (\mathrm{\Lambda },\widehat{𝐤})`$ are given in  The superselection rule that was mentioned above makes it impossible to define a reduced density matrix in the usual way. We can however define an “effective” reduced density matrix for polarization, as follows. The labelling of polarization states by Euclidean vectors $`𝐞_𝐤^n`$, and the fact that photons are spin-1 particles, suggest the use of a $`3\times 3`$ matrix with entries labelled $`x`$, $`y`$ and $`z`$. Classically, they correspond to different directions of the electric field. For example, when $`𝐤=k_A\widehat{𝐳}`$, only $`\rho _{xx}`$, $`\rho _{xy}`$, $`\rho _{yy}`$ are non-zero.

For a generic photon state $`|\mathrm{\Psi }`$, let us try to construct a reduced density matrix $`\rho _{xx}`$ that gives the expectation value of an operator representing the polarization in the $`x`$ direction, irrespective of the particle’s momentum. To have a momentum-independent polarization is to tacitly admit longitudinal photons. Therefore, in terms of real transversal photons of a momentum $`𝐤`$, this problem suggests to find the direction that is perpendicular to $`𝐤`$ and closest to an arbitrary unit vector $`\widehat{𝐱}`$. That is, we are looking for a unit Euclidean complex vector $`𝐞_x(𝐤)`$ such that $`𝐤𝐞_x(𝐤)=0`$ and $`\widehat{𝐱}𝐞_x(𝐤)=\mathrm{max}`$.

Momentum-independent polarization states thus consist of physical (transversal) and unphysical (longitudinal) parts with a polarization vector $`\mathit{ϵ}^{\mathrm{}}=\widehat{𝐤}`$. For example, a polarization state along the $`x`$-axis is

$$|\widehat{𝐱}=x_+(𝐤)|\mathit{ϵ}_𝐤^++x_{}(𝐤)|\mathit{ϵ}_𝐤^{}+x_{\mathrm{}}(𝐤)|\mathit{ϵ}_𝐤^{\mathrm{}},$$
(18)
where $`x_\pm (𝐤)=\mathit{ϵ}_𝐤^\pm \widehat{𝐱}`$, and $`x_{\mathrm{}}(𝐤)=\widehat{𝐱}\widehat{𝐤}=\mathrm{sin}\theta \mathrm{cos}\varphi `$. It follows that $`|x_+|^2+|x_{}|^2+|x_{\mathrm{}}|^2=1`$, and

$$𝐞_x(𝐤)=\frac{x_+(𝐤)\mathit{ϵ}_𝐤^++x_{}(𝐤)\mathit{ϵ}_𝐤^{}}{\sqrt{x_+^2+x_{}^2}}.$$
(19)
Note that $`\widehat{𝐱}|\widehat{𝐲}=\widehat{𝐱}\widehat{𝐲}=0`$, whence

$$|\widehat{𝐱}\widehat{𝐱}|+|\widehat{𝐲}\widehat{𝐲}|+|\widehat{𝐳}\widehat{𝐳}|=\text{1l}.$$
(20)
To the direction $`\widehat{𝐱}`$ corresponds a projection operator

$$P_{xx}=|\widehat{𝐱}\widehat{𝐱}|\text{1l}_p=|\widehat{𝐱}\widehat{𝐱}|𝑑\mu (𝐤)|𝐤𝐤|,$$
(21)
where $`\text{1l}_p`$ is the unit operator in momentum space. The action of $`P_{xx}`$ on $`|\mathrm{\Psi }`$ follows from Eq. (18) and $`\mathit{ϵ}_𝐤^\pm |\mathit{ϵ}_𝐤^{\mathrm{}}=0`$. Only the transversal part of $`|\widehat{𝐱}`$ appears in the expectation value:

$$\mathrm{\Psi }|P_{xx}|\mathrm{\Psi }=𝑑\mu (𝐤)|f(𝐤)|^2|x_+(𝐤)\alpha _+^{}(𝐤)+x_{}(𝐤)\alpha _{}^{}(𝐤)|^2.$$
(22)
Define the transversal part of $`|\widehat{𝐱}`$:

$$|𝐛_x(𝐤)(|\mathit{ϵ}_𝐤^+\mathit{ϵ}_𝐤^+|+|\mathit{ϵ}_𝐤^{}\mathit{ϵ}_𝐤^{}|)|\widehat{𝐱}=x_+(𝐤)|\mathit{ϵ}_𝐤^++x_{}(𝐤)|\mathit{ϵ}_𝐤^{}.$$
(23)
Likewise define $`|𝐛_y(𝐤)`$ and $`|𝐛_z(𝐤)`$. These three state vectors are neither of unit length nor mutually orthogonal. For $`𝐤=(\mathrm{sin}\theta \mathrm{cos}\varphi ,\mathrm{sin}\theta \mathrm{sin}\varphi ,\mathrm{cos}\theta )`$ we have

$`|𝐛_x(𝐤)`$ $`=`$ $`[(\mathrm{cos}\theta \mathrm{cos}\varphi +i\mathrm{sin}\varphi )|\mathit{ϵ}_𝐤^++(\mathrm{cos}\theta \mathrm{cos}\varphi i\mathrm{sin}\varphi )|\mathit{ϵ}_𝐤^{}]/\sqrt{2},`$ (24)
$`=`$ $`c(\theta ,\varphi )|𝐤,𝐞_x(𝐤),`$ (25)
where $`𝐞_x(𝐤)`$ is given by Eq. (19), and $`c(\theta ,\varphi )=\sqrt{x_+^2+x_{}^2}`$.

Finally, a POVM element $`E_{xx}`$ which is the physical part of $`P_{xx}`$, namely is equivalent to $`P_{xx}`$ for physical states (without longitudinal photons) is

$$E_{xx}=𝑑\mu (𝐤)|𝐤,𝐛_x(𝐤)𝐤,𝐛_x(𝐤)|,$$
(26)
and likewise for other directions. The operators $`E_{xx}`$, $`E_{yy}`$ and $`E_{zz}`$ indeed form a POVM in the space of physical states, owing to Eq. (20). It then follows from Eq. (23) and similar definitions for other directions that, for any $`𝐤`$,

$$|𝐛_x(𝐤)𝐛_x(𝐤)|+|𝐛_y(𝐤)𝐛_y(𝐤)|+|𝐛_z(𝐤)𝐛_z(𝐤)|=\text{1l}_𝐤,$$
(27)
where $`\text{1l}_𝐤`$ is the identity operator in the subspace of polarizations orthogonal to $`𝐤`$.

To complete the construction of the density matrix, we introduce additional directions. Following a standard practice of state reconstruction , we consider $`P_{x+z,x+z}`$, $`P_{xiz,xiz}`$ and similar combinations. For example,

$$P_{x+z,x+z}=\frac{1}{2}(|\widehat{𝐱}+|\widehat{𝐳})(\widehat{𝐱}|+\widehat{𝐳}|)\text{1l}_p$$
(28)
The diagonal elements of the new polarization density matrix are defined as

$$\rho _{mm}=\mathrm{\Psi }|E_{mm}|\mathrm{\Psi },m=x,y,z,$$
(29)
and the off-diagonal elements are recovered by combinations such as

$$\rho _{xz}=\mathrm{\Psi }(|\widehat{𝐱}\widehat{𝐳}|\text{1l}_p)|\mathrm{\Psi }=\mathrm{\Psi }|E_{x+z,x+z}+E_{xiz,xiz}E_{xx}E_{zz}|\mathrm{\Psi },$$
(30)
where we denote $`|\widehat{𝐱}\widehat{𝐳}|\text{1l}_p`$ as $`P_{xz}`$, and by $`E_{xz}`$ its “physical” part. We then get a simple expression for the reduced density matrix corresponding to the polarization state $`|𝜶(𝐤)`$:

$$\rho _{mn}=\mathrm{\Psi }|E_{mn}|\mathrm{\Psi }=d\mu (𝐤)|f(𝐤)|^2𝜶(𝐤)|𝐛_m(𝐤)𝐛_n(𝐤)|𝜶(𝐤),m,n,=x,y,z.$$
(31)
Our basis states $`|𝐤,\mathit{ϵ}_𝐤`$ are direct products of momentum and polarization. Owing to the transversality requirement $`\mathit{ϵ}(𝐤)𝐤=0`$, they remain direct products under Lorentz transformations. All the other states have their polarization and momentum degrees of freedom entangled. As a result, if one is restricted to polarization measurements as described by the above POVM, there do not exist two orthogonal polarization states. In general, any measurement procedure with finite momentum sensitivity will lead to the errors in identification. This can be seen as follows.

Let two states $`|\mathrm{\Phi }`$ and $`|\mathrm{\Psi }`$ be of the form in Eq. (12). Their reduced polarization density matrices, $`\rho _\mathrm{\Phi }`$ and $`\rho _\mathrm{\Psi }`$, respectively, are calculated using Eq. (31). Since the states are entangled, the von Neumann entropies of the reduced density matrices, $`S=\mathrm{tr}(\rho \mathrm{ln}\rho )`$, are positive . Therefore, both matrices are at least of rank two. Since the overall dimension is 3, it follows that $`\mathrm{tr}(\rho _\mathrm{\Phi }\rho _\mathrm{\Psi })>0`$ and these states are not perfectly distinguishable. An immediate corollary is that photon polarization states cannot be cloned perfectly. This is because no-cloning theorem, in its various versions , forbids an exact copying of unknown non-orthogonal states.

To quantify the distinguishability of a pair of quantum states, we shall use the simplest criterion, namely the probability of error $`P_E`$, defined as follows: an observer receives a single copy of one of the two known states and performs any operation permitted by quantum theory, in order to decide which state was supplied. The probability of a wrong answer for an optimal measurement is

$$P_E(\rho _1,\rho _2)=\frac{1}{2}+\frac{1}{4}\mathrm{tr}|\rho _1\rho _2|,$$
(32)
where, for any operator $`O`$, the operator $`|O|`$ is defined as $`\sqrt{(O^{}O)}`$. As shown below, the distinguishability of polarization density matrices depends on the observers’ motion. First we present some general considerations and then illustrate them with a simple example.

Let us take the $`z`$-axis to coincide with the average direction of propagation so that the mean photon momentum is $`k_A\widehat{𝐳}`$. Typically, the spread in momentum is small, but not necessarily equal in all directions. Usually the intensity profile of laser beams has cylindrical symmetry, and we may assume that $`\mathrm{\Delta }_x\mathrm{\Delta }_y\mathrm{\Delta }_r`$ where the index $`r`$ means radial. We may also assume that $`\mathrm{\Delta }_r\mathrm{\Delta }_z`$, since the momentum spread along the average propagation direction does not necessarily increase the mixness of the reduced density matrix, as we shell see below. We then have

$$f(𝐤)f_1[(k_zk_A)/\mathrm{\Delta }_z]f_2(k_r/\mathrm{\Delta }_r).$$
(33)
We approximate

$$\theta \mathrm{tan}\theta k_r/k_zk_r/k_A.$$
(34)
In pictorial language, polarization planes for different momenta are tilted by angles up to $`\mathrm{\Delta }_r/k_A`$, so that we expect an error probability of the order $`\mathrm{\Delta }_r^2/k_A^2`$. In the density matrix $`\rho _{mn}`$ all the elements of the form $`\rho _{mz}`$ should vanish when $`\mathrm{\Delta }_r0`$. Moreover, if $`\mathrm{\Delta }_z0`$, the non-vanishing $`xy`$ block goes to the usual (monochromatic) polarization density matrix,

$$\rho _{\mathrm{pure}}=\left(\begin{array}{ccc}|\alpha |^2& \beta & 0\\ \beta ^{}& 1|\alpha |^2& 0\\ 0& 0& 0\end{array}\right).$$
(35)
As an example we consider two states which, if the momentum spread could be ignored, would be $`|k_A\widehat{𝐳},\mathit{ϵ}_{k_A\widehat{𝐳}}^\pm `$. To simplify the calculations we assume a Gaussian distribution:

$$f(𝐤)=Ne^{(k_zk_A)^2/2\mathrm{\Delta }_z^2}e^{k_r^2/2\mathrm{\Delta }_r^2},$$
(36)
where $`N`$ is a normalization factor and $`\mathrm{\Delta }_z\mathrm{\Delta }_r`$. Moreover, we take the polarization components to be $`\mathit{ϵ}_𝐤^\pm R(\widehat{𝐤})\mathit{ϵ}_p^\pm `$. That means we have to analyze the states

$$|\mathrm{\Psi }_\pm =𝑑\mu (𝐤)f(𝐤)|\mathit{ϵ}_𝐤^\pm ,𝐤,$$
(37)
where $`f(𝐤)`$ is given above.

We expand $`R(\widehat{𝐤})`$ up to second order in $`\theta `$. Reduced density matrices are calculated by techniques similar to those for massive particles , using rotational symmetry around the $`z`$-axis and normalization requirements. The leading order in $`\mathrm{\Omega }\mathrm{\Delta }_r/k_A`$ gives

$$\rho _+=\frac{1}{2}(1\frac{1}{2}\mathrm{\Omega }^2)\left(\begin{array}{ccc}1& i& 0\\ i& 1& 0\\ 0& 0& 0\end{array}\right)+\frac{1}{2}\mathrm{\Omega }^2\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right),$$
(38)
and $`\rho _{}=\rho _+^{}`$. At the same level of precision,

$$P_E(\rho _+,\rho _{})=\mathrm{\Delta }_r^2/4k_A^2.$$
(39)
It is interesting to note that an optimal strategy for distinguishing between these two states is a polarization measurement in the $`xy`$-plane. Then the effective $`2\times 2`$ density matrices are perfectly distinguishable, but there is a probability $`\mathrm{\Omega }^2/2`$ that no photon be detected at all. The above result was valid due to the special form of the states that we chose. Potential errors in the upper $`2\times 2`$ blocks were averaged out in the integration over $`\varphi `$.

We now calculate Bob’s reduced density matrix. We again assume that Bob moves along the $`z`$-axis with a velocity $`v`$. Recall that reduced density matrices have no transformation law , (only the complete density matrix has one) except in the limiting case of sharp momenta. The only way to calculate Bob’s reduced density matrix is to transform the complete state, and only then take a partial trace. According to Eqs. (31) and (37), reduced density matrices in both frames are given by the expression

$$(\rho _\pm )_{mn}=𝑑\mu (𝐤)|f(𝐤)|^2R(\widehat{𝐤})\mathit{ϵ}_p^\pm |𝐛_m(𝐤)𝐛_n(𝐤)|R(\widehat{𝐤})\mathit{ϵ}_p^\pm .$$
(40)
Note that pure boosts preserve the orientation of the coordinate axes in 3-space, and therefore do not affect the indices of $`\rho _{mn}`$. The measure $`\mu (𝐤)`$ is Lorentz-invariant and $`\mathit{ϵ}_p^\pm `$ are constant by definition. Since $`f`$ is a scalar function, it transforms as $`f^{}(𝐤)=f(𝐤_{\mathrm{\Lambda }^1})`$, where primes indicate Bob’s frame, as in Eqs. (69). This is the only frame dependent expression in (40), since the phase factors $`e^{\pm i\xi }`$ cancel out for the helicity eigenstates. Therefore, there are two equivalent methods to calculate Bob’s polarization density matrix. One is to change the argument of $`f(𝐤)`$ to $`𝐤_{\mathrm{\Lambda }^1}`$, and another is to change the argument of the rotation matrix $`R(\widehat{𝐤})`$ to $`\widehat{𝐤}_\mathrm{\Lambda }`$. Using the second option, we obtain

$$(\rho _\pm ^{})_{mn}=𝑑\mu (𝐤)|f(𝐤)|^2R(\widehat{𝐤}_\mathrm{\Lambda })\mathit{ϵ}_p^\pm |𝐛_m(𝐤)𝐛_n(𝐤)|R(\widehat{𝐤}_\mathrm{\Lambda })\mathit{ϵ}_p^\pm .$$
(41)
A boost along the $`z`$-axis preserves $`k_r`$ and $`\varphi `$. On the other hand, from Eq. (7) it follows that

$$k_z^{}k_A\sqrt{\frac{1v}{1+v}}.$$
(42)
Thus, at leading order in $`\theta `$ we have

$$\theta ^{}\sqrt{\frac{1+v}{1v}}k_r/k_A\sqrt{\frac{1+v}{1v}}\theta ,$$
(43)
which is substituted into $`R(𝐤_\mathrm{\Lambda })`$. Since everything else in the integral remains the same, the effect of relative motion is given by a substitution

$$\mathrm{\Omega }\sqrt{\frac{1+v}{1v}}\mathrm{\Omega }.$$
(44)
It follows that

$$P_E^{}=\frac{1+v}{1v}P_E,$$
(45)
which may be either larger or smaller than $`P_E`$. As expected, we obtain for one-photon states the same Doppler effect as in the preceding classical calculations.

Although reduced polarization density matrices have no general transformation rule, the above results, as well as the analysis of massive particles , show that such rules can be derived for particular classes of experimental procedures. We can then ask how these effective transformation rules $`\rho ^{}=T[\rho ]`$ fit into the framework of general state transformations. A general state transformation $`T`$ is usually required to be completely positive (CP), namely ,

$$T[\rho ]=\underset{i}{}M_i\rho M_i^{},$$
(46)
where the $`M_i`$ are bounded operators. It can be proved that distinguishability, as expressed by natural measures like $`P_E`$, cannot be improved by any CP transformation . It is also known that the CP requirement may fail if there is a prior entanglement of $`\rho `$ with another system . Since from  and Eq. (45) it follows that in our case distinguishability can be improved, we conclude that these transformations are not completely positive. The reason is that the Lorentz transformation acts not only on the “interesting” polarization variables, but also on the “hidden” momentum variables that we elected to ignore and to trace out.

This technicality has one important consequence. In quantum information theory quantum channels are described as completely positive maps  that act on qubit states. Qubits themselves are realized as particles’ discrete degrees of freedom. If relativistic motion is important, then not only does the vacuum behave as a noisy quantum channel, but the very representation of a channel by a CP map fails.

Acknowledgments

Work by AP was supported by the Gerard Swope Fund and the Fund for Encouragement of Research. DRT was supported by a grant from the Technion Graduate School.


format: text

text:
# Variational analysis of the deconfinement phase transition

## I Introduction

Attempts to understand the nature of the deconfining phase transition in QCD date back almost 25 years. Since the pioneering work of Polyakov polyakov , much effort has been made to study the basic physics as well as the quantitative characteristics of the transition. The high temperature phase of QCD is widely believed to resemble an almost free plasma of quarks and gluons. At asymptotically high temperatures this is confirmed by explicit perturbative calculations of the free energy perturb . Perturbation theory in its simplest form, however, is valid only at unrealistically high temperatures. In recent years a different promising avenue has been explored. This incorporates analytical resummation of the effects of the gluon screening mass into the 3D effective Lagrangian, which is then solved numerically by 3D lattice gauge theory methods nieto ; mikko . The results of this approach seem to be in good agreement hart  with direct 4D lattice gauge theory calculations karsch  all the way down to $`2𝖳_𝖼`$. The numerical results indicate that, interestingly enough, although asymptotically the free energy does approach that of the free partonic plasma, the deviations from the Stefan-Boltzmann law even at temperatures of order $`10𝖳_𝖼`$ are quite sizable, of order of $`15\%`$. One can interpret this as an indication that the interesting physics of the transition region remains important even at these high temperatures.

Although we are quite advanced in the understanding of the high temperature phase, the transition region itself is understood very poorly. This region of temperatures $`𝖳_𝖼<𝖳<2𝖳_𝖼`$, is of course the most interesting one, since it is in this region that the transition between “hadronic” and “partonic” degrees of freedom occurs. In fact, it is quite possible that the dynamics of the transition region is dominated by the $`Z_N`$ magnetic vortices thooft , which are responsible for the abrupt increase in entropy kaks . The study of the transition region is a complicated and inherently nonperturbative problem which has not been tackled so far by analytical methods.

The purpose of this paper is to explore an analytical investigation of the phase transition region. Since the nature of the problem is nonperturbative, a nonperturbative method is called for. Unfortunately there are not many of those on the market. It does however seem to us that the variational method developed several years ago in kk  is quite suitable. We are encouraged by the fact that the variational approximation of kk  manages to reproduce all known infrared features in lower dimensional confining theories kk1 ; kk2 ; ks ; brogakoko1 , and also gives a stable picture of instanton dynamics in QCD brogakoko . Even though, as any variational method, it has an element of guesswork about it, it is a first principle calculation. For such a difficult physics problem as the one at hand, clearly even an exploratory investigation which provides a qualitative picture should be very valuable and is well worth attempting.

In this paper, therefore, we will study the deconfining phase transition in a pure $`SU(N)`$ Yang-Mills theory using the variational approach of kk  suitably extended to finite temperature. We will minimize the relevant thermodynamic potential at finite temperature, i.e. the Helmholtz free energy, on a set of gauge invariant density matrices.

The results of our calculation are very interesting. We find that the phase transition in gluodynamics is closely related to the symmetry breaking phase transition in the effective $`\sigma `$-model for the Polyakov loop variable $`U`$. This $`\sigma `$-model arises explicitly in our calculation, and its parameters depend on the variational parameters of the trial density matrix. At low temperature the effective $`\sigma `$-model is in the disordered phase, and $`U=0`$. The entropy of the Yang-Mills thermal density matrix vanishes in the low temperature phase. At the critical temperature of order $`𝖳_𝖼450\mathrm{M}\mathrm{e}\mathrm{v}`$, the best variational density matrix changes abruptly, and above $`𝖳_𝖼`$ corresponds to the ordered phase of the $`\sigma `$-model. Thus at $`𝖳>𝖳_𝖼`$ we find $`U0`$. The transition itself turns out to be strongly first order for large $`N`$. This is consistent with recent lattice gauge theory calculations mike . These results are not surprising in themselves. They are expected from the general relation between the average value of the Polyakov loop and the (de)confining properties of the density matrix. Moreover, the order of the transition is predicted by the Svetitsky-Yaffe conjecture ben . The vanishing of the entropy in the confining phase is also very natural, since the glueballs are heavy and their contribution to the entropy should be suppressed by the Boltzmann factor $`\mathrm{exp}[M_g/T]`$. Nevertheless, we are not aware of any first principle analytic calculation in QCD which provides an explicit realization of these general arguments.

This paper is structured as follows. In sec. 2 we recap the variational ansatz of kk  for energy minimization at zero temperature and recall the methods of evaluating the expectation values in the variational wave function as well as the results of the minimization procedure. In sec. 3 we generalize the variational ansatz so that it incorporates not only pure states but also general density matrices. We introduce the two variational parameters with respect to which the free energy is to be minimized. In sec. 4 we recast the problem of the minimization of free energy in terms of an effective nonlinear $`\sigma `$-model. In sec. 5 we calculate the Helmholtz free energy of the trial density matrix and perform the minimization. We conclude in sec. 6 by discussing our results.

## II The Variational Calculation at Zero Temperature

The variational method proposed in kk  is based on minimization of the expectation value of the Hamiltonian on a set of gauge invariant states. The set of states is chosen so that their functional form is a relatively simple generalization of the vacuum of the free theory. Ideally, one would like to choose a set of variational states that is both simple enough to perform explicit calculations and rich enough to allow for variations of its parameters to span interesting physics. One starts then with a Gaussian state

$`\mathrm{\Psi }_0[A_i^a]=\mathrm{exp}\left\{{\displaystyle \frac{1}{2}}{\displaystyle d^3xd^3yA_i^a(x)(G^1)_{ij}^{ab}(x,y)A_j^b(y)}\right\},`$ (1)
where the set of functions $`G_{ab}^{ij}(x)`$ are variational parameters. Expectation values of various operators are easily calculable in this set of states. One also expects that the freedom allowed by variation of $`G`$ is wide enough to probe the nonperturbative physics of confinement. The problem is that these states are not gauge invariant and, therefore, as such do not belong to the physical Hilbert space of gluodynamics. To remedy this problem one projects the states onto the gauge invariant subspace by gauge transforming them and integrating over the whole gauge group:

$$\mathrm{\Psi }[A_i^a]=DU(x)\mathrm{exp}\left\{\frac{1}{2}d^3xd^3yA_i^{Ua}(x)G_{ij}^{1ab}(xy)A_j^{Ub}(y)\right\}$$
(2)
with $`A_i^{Ua}`$ defined as

$$A_i^{Ua}(x)=S^{ab}(x)A_i^b(x)+\lambda _i^a(x),$$
(3)
and

$`S^{ab}(x)={\displaystyle \frac{1}{2}}tr\left(\tau ^aU^{}\tau ^bU\right);\lambda _i^a(x)={\displaystyle \frac{i}{g}}tr\left(\tau ^aU^{}_iU\right)`$ (4)
with $`\tau `$ — traceless $`N\times N`$ hermitian matrices — normalized by

$$\mathrm{tr}(\tau ^a\tau ^b)=2\delta ^{ab},$$
(5)
and obeying the completeness condition for $`SU(N)`$

$$\tau _{ij}^a\tau _{kl}^a=2\left(\delta _{il}\delta _{jk}\frac{1}{N}\delta _{ij}\delta _{kl}\right).$$
(6)
The integration in eq. (2) is performed over the space of special unitary matrices with the $`SU(N)`$ group invariant measure. This integration projects the original Gaussian state onto a colour singlet. Due to the projection operation, the calculation of expectation values in this state is much more involved than in the case of a simple gaussian. Thus it is a very difficult task to perform the full functional minimization with respect to $`G_{ab}^{ij}(x)`$. This can be done and was indeed performed in simpler calculations in 2+1 dimensional models kk1 ; kk2 . In these models the approach was also extended to provide for the calculation of the confining potential between the external charges ks . It has also been applied to 1+1 dimensional models with fermions brogakoko1 . In these simple models the variational approximation has performed very well, reproducing all known results.

In nonabelian theories, unfortunately, full functional minimization is beyond our calculational abilities. Therefore to make headway, the width of the Gaussian $`G`$ was restricted in kk  to a one parameter form as follows

$$G_{ij}^{ab}(xy)=\delta ^{ab}\delta _{ij}G(xy)$$
(7)
with the Fourier transform of $`G`$ taken to be

$`G^1(k)=\{\begin{array}{cc}\sqrt{k^2}\hfill & \text{ if }k^2>M^2\hfill \\ M\hfill & \text{ if }k^2<M^2\hfill \end{array}.`$ (10)
The dimensional parameter $`M`$ is thus the only variational parameter in the calculation.

The expectation value of any gauge invariant operator $`O`$ in the variational state eq. (2) is given by the following integral

$`O`$ $`=`$ $`{\displaystyle \frac{1}{Z}}{\displaystyle DUO_A},`$ (11)
$`O_A`$ $`=`$ $`{\displaystyle DAe^{\frac{1}{2}{\scriptscriptstyle 𝑑x𝑑yA_i^{Ua}(x)G^1(xy)A_i^{Ua}(y)}}Oe^{\frac{1}{2}{\scriptscriptstyle 𝑑x^{}𝑑y^{}A_j^b(x^{})G^1(x^{}y^{})A_j^b(y^{})}}}.`$
Note that only one group integral is present. This is because for a gauge invariant operator $`O`$ only one of the states has to be gauge projected.

Since the gauge transform of a vector potential is a linear function of $`A`$ eq. (3), for fixed $`U(x)`$ this is a Gaussian integral, and can therefore be performed explicitly for any reasonable operator $`O`$. The nontrivial part of the calculation is the path integral over the group variable $`U(x)`$. Consider first the normalization factor $`Z`$. After integrating over the vector potential we obtain:

$$Z=DU\mathrm{exp}\{𝒮[U]\}$$
(12)
with the action

$$𝒮[U]=\frac{1}{2}Tr\mathrm{ln}+\frac{1}{2}\lambda [G+SGS^T]^1\lambda .$$
(13)
Here

$$S_{ij}^{ab}(x,y)=S^{ab}(x)\delta _{ij}\delta (xy),_{ij}^{ab}(x,y)=[S^{Tac}(x)S^{cb}(y)+\delta ^{ab}]G^1(xy)\delta _{ij},$$
(14)
and we use matrix product notation, which implies summation over the colour and Lorentz indices, as well as integration over the spatial coordinates. In the rest of this paper $`\mathrm{Tr}`$ stands for summation over all indices, including the space coordinates, while $`\mathrm{tr}`$ stands for summation over the colour indices alone. The path integral eq. (12) defines a partition function of a nonlinear $`\sigma `$-model with the target space $`SU(N)/Z_N`$ in three dimensional Euclidean space. The fact that the target space is $`SU(N)/Z_N`$ rather than $`SU(N)`$, follows from the observation that the action eq. (13) is invariant under local transformations belonging to the centre of $`SU(N)`$. This can be trivially traced back to invariance of $`A_i^a`$ under gauge transformations that belong to the centre of the gauge group.

The matrix $`U(x)`$ has a well defined gauge invariant meaning. It arises in eq. (II) as a relative gauge transformation between the bra and the ket gaussian wave functionals. As such it plays exactly the same role as Polyakov’s loop $`P`$ at finite temperature. We will therefore refer to it as the Polyakov loop in the following note1 . The contribution of a given $`U(x)`$ to the partition function eq. (12) and to other expectation values corresponds to the contribution due to the off diagonal matrix element between the initial Gaussian and the Gaussian gauge rotated by $`U(x)`$. Therefore, if matrices $`U(x)`$ which are far from unity give a significant contribution to the partition function, it means that the off diagonal contributions are large and gauge projection is indeed mandatory to get physically sensible results.

The action of this $`\sigma `$-model is rather complicated. It is a nonlocal and a nonpolynomial functional of $`U(x)`$. The following set of approximations was employed in kk  to estimate this path integral. First, since the bare coupling constant of the Yang Mills theory is small and it also enters as a coupling constant of the nonlinear $`\sigma `$-model, the high momentum modes of the field $`U`$ are integrated out perturbatively. The result of the perturbative one loop integration of the modes with spatial momenta $`k>M`$ is a much simpler $`\sigma `$-model which involves low momentum modes only. The action of this $`\sigma `$-model is

$$𝒮_L[U]=\frac{M}{2g^2(M)}trd^3x_iU^{}(x)_iU(x),$$
(15)
where $`M`$ also serves as the ultraviolet cutoff. To be more precise, due to the $`Z_N`$ local symmetry of the original theory eq. (13), the action for the low momentum modes is slightly different. The derivatives should be understood as $`Z_N`$ covariant derivatives. The most convenient way to write this action, would be to understand $`U(x)`$ as belonging to $`U(N)`$ rather than $`SU(N)`$, and to introduce a $`U(1)`$ gauge field by

$$𝒮_L=\frac{1}{2}\frac{M}{g^2}trd^3x(_iiA_i)U^{}(x)(_i+iA_i)U(x).$$
(16)
This defines a $`\sigma `$-model on the target space $`U(N)/U(1)`$, which is isomorphic to $`SU(N)/Z_N`$. However, since the difference between the two actions is not important at large number of colours, we will disregard this subtlety in the following and use eq. (15).

As a result of the running of the QCD coupling constant, the effective coupling $`g^2(M)`$ is given by the one loop QCD formula brko . Clearly for the perturbative integration to be consistent we should find that the energy is eventually minimized for high enough value of the scale $`M`$, so that $`\alpha _s(M)1`$. The actual value found in kk  is $`\alpha _s(M).25`$, which is indeed reasonably small. To solve the resulting low energy $`\sigma `$-model the mean field approximation and perturbation theory are used. We will describe this procedure here in some detail, since we will use it in later sections to analyze the finite temperature case.

The Hamiltonian of the pure Yang-Mills theory is

$$𝖧=d^3x\left[\frac{1}{2}E_i^{a2}+\frac{1}{2}B_i^{a2}\right],$$
(17)
where

$`E_i^a(x)`$ $`=`$ $`i{\displaystyle \frac{\delta }{\delta A_i^a(x)}},`$
$`B_i^a(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}ϵ_{ijk}\{_jA_k^a(x)_kA_j^a(x)+gf^{abc}A_j^b(x)A_k^c(x)\}.`$ (18)
To calculate the expectation value of $`𝖧`$ in the variational state we first perform the integral over the gauge fields $`A`$ at fixed value of the Polyakov loop $`U`$. The result of this integration (for details see kk ) for the energy density is

$`{\displaystyle \frac{2𝖧}{V}}`$ $`=`$ $`{\displaystyle \frac{3(N^21)}{2}}G^1(x,x)+(N^21)_i^x_i^yG(xy)|_{x=y}`$ (19)
$``$ $`{\displaystyle \frac{1}{4V}}{\displaystyle d^3xd^3y\lambda _{iH}^a(x)G^2(xy)\lambda _{iH}^a(y)_U}+{\displaystyle \frac{1}{4}}(ϵ_{ijk}_j\lambda _{kH}^a)^2_U`$
$``$ $`{\displaystyle \frac{M^2}{4V}}{\displaystyle d^3x\lambda _{iL}^a(x)\lambda _{iL}^a(x)_U},`$
where the averaging over the $`U`$ field should be performed with the $`\sigma `$-model action, eq. (13). As the high momentum modes are only considered to one loop order, the $`\sigma `$-model action reduces to

$$𝒮=\frac{1}{4}𝑑x𝑑y\lambda _{iH}^a(x)G^1(xy)\lambda _{iH}^a(y)+\frac{M}{2g^2(M)}trd^3x_iU_L^{}(x)_iU_L(x).$$
(20)
The subscripts $`H`$ and $`L`$ indicate that the corresponding expressions contain only high and low momentum components of the fields respectively. Integration over the high momentum modes gives

$$\frac{2𝖧}{V}=\frac{N^21}{10\pi ^2}M^4\frac{M^2}{4V}d^3x\lambda _{iL}^a(x)\lambda _{iL}^a(x)_U,$$
(21)
where now only the contribution of the low momentum modes remains to be evaluated.

This last step of the calculation is the most difficult one. The nonlinear $`\sigma `$-model eq. (15) possesses a global $`SU_L(N)SU_R(N)`$ symmetry. Thus it has a symmetry restoring phase transition as a function of $`M`$. The coupling constant $`g^2(M)`$ enters this action in the same way as temperature enters in the classical statistical problem of the $`SU(N)`$ spin system. High values of $`M`$ correspond to low values of $`g^2(M)`$ and thus to the low temperature regime in the $`\sigma `$-model. In this regime the global $`SU_L(N)SU_R(N)`$ symmetry is spontaneously broken down to the vector subgroup $`SU_V(N)`$. The model itself is weakly interacting and thus one can use perturbation theory. However, it is well known that the model undergoes a phase transition from the ordered low temperature phase to the disordered phase at some critical temperature $`g_c^2`$. The disordered phase is clearly nonperturbative. To treat it quantitatively the mean field approximation was used in kk . The partition function of the $`\sigma `$-model is rewritten by introducing a (hermitian matrix) auxiliary field $`\alpha `$ which imposes a unitarity constraint on $`U(x)`$

$`Z`$ $`=`$ $`{\displaystyle DUD\alpha \mathrm{exp}\left(𝒮[U,\alpha ]\right)},`$ (22)
$`𝒮[U,\alpha ]`$ $`=`$ $`{\displaystyle \frac{M}{2g^2(M)}}tr{\displaystyle d^3x\left[_iU^{}(x)_iU(x)+\alpha \left(U^{}U1\right)\right]}.`$
In this step we have disregarded the fact that $`U`$ is an $`SU(N)`$ matrix rather than $`U(N)`$. Thus, strictly speaking, our calculation applies only at large $`N`$. The mean field equations for this action are

$$U^{}U=1,$$
(23)
$$\alpha U=0.$$
(24)
From eq. (24) it follows that either $`\alpha =0`$, $`U0`$ (the ordered, broken symmetry phase with massless Goldstone bosons), or $`\alpha 0`$, $`U=0`$ (the disordered, unbroken phase with massive excitations).

In the disordered phase the expectation value of $`\alpha `$ is proportional to the unit matrix

$$\alpha _{\alpha \beta }=\alpha ^21_{\alpha \beta }.$$
(25)
Eq. (23) then becomes

$$2N^2\frac{g^2(M)}{M}_0^M\frac{d^3k}{(2\pi )^3}\frac{1}{k^2+\alpha ^2}=\frac{N^2g^2(M)}{\pi ^2}\left(1\frac{\alpha }{M}\mathrm{arctan}\frac{M}{\alpha }\right)=N.$$
(26)
The gap equation, eq. (26), has a solution only for couplings $`g^2(M)`$ larger than the critical coupling $`g_c^2`$, which is determined by the condition that $`\alpha =0`$

$$\alpha _s^c=\frac{g_c^2}{4\pi }=\frac{\pi }{4}\frac{1}{N}.$$
(27)
The low momentum mode contribution to the ground state energy is

$$N^2M_0^M\frac{d^3k}{(2\pi )^3}\frac{k^2}{k^2+\alpha ^2}=\frac{N^2}{2\pi ^2}M\left[\frac{1}{3}M^3\alpha ^2M+\alpha ^3\mathrm{arctan}\frac{M}{\alpha }\right].$$
(28)
The final mean field expression for the ground state energy density is

$$𝖤=\frac{N^2}{4\pi ^2}M^4\left[\frac{2}{15}+\frac{\alpha ^2}{M^2}\frac{\alpha _s^c}{\alpha _s(M)}\right],$$
(29)
where $`\alpha _s(M)`$ is the QCD coupling at the scale $`M`$, $`\alpha _s^c`$ is given by eq. (27), and $`\alpha `$ is determined by

$$\frac{\alpha }{M}\mathrm{arctan}\frac{M}{\alpha }=\frac{\alpha _s(M)\alpha _s^c}{\alpha _s(M)}.$$
(30)
This function has the minimum at the critical point $`\alpha _s(M)=\alpha _s^c`$. With the one loop Yang-Mills $`\beta `$ function and $`\mathrm{\Lambda }_{QCD}=150\mathrm{M}\mathrm{e}\mathrm{v}`$, we find for $`N=3`$

$$M_c=\mathrm{\Lambda }_{QCD}e^{\frac{24}{11}}=8.86\mathrm{\Lambda }_{QCD}=1.33\mathrm{Gev}.$$
(31)
Eq. (29) is valid only in the disordered phase. Here we expect the mean field approximation to be reasonably good, since the spectrum of the $`\sigma `$-model is massive and fluctuations are not expected to be overwhelmingly large.

In the ordered phase, on the other hand, the situation is different. The spectrum contains $`O(N^2)`$ “Goldstone bosons”. In the mean field approximation in addition it also has $`O(N^2)`$ massive particles — the analogs of $`\delta `$ particles in QCD. In the mean field the $`\delta `$’s are all stable. In reality of course, since there is nothing to prevent them from decaying into the Goldstone bosons, they are not at all stable. In fact, the further into the ordered phase one gets, the less stable they are due to the increase of the phase space available for decay. Thus we expect the mean field to be unreliable in the ordered phase. Asymptotically at very large $`M`$ the mass of $`\delta `$ becomes much higher than the ultraviolet cutoff of the $`\sigma `$-model. For these values of $`M`$ the massive modes decouple even within the mean field approximation, and one recovers the perturbative result. This decoupling, however, occurs extremely slowly. We have checked that the mean field formulae tend to the perturbative ones only when the mass of the would be stable massive states is about 1000 times the ultraviolet cutoff. This is, of course, unphysical, confirming our suspicion that the mean field approximation cannot be used on the ordered side of the transition.

The simplest option to calculate the contribution of the low momentum modes to the energy in the ordered phase is to use perturbation theory. The raison d’etre for this is the following. Perturbation theory is certainly appropriate for large enough values of $`M`$, where the expectation value of the $`U`$ field is of order unity. On the other hand, we know from the numerical studies of kogut  that the phase transition in our $`\sigma `$-model is strongly first order. More specifically, according to kogut  the transition occurs when the expectation value of $`U`$ is greater than $`.5`$. We thus expect perturbation theory to be qualitatively reliable all the way down to the transition point. Since we have neglected previously many perturbative contributions of order $`g^2`$, in order to be consistent we have to limit ourselves to the leading order contribution only. Calculating with this precision the energy of the trial states for $`MM_c`$ we get

$$\frac{𝖤(M)}{V}=\frac{N^21}{120\pi ^2}M^4.$$
(32)
Thus for $`MM_c`$ the energy is a monotonically increasing function of $`M`$.

Eq. (31) is thus the result of the minimization of the energy over the whole range of the variational parameter $`M`$.

## III The variational Ansatz for the density matrix

In order to extend the variational analysis to finite temperature we have to generalize our ansatz so that it includes mixed states. In scalar theories the Gaussian approximation has a long history of applications at finite temperature eboli ; moshe . We generalize our ansatz along the same lines.

We start by considering the density matrices which in the field basis have Gaussian matrix elements

$$\varrho [A,A^{}]=\mathrm{exp}\left\{\frac{1}{2}_{x,y}A_i^a(x)G_{ij}^{1ab}(x,y)A_j^b(y)+A_i^{}_{}{}^{}a(x)G_{ij}^{1ab}(x,y)A_j^{}_{}{}^{}b(y)2A_i^a(x)H_{ij}^{ab}(x,y)A_j^{}_{}{}^{}b(y)\right\}.$$
(33)
As before, we take the variational functions diagonal in both colour and Lorentz indices, and translationally invariant

$`G_{ij}^{1ab}(x,y)`$ $`=`$ $`\delta ^{ab}\delta _{ij}G^1(xy),`$ (34)
$`H_{ij}^{ab}(x,y)`$ $`=`$ $`\delta ^{ab}\delta _{ij}H(xy).`$
Then

$$\stackrel{~}{\varrho }[A,A^{}]=\mathrm{exp}\{\frac{1}{2}_{x,y}(AG^1A+A^{}G^1A^{}2AHA^{}\left)\right\}.$$
(35)
Note that for $`H=0`$ this density matrix represents a pure state, since it can be written in the form

$$\stackrel{~}{\varrho }=|\mathrm{\Psi }[A]><\mathrm{\Psi }[A]|$$
(36)
with $`\mathrm{\Psi }[A]`$ a gaussian wave function. At nonzero $`H`$ the density matrix is, however, mixed. The magnitude of $`H`$, therefore, determines the entropy of this trial density matrix.

We now make an additional simplification in our ansatz. First, we restrict the functions $`G^1(x)`$ to the same functional form as at zero temperature, eq. (10). Further, we will take $`H(k)`$ to be small and nonvanishing only at low momenta

$`H(k)=\{\begin{array}{cc}0\hfill & \text{ if }k^2>M^2\hfill \\ HM\hfill & \text{ if }k^2<M^2\hfill \end{array}.`$ (39)
The logic behind this choice of ansatz is the following. At finite temperature we expect $`H(k)`$ to be roughly proportional to the Bolzmann factor $`\mathrm{exp}\{𝖤(k)\beta \}`$. In our ansatz, the role of one particle energy is played by the variational function $`G^1(k)`$. We will be interested only in temperatures close to the phase transition, and those we anticipate to be small, $`𝖳_𝖼M`$. For those temperatures one particle modes with momenta $`kM`$ are not populated, and we thus put $`H(k)=0`$. For $`kM`$ the Bolzmann factor is nonvanishing, but small. Further, it depends only very weakly on the value of the momentum. We will have, of course, to verify a posteriori that our assumptions about the smallness of $`𝖳_𝖼`$ and $`H`$ are justified. As we will see later, this turns out indeed to be the case with reasonable accuracy.

As before, we explicitly impose gauge invariance by projecting $`\varrho `$ onto a gauge invariant sector

$$\varrho [A,A^{}]=𝒟U^{}𝒟U^{\prime \prime }\mathrm{exp}\left\{\frac{1}{2}_{x,y}A^U^{}G^1A^U^{}+A_{}^{}{}_{}{}^{U^{\prime \prime }}G^1A_{}^{}{}_{}{}^{U^{\prime \prime }}2A^U^{}HA_{}^{}{}_{}{}^{U^{\prime \prime }}\right\}.$$
(40)
One of the group integrations in eq. (40) is redundant, since we will only calculate the quantities of the form $`\mathrm{Tr}\varrho O`$, with $`O`$ being gauge invariant. We thus have the following ansatz for the density matrix

$$\varrho [A,A^{}]=𝒟U\mathrm{exp}\left\{\frac{1}{2}_xAG^1A+A_{}^{}{}_{}{}^{U}G^1A_{}^{}{}_{}{}^{U}2AHA_{}^{}{}_{}{}^{U}\right\}.$$
(41)
This expression is not explicitly normalized to unity. Nevertheless, we find it convenient to refer to it as density matrix while explicitly inserting a normalization factor whenever necessary. Thus the average of a gauge invariant operator $`𝒪`$ is given by

$`𝒪_{A,U}`$ $`=`$ $`Z^1\mathrm{Tr}(\varrho 𝒪)`$ (42)
$`=`$ $`Z^1{\displaystyle 𝒟U𝒟A𝒪(A,A^{})\mathrm{exp}\left\{\frac{1}{2}_xAG^1A+A_{}^{}{}_{}{}^{U}G^1A_{}^{}{}_{}{}^{U}2AHA_{}^{}{}_{}{}^{U}\right\}}|_{A^{}=A},`$
where $`Z`$ is the normalization of the trial density matrix $`\varrho `$, i.e.

$`Z`$ $`=`$ $`\mathrm{Tr}\varrho ={\displaystyle 𝒟U𝒟A\mathrm{exp}\left\{\frac{1}{2}_xAG^1A+A^UG^1A^U2AHA^U\right\}}`$ (43)
$`=`$ $`{\displaystyle 𝒟U𝒟\stackrel{~}{A}\mathrm{exp}\left\{\frac{1}{2}_x\stackrel{~}{A}\mathrm{\Delta }\stackrel{~}{A}+\lambda \left(G^1\omega \mathrm{\Delta }^1\omega ^T\right)\lambda \right\}}`$
with

$`\stackrel{~}{A}`$ $`=`$ $`A+\lambda \omega \mathrm{\Delta }^1,`$ (44)
$`\mathrm{\Delta }`$ $`=`$ $`2G^1\left(1{\displaystyle \frac{HG}{2}}(S+S^T)\right),`$ (45)
$`\omega `$ $`=`$ $`(G^1SH).`$ (46)
The $`\stackrel{~}{A}`$ integration can then be easily performed to yield

$$\mathrm{Tr}\varrho =𝒟U\mathrm{exp}\left\{\frac{1}{2}\lambda \left(G^1\omega \mathrm{\Delta }^1\omega ^T\right)\lambda \frac{3}{2}\mathrm{Tr}\mathrm{ln}\frac{\mathrm{\Delta }}{2}\right\}.$$
(47)
We now adopt the same strategy for treating the high momentum modes of $`U`$ as at $`𝖳=0`$. Namely, they are integrated perturbatively to one loop accuracy. The result is the effective $`\sigma `$-model for the matrices $`U`$ with momenta below $`M`$. The coupling constant $`g`$ of this $`\sigma `$-model gets renormalized as before according to the one loop Yang-Mills $`\beta `$-function, and thus has to be understood as $`g(M)`$. Additionally, due to independence of $`H`$ on momentum, for low momentum modes of $`U`$ the function $`H(xy)`$ is equivalent to $`H\delta ^3(xy)`$.

The final approximation has to do with the fact that $`H`$ is assumed to be small. For arbitrarily large $`H`$ the variational calculation is forbiddingly complicated even with all the above mentioned simplifications. This is because the gauge projection renders the calculation of entropy in the general case unfeasible. However, at small $`H`$ we only need to calculate the leading term in entropy. This calculation can indeed be done, and is described in the next section. Since we are only calculating the leading order contribution in $`H`$, we only have to consider corrections to the $`\sigma `$-model action of first order in $`H`$. With this in mind, the normalization factor becomes

$$\mathrm{Tr}\varrho =𝒟U\mathrm{exp}\left\{\frac{1}{2}\lambda \left(\frac{G^1}{2}+\frac{H}{4}(S+S^T)\right)\lambda +\frac{3}{4}HG\mathrm{tr}(S+S^T)\right\}.$$
(48)
## IV The effective $`\sigma `$-model

The normalization $`Z`$ can be interpreted as the generating functional for a theory defined by the action $`𝒮(U)`$

$$Z=\mathrm{Tr}\varrho =𝒟Ue^{𝒮(U)},$$
(49)
where

$$𝒮(U)=\frac{M}{4}\lambda \lambda +\frac{1}{8}\lambda H(S+S^T)\lambda \frac{1}{4\pi ^2}HM^2\mathrm{tr}S.$$
(50)
We simplify this expression using

$`\lambda \lambda `$ $`=`$ $`{\displaystyle \frac{2}{g^2}}\mathrm{tr}(UU^{}),`$ (51)
$`\lambda S^T\lambda `$ $`=`$ $`\lambda S\lambda ={\displaystyle \frac{1}{2g^2}}\mathrm{tr}\left[(U^{}UU^{}U)(UU^{}UU^{})\right],`$ (52)
$`\mathrm{tr}S`$ $`=`$ $`\mathrm{tr}S^T=\mathrm{tr}U^{}\mathrm{tr}U1.`$ (53)
Inserting these into the action we get

$$𝒮(U)=\frac{M}{2g^2}\mathrm{tr}(UU^{})\frac{H}{8g^2}\mathrm{tr}\left[(U^{}UU^{}U)(UU^{}UU^{})\right]\frac{1}{4\pi ^2}HM^2\mathrm{tr}U^{}\mathrm{tr}U,$$
(54)
where $`U`$ independent pieces have been dropped.

As noted before, this effective $`\sigma `$-model is nothing but the effective theory for the low momentum modes of the Polyakov loop variable. Its status and applicability region are different from the usual perturbative effective actions, see e.g. chris . The standard effective action is calculated in perturbation theory and is valid at high temperature. Our effective action eq. (54) depends on the variational parameters $`M`$ and $`H`$, and in a sense is a variational effective action. Also due to our restrictions to small values of $`H`$, a priori we do not expect it to be valid at high temperatures but, it rather, should represent correctly the physics in the phase transition region.

Another important difference is that our effective $`\sigma `$-model does not have the local gauge invariance $`U(x)V^{}(x)U(x)V(x)`$ which is usually associated with the effective action for the Polyakov loop. The reason for this is that our setup is different from that of the standard finite temperature calculation. The way this gauge invariance usually appears is the following. Consider the calculation of any gauge invariant observable

$$𝒪=𝒟U\mathrm{Tr}[\mathrm{exp}\{\beta 𝖧\}𝒪g(U)],$$
(55)
where $`g(U)`$ is the quantum mechanical operator of the gauge transformation represented by the matrix $`U`$. This expression for fixed $`U`$ can be compared to the same expression but with $`U`$ gauge transformed

$$\mathrm{Tr}[\mathrm{exp}\{\beta 𝖧\}𝒪g(V^{}UV)]=\mathrm{Tr}[\mathrm{exp}\{\beta 𝖧\}𝒪g(V^{})g(U)g(V)]=\mathrm{Tr}[\mathrm{exp}\{\beta 𝖧\}𝒪g(U)].$$
(56)
The last equality here follows from the fact that both $`𝒪`$ and $`\mathrm{exp}\{\beta 𝖧\}`$ are gauge invariant, and thus the operator $`g(V^{})`$ can be commuted all the way to the left. The only effect of the transformation is then to change the basis over which the trace is being taken, which obviously leaves the trace invariant.

Our setup is somewhat different. Expectation values are calculated as

$$𝒟U\mathrm{Tr}[\stackrel{~}{\varrho }g(U)𝒪]$$
(57)
with $`\stackrel{~}{\varrho }`$ defined in eq. (35). This expression is altogether gauge invariant, since the integral over $`U`$ correctly projects only the contribution of gauge singlet states. However the operator $`\stackrel{~}{\varrho }`$ is not itself explicitly gauge invariant. For that reason the gauge transformation operator $`g(V^{})`$ cannot be commuted through it, and thus

$$\mathrm{Tr}[\stackrel{~}{\varrho }𝒪g(V^{}UV)]\mathrm{Tr}[\stackrel{~}{\varrho }𝒪g(U)]$$
(58)
even for gauge invariant operators $`𝒪`$. This manifests itself as absence of local gauge invariance in the action of the effective $`\sigma `$-model, eq. (54).

Nevertheless, we stress again that the meaning of the $`SU(N)`$ valued field $`U`$ is precisely the same as that of the Polyakov loop.

## V The Calculation of the Free Energy

To find the best variational density matrix we have to minimize the free energy with respect to the variational parameters $`M`$ and $`H`$. The Helmholtz free energy $`𝖥`$ of the density matrix $`\varrho `$ is given by

$$𝖥=𝖧\mathrm{𝖳𝖲},$$
(59)
where $`𝖧`$ is the standard Yang-Mills Hamiltonian eq. (17), $`𝖲`$ is the entropy, and $`𝖳`$ is the temperature.

Thus

$$𝖥=\frac{1}{2}\left(\mathrm{Tr}(E^2\varrho )+\mathrm{Tr}(B^2\varrho )\right)+𝖳\mathrm{Tr}(\varrho \mathrm{ln}\varrho ).$$
(60)
First of all we need to perform the integration over the gauge fields, and reduce this expression to the average of a $`U`$-dependent operator in the effective $`\sigma `$-model. In fact, as we shall see soon, to leading order in $`H`$ the only nontrivial calculation we need to perform is that of the entropy.

We will calculate the entropy up to the first nontrivial order in $`H`$. As we now show, the leading term at small $`H`$ is $`O(H\mathrm{ln}H)`$.

Let us denote by $`\varrho _0`$ the density matrix of the pure state with $`H=0`$:

$$\varrho _0=|00|.$$
(61)
Here $`|0`$ does not denote necessarily the actual ground state, but rather a projected Gaussian state with arbitrary $`M`$. Now, since the matrix elements of the density matrix can be expanded in powers of $`H`$, to leading order we can write

$$\varrho =\varrho _0+\delta \varrho ,$$
(62)
where $`\delta \varrho `$ is $`O(H)`$.

Imagine that we have diagonalized $`\varrho `$. It will have one large eigenvalue $`\alpha _0=1O(H)`$, which corresponds to the eigenstate

$$|0^{}=|0+O(H).$$
(63)
All the rest of the eigenvalues $`\alpha _i`$ are at most $`O(H)`$. Then the entropy can be written as

$$𝖲=\mathrm{Tr}(\varrho \mathrm{ln}\varrho )=\alpha _0\mathrm{ln}\alpha _0\underset{i=1}{\overset{\mathrm{}}{}}\alpha _i\mathrm{ln}\alpha _i.$$
(64)
The second term is $`O(H\mathrm{ln}H)`$, and it is the coefficient of this term that we will now calculate. Neglecting $`O(H)`$ corrections, we can substitute $`\alpha _i=H/M`$ under the logarithm. Thus to leading logarithmic order

$$𝖲=\underset{i}{}\alpha _i\mathrm{ln}H/M.$$
(65)
So that all we have to calculate is $`_i\alpha _i`$. Let

$$|0^{}=|0+H|x.$$
(66)
Then

$$\varrho =\alpha _0|0^{}0^{}|+\underset{i=1}{\overset{\mathrm{}}{}}\alpha _i|x_ix_i|$$
(67)
with

$$x_i|0^{}=0.$$
(68)
Note that $`0|x0`$, but

$$0|x+x|0=0,$$
(69)
since $`|0^{}`$ has to be normalized at $`O(H)`$. Also

$$x_i|0+Hx_i|x=0.$$
(70)
Thus the overlap $`x_i|0`$ is $`O(H)`$, and we have

$$\varrho =\alpha _0|00|+H\left(|0x|+|x0|\right)+\alpha _i|x_ix_i|.$$
(71)
Multiplying this by $`\varrho _0`$ we get

$`\varrho _0\varrho `$ $`=`$ $`\alpha _0\varrho _0+H|0x|+H0|x|00|,`$ (72)
$`\varrho \varrho _0`$ $`=`$ $`\alpha _0\varrho _0+H|x0|+Hx|0|00|.`$
Thus,

$$\varrho _0\varrho +\varrho \varrho _0\varrho =\alpha _0\varrho _0\alpha _i|x_ix_i|.$$
(73)
Multiplying by $`\varrho _0`$ again, we get rid of $`|x_ix_i|`$ to $`O(H)`$

$$\alpha _0\varrho _0=\varrho _0\varrho +\varrho _0\varrho \varrho _0\varrho _0\varrho =\varrho _0\varrho \varrho _0.$$
(74)
Then,

$$\alpha _0=\mathrm{Tr}(\varrho _0\varrho ).$$
(75)
Since $`\mathrm{Tr}\varrho =1`$ we have

$$\underset{i}{}\alpha _i=1\alpha _0=\mathrm{Tr}(\varrho _0(1\varrho ))$$
(76)
which, inserted into eq. (64), gives

$$𝖲=(1\mathrm{Tr}(\varrho _0\varrho ))\mathrm{ln}H/M.$$
(77)
The derivation has been given for the normalized density matrices $`\varrho _0`$ and $`\varrho `$. In terms of our Gaussian matrices we should restore the normalization factors $`Z`$ and $`Z_0`$, so that finally we have

$$𝖲=\left(\frac{\mathrm{Tr}(\varrho _0\varrho )}{\mathrm{Tr}\varrho _0\mathrm{Tr}\varrho }1\right)\mathrm{ln}H/M.$$
(78)
It is easy to check that to $`O(H)`$

$$\mathrm{Tr}(\varrho _0\varrho )=(\mathrm{Tr}\varrho _0)^2.$$
(79)
Then

$$𝖲=\left(\frac{\mathrm{Tr}\varrho _0}{\mathrm{Tr}\varrho }1\right)\mathrm{ln}H/M.$$
(80)
From eq. (48), it is clear that

$$\mathrm{Tr}\varrho =\left[1+H\left(\frac{1}{4\pi ^2}M^2\mathrm{tr}S\frac{1}{4}\lambda S\lambda \right)\right]\mathrm{Tr}\varrho _0$$
(81)
Using eqs. (52,53) we finally get

$$𝖲=\left[<\frac{1}{8g^2}\mathrm{tr}(U^{}UU^{}U)(UU^{}UU^{})+\frac{1}{4\pi ^2}M^2(\mathrm{tr}U^{}\mathrm{tr}U1)>_U\right]H\mathrm{ln}H/M.$$
(82)
To leading order in $`H`$, the averaging over $`U`$ in this expression has to be performed with the $`\sigma `$-model action at vanishing $`H`$.

Expression eq. (82) has the following striking property. For $`M<M_c`$ it vanishes identically. The reason is very simple. The first term in eq. (82) is the product of the left handed $`SU(N)`$ current and the right handed $`SU(N)`$ current in the $`\sigma `$-model. Thus it transforms as an adjoint representation under each one of the $`SU(N)`$ factors of the $`SU_L(N)SU_R(N)`$ transformation. The same is also true for the second term in eq. (82). The $`\sigma `$-model action itself at $`H=0`$ obviously is invariant under the whole $`SU_L(N)SU_R(N)`$ group. Now, at $`M<M_c`$, the symmetry group is not spontaneously broken, and thus any operator which is not a scalar has a vanishing expectation value. It follows immediately that the entropy has an $`O(H\mathrm{ln}H/M)`$ contribution only for $`M>M_c`$, when the $`SU_L(N)SU_R(N)`$ group is spontaneously broken down to $`SU_V(N)`$.

This observation makes our task considerably simpler. Since for $`M<M_c`$ the entropy is zero, we do not have to consider at all the disordered phase of the effective $`\sigma `$-model. In this disordered phase the free energy coincides with energy, and thus the calculation is identical to the calculation at zero temperature. There is one slight subtlety here. We expect of course that since the energy alone must have a minimum on a pure state, the nonvanishing $`H`$ should give always a positive contribution to the energy. Thus in the absence of the entropy contribution, $`H`$ should vanish in the minimal energy state. Nevertheless since in kk  the energy at $`𝖳=0`$ was only minimized with respect to $`M`$, we have to check this property explicitly. We have indeed performed this check and found that for small $`H`$ the derivative of the expectation value of the energy with respect to $`H`$ is strictly positive.

Thus we only need to consider the effective $`\sigma `$-model in the ordered phase. As at $`𝖳=0`$ we perform the calculations in the ordered phase to leading order in $`\alpha _s`$. Since there are no $`O(H\mathrm{ln}H/M)`$ corrections to energy at this order, the result for the energy in the disordered phase is identical to the result at zero temperature, eq. (32). Thus our expression for the free energy is

$$𝖥=\frac{N^21}{120\pi ^2}M^4+𝖳\left(\frac{1}{8g^2}\mathrm{tr}(U^{}UU^{}U)(UU^{}UU^{})+\frac{1}{4\pi ^2}M^2(\mathrm{tr}U^{}\mathrm{tr}U1)_U\right)H\mathrm{ln}H/M.$$
(83)
We now average over $`U`$ in the leading order perturbation theory.

## VI The $`\sigma `$-model perturbation theory — minimization of the free energy and the Debye mass

We parameterize the $`U`$ matrices as

$$U=\mathrm{exp}\left\{\frac{i}{2}\varphi ^a\tau ^a\right\}.$$
(84)
Although we only need the leading order, it is instructive to check that the order $`g^2`$ term in the expansion is indeed small. To this order we have

$$U\left(1\frac{i}{2}\varphi ^a\tau ^a\frac{1}{8}\varphi ^a\varphi ^b\tau ^a\tau ^b+\frac{i}{48}\varphi ^a\varphi ^b\varphi ^c\tau ^a\tau ^b\tau ^c\right).$$
(85)
So that the $`\sigma `$-model action becomes

$`𝒮`$ $`=`$ $`{\displaystyle \frac{M}{2g^2}}\mathrm{tr}(UU^{})`$
$`=`$ $`{\displaystyle \frac{M}{4g^2}}\varphi \varphi +{\displaystyle \frac{M}{192g^2}}(\varphi ^a)(\varphi ^c)\varphi ^b\varphi ^d\mathrm{tr}\left[\tau ^a\tau ^b\tau ^c\tau ^d\tau ^a\tau ^c\tau ^b\tau ^d\right].`$
The propagator of the phase field $`\varphi `$ is thus

$$\varphi ^a\varphi ^b=\frac{2g^2}{Mk^2}\delta ^{ab}.$$
(87)
To get the idea of the quality of this perturbative expansion we can calculate for example $`𝒮`$. In this calculation one has to take into account the fact that the measure in the path integral over the phase $`\varphi ^a`$ is not the simple $`𝒟\varphi `$, but rather the group invariant $`U(N)`$ measure $`\mu `$ . To first order in $`g^2`$ it is

$$\mu =𝒟\varphi ^a\mathrm{exp}\left\{\frac{M^3N}{144\pi ^2}d^3x\varphi ^2(x)\right\}.$$
(88)
Taking this into account we find that $`𝒮`$ gets no correction of order $`g^2`$. We thus feel confident that the use of the perturbation theory in the ordered phase of the $`\sigma `$-model is an admissible approximation. In the following we will only keep leading order expressions.

Calculating to leading order the entropy eq. (82) and keeping only the $`O(N^2)`$ terms we find

$$𝖲=\frac{N^2}{6\pi ^2}M^2H\mathrm{ln}\frac{H}{M}.$$
(89)
Introducing the dimensionless quantity

$$h=\frac{H}{M}$$
(90)
we can write the expression for the free energy as

$`𝖥`$ $`=`$ $`𝖧𝖳𝖲`$ (91)
$`=`$ $`{\displaystyle \frac{N^2}{120\pi ^2}}M^4+𝖳{\displaystyle \frac{N^2}{6\pi ^2}}M^3h\mathrm{ln}h.`$
We now have to minimize this expression with respect to $`h`$ and $`M`$. It is convenient to first perform the minimization with respect to $`h`$ at fixed $`M`$. This obviously gives

$$\frac{𝖥}{h}=0h=\frac{1}{e}.$$
(92)
Thus as a function of $`M`$ only, the free energy becomes

$$𝖥=\frac{N^2}{120\pi ^2}M^4\frac{𝖳}{e}\frac{N^2}{6\pi ^2}M^3.$$
(93)
Now minimizing with respect to $`M`$ we find

$$\frac{𝖥}{M}=0M=\frac{15𝖳}{e}.$$
(94)
Thus for $`MM_c`$ the free energy of the best variational density matrix as a function of temperature is

$$𝖥_{MM_c}=\frac{N^2}{360\pi ^2}\left(\frac{15𝖳}{e}\right)^4.$$
(95)
We now have to compare this value with the free energy for $`MM_c`$. As we have discussed above, this is given by the expectation value of the Hamiltonian alone, and is minimized at $`M=M_c`$. Its value is

$$𝖥_{MM_c}=\frac{N^2}{30\pi ^2}M_c^4.$$
(96)
Comparing the two expressions we find

$$𝖳_𝖼=\frac{12^{\frac{1}{4}}e}{15}M_c.$$
(97)
Using the value of $`M_c`$ from eq. (31) we have

$$𝖳_𝖼=450\mathrm{M}\mathrm{e}\mathrm{v}.$$
(98)
For $`𝖳𝖳_𝖼`$ the free energy is minimized in the variational state with $`M=M_c`$. In our approximation this state is the same as at zero temperature. Its entropy vanishes, and the effective $`\sigma `$-model is in the disordered phase. The Polyakov loop vanishes, $`U=0`$ and according to the standard wisdom this is a confining state.

For $`𝖳𝖳_𝖼`$ the best variational state is very different. The entropy of this state is nonzero,

$$𝖲=\frac{N^2}{6\pi ^2e}\left(\frac{15𝖳}{e}\right)^3.$$
(99)
The Polyakov loop is nonzero $`U0`$ and thus the high temperature density matrix describes a deconfined phase.

Finally, we note that in the deconfined phase our best variational density matrix has a nonvanishing “electric screening” or “Debye” mass. The Debye mass is conveniently defined as the “mass” of the phase of the Polyakov loop. This mass is nonvanishing in our calculation for the following reason. As long as $`H=0`$, the effective $`\sigma `$-model action has a global $`SU_L(N)SU_R(N)`$ symmetry. Thus in the ordered phase of the $`\sigma `$ model the phases $`\varphi ^a`$ are massless. However, as discusses above, the terms of order $`H`$ in eq.(54) break this symmetry explicitly down to the diagonal $`SU_V(N)`$. As a result the would be “Goldstone” phases $`\varphi ^a`$ acquire mass. To calculate this mass it is convenient first to note that to $`O(g^2)`$

$$\mathrm{tr}(U^{}UU^{}U)(UU^{}UU^{})=4\mathrm{t}\mathrm{r}U^{}U\frac{1}{4}\varphi ^a\varphi ^c\varphi ^b\varphi ^d\mathrm{tr}\left(\tau ^a\tau ^b\tau ^c\tau ^d\tau ^a\tau ^c\tau ^b\tau ^d\right).$$
(100)
The contribution of the $`SU_L(N)SU_R(N)`$ term to the mass cancels against the contribution of the measure eq.(88). Using eqs.(54, 100) we then find to $`O(g^2)`$ and to leading order in $`H`$

$$M_D^2=\frac{4}{3\pi }\alpha _s(M)NMH.$$
(101)
As a function of temperature we have

$$M_D^2=\alpha _s\left(\frac{15}{e}T\right)N\frac{300}{\pi e^3}𝖳^2.$$
(102)
## VII Discussion

Let us summarize the results of our variational calculation. We find the phase transition at a temperature of about $`𝖳_𝖼450\mathrm{M}\mathrm{e}\mathrm{v}`$. The transition is strongly first order at large $`N`$. The latent heat is $`\mathrm{\Delta }E=\frac{N^2}{90\pi ^2}\left(\frac{15𝖳}{e}\right)^4`$. Below the transition the entropy is zero, the best variational state is the same as at zero temperature, and the average value of the Polyakov loop is zero. Above the transition, the entropy is nonzero and proportional to the number of “coloured” degrees of freedom, $`𝖲N^2`$. The average value of the Polyakov loop is nonzero and the phase is deconfined.

It is quite interesting that at high temperature our formulae numerically are quite close to the predictions of free gluon plasma. In particular, our value for the free energy, eq. (95), should be compared to the free gluon plasma expression

$$𝖥_{free}=\frac{(N^21)\pi ^2}{45}𝖳^4.$$
(103)
The ratio between the two is

$$\frac{𝖥_{free}}{𝖥_{var}}0.85.$$
(104)
The ratio of the entropies is the same.

Interestingly we get the same ratio comparing our value for the Debye mass eq.(102) with the leading order perturbative one, $`M_{pert}^2=\frac{4\pi }{3}\alpha _sNT^2`$,

$$\frac{M_{pert}^2}{M_D^2}0.85.$$
(105)
The pressure approaches its asymptotic value according to the simple formula

$$\frac{𝖯(𝖳)}{𝖯_{\mathrm{asympt}}}=1\frac{𝖳_{𝖼}^{}{}_{}{}^{4}}{𝖳^4}.$$
(106)
One has to take the comparison eqs.(104,105) with a grain of salt. As explained above, our calculations were performed assuming small $`H`$. A priori we expect that this restriction should confine us to not too large temperatures. Interestingly, however, the minimization of the free energy resulted in the value $`H/M=1/e`$ independently of temperature. Thus, we feel ourselves justified to consider the comparison eq.(104) as meaningful.

The main features of these results are indeed what we expect from the deconfinement phase transition on general grounds. It is extremely gratifying that a simple minded calculation such as ours does qualitatively so well in such a complicated problem. It therefore appears that the projection of the trial density matrix on the gauge invariant Hilbert space is, just like at zero temperature, the crucial feature that dictates most if not all the important aspects of the low energy and low temperature physics. In the context of the present calculation the most important effect of the gauge projection is obviously vanishing of the entropy in the low temperature phase. We stress that this feature was not at all built into our initial ansatz, but followed naturally and unavoidably in the disordered phase of the effective $`\sigma `$-model.

Quantitatively, our calculation of course should be taken for what it is — an approximate implementation of the variational principle. The projection over the gauge group, which as we saw is so physically important, is what makes the calculational task difficult. The most severe simplifications that we had to impose are the perturbation theory in the ordered phase of the $`\sigma `$-model and the assumption of smallness of $`H`$.

We believe that both these approximations affected the quality of our results. In particular, in the leading order of the perturbation theory the expectation value of the Polyakov loop $`U`$ is equal to unity. The actual value of $`U`$ on the ordered side of the transition according to kogut  is close to one half. Thus our perturbative calculation is rather more reliable somewhat further away from the transition. In line with this we expect that the estimate for the critical temperature we obtained here is somewhat higher than we would get, had we treated the $`\sigma `$-model more accurately in the transition region. This is consistent with the fact that our result for $`𝖳_𝖼`$ is by about $`50\%`$ higher than the lattice value of $`270\mathrm{M}\mathrm{e}\mathrm{v}`$.

The smallness of $`h`$ is also quite important. The value of $`h=1/e`$ that we obtain is in fact a reasonably small number, so omitting the corrections in powers of $`h`$ is fairly safe. On the other hand, the terms linear in $`h`$ but not enhanced by $`\mathrm{ln}h`$, which we have ignored in the present calculation, have to be accounted for in a more careful way. Once these terms are taken care of our calculation can be extended to high temperatures.

We believe that with some effort both these limitations can be overcome, at least to some extent. It should be possible to treat the $`\sigma `$-model action in a better way, perhaps along the lines of a continuum version of kogut . As for terms linear in $`h`$, those should be accessible in perturbative expansion in $`h`$, once the nonanalytic $`O(h\mathrm{ln}h)`$ terms have been understood. We plan to continue investigations along these lines and hope to improve on the quality of our results.

Notwithstanding this critique, we are very much encouraged by the results of this exploratory investigation. We believe that a reasonably simple improvement on the calculational methods has a good chance to significantly improve the results and bring them into a better quantitative agreement with the lattice data.

###### Acknowledgements.

We thank R. Buffa and J.P. Garrahan for their contribution during the early stages of this work. I.K. and A.K. thank the Royal Society for the joint project RS-CNRS grant. J.G.M. thanks the Department of Mathematics and Statistics, University of Plymouth, for hospitality while part of this work was done.