Introduction to Thermal Field Theory: From First Principles to Applications

This review provides an introduction to the fundamentals and some important applications of thermal field theory (TFT), namely, the combination of statistical mechanics and relativistic quantum field theory.

The combination of statistical mechanics and quantum field theory was first developed in the non-relativistic setting in the late 1950s to theoretically describe condensed and nuclear matter in standard laboratory experiments. This was important to bridge the gap between the theoretical description of microscopic and macroscopic systems. The formalism of second quantization, which leads to quantum field theory, can be applied within both Galilean and Einsteinian relativity, and is the most effective language to describe many-body systems, where the number of particles is allowed to change. In 1955, Matsubara [1] showed how a statistical equilibrium description of quantum field theory can be obtained by formally substituting time with an imaginary quantity. This paved the way for the imaginary-time formalism in equilibrium TFT, where the theory is described as an Euclidean quantum field theory with periodic (for bosonic fields) or antiperiodic (for fermionic fields) boundary conditions with respect to the imaginary time.

In 1965, Fradkin [2] pioneered the statistical study of relativistic quantum field theory (i.e., TFT). Several years later, it was realized that TFT is an essential tool to describe processes in the cosmological plasma. Indeed, in that context, the presence of a huge number of particles prevents the exact (non-statistical) description, and relativistic effects are non-negligible. During the subsequent decades, much work was performed to understand the nature of the electroweak phase transition using TFT, which is of great interest both in particle physics and cosmology.

Moreover, at some point during the history of the universe, when the temperature was sufficiently high, quarks and gluons behaved as weakly interacting particles, almost free to move in their plasma, the quark–gluon plasma. TFT was and still is an essential tool to describe such a physical situation. Experimental activities conducted in laboratories, such as CERN, aimed to reproduce the quark–gluon plasma through collisions of heavy ions, such as gold or lead nuclei. The application of TFT to heavy-ion collisions further boosted the interest in this theory. The study of the frontier (the QCD transition) between the regime where quarks and gluons form such a plasma and the regime where they are confined in colorless bound particles needs non-perturbative methods. It was then natural to combine lattice gauge theory and TFT. The first studies of this type were conducted in the beginning of the 1980s. In recent years, several lattice studies have agreed on the evidence that the QCD transition temperature is about 156.5 MeV at negligible densities [3].

Today, TFT remains the standard theoretical tool to study particle physics processes (decays, scattering processes, particle production, phase transitions, etc.) in a thermal medium. In this review, we will provide several examples of works along these lines. In cosmology, it is believed that a thermal medium was originally formed shortly after a primordial exponential expansion of the spacetime volume (i.e., inflation) through a mechanism called reheating.

The aim of this review is not to give an exhaustive overview of the TFT literature but is rather to introduce in a constructive way the fundamentals of TFT starting from the basic principles of quantum mechanics and relativity. Once this is performed, it appears natural to provide some of the most important physical applications of TFT that can already been understood. Here, two main applications are discussed: particle production in a thermal bath and the theory of phase transitions in field theory. To achieve these goals, a detailed treatment of the thermal Green’s functions is provided. These functions, whose importance in TFT goes even beyond the applications considered here, are defined as the statistical average of the expectation values of the time-ordered product of fields taken on a complete set of states. Both the real-time and imaginary-time formalisms are explained, as they are, respectively, the most suitable approaches for describing particle production in a thermal bath and phase transitions.

As any other review (for other reviews and textbooks related to the subjects discussed here, see [4,5,6,7,8,9,10] (see also Ref. [11] for a related discussion)), most of the results presented here are not original, but actually some derivations and explanations are first provided here. These include, for example, the derivation and interpretation of the full equilibrium density matrix given in Section 2, the explanation, given in the same section, of the general relevance in TFT of the thermal Green’s functions, which are then studied in Section 4, the full derivation of the path-integral representation of these thermal functions in the fermionic case, and the description of vacuum decay for a generic first-order phase transition provided in Section 6.5.

The writing of this review originally began as part of the Ph.D. lectures on TFT given by the author in the spring of 2023 at the University of Rome Tor Vergata. No previous knowledge of TFT is required to understand this monograph, but a good knowledge of quantum field theory without statistical mechanics is assumed.

In quantum field theory (QFT), one usually focuses on processes with a small number of particles. However, there are many physical situations where one encounters processes in a medium such as, for example, the propagation of particles and scattering processes in a gas at finite temperature and density. Of course, microscopically, such situations involve a very large number of particles.

If the system is initially in a pure state $|\alpha \rangle$, we can still write the transition amplitude for a scattering process as where $|\beta \rangle$ is the final state in question (we take as usual $|\alpha \rangle$ and $|\beta \rangle$ with unit norm without loss of generality) and $\widehat{S}$ is the scattering operator, which transforms $|\alpha \rangle$ into the state (represented in the interaction picture) at time $+\infty$. The evaluation of ${\mathcal{A}}_{\alpha \to \beta }$ may be computationally difficult, but Equation (1) provides the right formula for the transition amplitude in general.

However, if the initial state is only statistically specified, we are interested in averaging the rates of the processes over initial states with appropriate probabilities $p_{\alpha }$ (or over final states with probabilities $p_{\beta }$). It is important to note that these are the probabilities of choosing particular states in a statistical ensemble and are not the quantum probabilities of the collapse of the wave function onto a specific eigenstate of the observable being measured.

If we want to compute the rates of inclusive processes (such as decay or scattering rates summed over final states or production rates summed over initial states), we can use the optical theorem, and express them in terms of the imaginary part of the $\alpha \to \alpha$ amplitude, where $\widehat{T}\equiv -i(\widehat{S}-1)$, and the unitarity condition $1={\widehat{S}}^{†}\widehat{S}=1+i(\widehat{T}-{\widehat{T}}^{†})+{\widehat{T}}^{†}\widehat{T}$ has been used, or in terms of the $\beta \to \beta$ amplitude where the unitarity condition in the form $1=\widehat{S}{\widehat{S}}^{†}=1+i(\widehat{T}-{\widehat{T}}^{†})+\widehat{T}{\widehat{T}}^{†}$ has been used. If the system is only statistically specified, we should sum over $\alpha$ with weights $p_{\alpha }$ (using (2)) or equivalently sum over $\beta$ with weights $p_{\beta }$ (using (3))

From these results, it is clear that in TFT, important quantities are the so-called thermal Green’s functions for the field operators $\Phi (x)$, which depend on n spacetime variables $x_1,\dots ,x_n$. Here, $\mathcal{T}$ is the time-ordered product. The functions in (6) are also known as n-point functions and the one with $n=2$ is called the thermal propagator.

In quantum statistical mechanics, the above-mentioned probabilities are given by the density matrix (which is actually an operator) where the unit-norm states $|\alpha \rangle$ form an orthonormal basis, and $p_{\alpha }$ is the probability of finding $|\alpha \rangle$ in the statistical ensemble. The expectation value of an operator A with this statistical distribution is which, using the orthonormality of the basis $\left\{|\alpha \rangle \right\}$, can be rewritten as

Note that the thermal Green’s function, Equation (6), can be written as

Since the $p_{\alpha }$ represent a probability distribution, we have that all $p_{\alpha }$ are real and

The reality of $p_{\alpha }$ implies the hermiticity of $\rho$. Note that choosing a statistical probability distribution $p_{\alpha }$ for the states $|\alpha \rangle$ is equivalent to choosing a density matrix $\rho$. Moreover, note that $p_{\alpha }$ can be written as $exp(-h_{\alpha })$ where $h_{\alpha }$ is a non-negative real number. So we can write where $\mathcal{H}$ is a Hermitian operator with real eigenvalues (which are the $h_{\alpha }$). The trace is independent of the basis, so Equation (9) and $\mathrm{Tr}\rho =1$ hold in any basis. The case of a pure state $|\lambda \rangle$ corresponds to $p_{\lambda }=1$ and $p_{\alpha }=0$ for all $\alpha \ne \lambda$, and in this case ${\rho }^2=\rho$. The equation ${\rho }^2=\rho$ is not only a necessary but also a sufficient condition for the system to be in a pure state. Indeed, for ${\rho }^2=\rho$ where in the last step, the orthonormality of the basis $\left\{|\alpha \rangle \right\}$ has been used. Since the $p_{\alpha }$ form a probability distribution, Equation (13) tells us that one of the $p_{\alpha }$ is one, and all the others are zero, namely, the system is in a pure state. Therefore, for $\rho \ne {\rho }^2$, we necessarily have a non-trivial mixing of states.

From Equation (7), it follows that under a generic unitary symmetry operator U, the density matrix transforms as which resembles but is different from the usual transformation of a generic operator $\mathcal{O}$ introduced in basic courses on quantum mechanics, i.e., $\mathcal{O}\to U^{†}\mathcal{O}U$ because $U^{†}$ in (14) is on the right rather than on the left. Suppose now that at time $t_0$, the density matrix $\rho (t_0)$ is given by (7). Then, at time t, the density matrix is where $U(t)=exp(-iH(t-t_0))$ is the time-evolution operator. So $\rho (t)$ satisfies the time-evolution equation where a dot represents the derivative with respect to t.

Let us suppose now that the system is in thermal equilibrium. In this case, $\rho$ must be independent of time and so $[H,\rho ]=0$. This equation tells us that the density matrix is a function of the form $\rho =\rho ({\mathcal{O}}_1,{\mathcal{O}}_2\dots )$, where the ${\mathcal{O}}_i$ are operators that commute with the Hamiltonian, namely, they are conserved quantities. An important piece of information on this function can be obtained by considering additively conserved quantities ${\mathcal{O}}_i$, such as the energy itself, the (linear) momentum, the angular momentum, the charges, etc. Now, if one divides the system under study into independent subsystems I, II, …, the probability distribution is the product of the probability distribution of I, II, …, and so the density matrix is the direct product of the density matrices of I, II, …, where $\rho ={\rho }_{\mathrm{I}}\times {\rho }_{\mathrm{II}}\times \dots$. On the other hand, since ${\mathcal{O}}_i$ are additively conserved quantities, they can be expressed as the direct sum of conserved quantities for the subsystems, ${\mathcal{O}}_i={\mathcal{O}}_{i\mathrm{I}}+{\mathcal{O}}_{i\mathrm{II}}+\dots$. (To show that there no other solutions, consider the density matrix, ${\rho }_0(\mathcal{O})$, with another normalization, i.e., ${\rho }_0(0)=1$ (here, the index i that labels the conserved quantities is understood for simplicity). This is related to the density matrix $\rho (\mathcal{O})$ with normalization Tr$\rho$ = 1 through $\rho ={\rho }_0/$Tr${\rho }_0\equiv {\rho }_0/Z$. Now the conditions above give ${\rho }_0({\mathcal{O}}_I+{\mathcal{O}}_{II}+\dots )={\rho }_0({\mathcal{O}}_I)\times {\rho }_0({\mathcal{O}}_{II})\times \dots$, which implies ${\rho }_0=exp(-\beta \mathcal{O})$ and so $\rho =exp(-\beta \mathcal{O})/Z$. Indeed, for generic commuting variables x and h, the condition ${\rho }_0(x+h)={\rho }_0(x){\rho }_0(h)$ tells us which gives ${\rho }_0(x)=exp(-\beta x)$, where $\beta \equiv -{\rho }_0^{\prime }(0)$. The only function of the ${\mathcal{O}}_i$ with this property is the exponential (In these notes repeated indices understand a summation (unless otherwise stated).), where Z and the ${\beta }_i$ are constants. If the ${\mathcal{O}}_i$ are Hermitian, the ${\beta }_i$ and thus Z are real because, as we have seen, we can write $\rho$ as in (12). From $\mathrm{Tr}\rho =1$ it follows which is known as the partition function.

(In this review, we will only consider special relativity, leaving aside the (nevertheless interesting) study of gravity and non-inertial frames. These can lead to interesting effects such as the Unruh one [12].) In the relativistic case, the conserved quantities are H, the (linear) momentum $\overrightarrow{P}$, the angular momentum $\overrightarrow{J}$ and, possibly, a set of charges $Q_a$, which generate possible continuous internal symmetry groups. Thus, the equilibrium density matrix can be written as where $P^{\mu }$ is the four-momentum, $J_{ij}={ϵ}_{ijk}{J}_k$, the $J_k$ are the three components of the angular momentum $\overrightarrow{J}$, and ${ϵ}_{ijk}$ is the totally antisymmetric Levi–Civita symbol (we choose ${ϵ}_{123}=1$). We take ${\beta }_{ij}=-{\beta }_{ji}$ without loss of generality because $J_{ij}=-J_{ji}$ (a symmetric part of ${\beta }_{ij}$ would not contribute to (20)).

Let us consider now a covariant generalization of (20) where the $J^{\mu \nu }$ are the generators of the full Lorentz group. Recalling $J_{\mu \nu }=-J_{\nu \mu }$, we take ${\beta }_{\mu \nu }=-{\beta }_{\nu \mu }$ without loss of generality. The density matrix in (21) does not generically describe a system at equilibrium because some generators of the Lorentz group, the boosts $J_{0i}$, do not commute with H. However, this density matrix allows us to easily obtain useful information on $\langle P^{\mu }\rangle$ and $\langle J^{\mu \nu }\rangle$ and so, as particular cases, on $\langle \overrightarrow{P}\rangle$ and $\langle \overrightarrow{J}\rangle$. Let us consider a proper (in these notes, a proper Lorentz transformation has by definition ${\Lambda }_0^0\ge 1$ and det$\Lambda =1$, where $\Lambda$ is the Lorentz matrix (with elements in the $\mu$th line and $\nu$th column given by ${\Lambda }_{\nu }^{\mu }$)) Lorentz transformation $x^{\mu }\to {\Lambda }_{\nu }^{\mu }{x}^{\nu }$ (the $x^{\mu }$ are the spacetime coordinates), and let $U(\Lambda )$ be the corresponding unitary operator on the Hilbert space. According to the general rule (14), the density matrix in (21) transforms as where the Lorentz invariance of charges, namely, has been used. Now $P^{\mu }$ and $J^{\mu \nu }$ transform respectively as a vector and a rank-two tensor under proper Lorentz transformations, i.e.,

But in (22), $U^{†}$ appears on the right rather than on the left as a consequence of the general rule in (14), so we need the inverse transformation rules where we use that the elements in the $\mu$th line and $\nu$th column of the inverse Lorentz matrix ${\Lambda }^{-1}$ are and ${\eta }_{\mu \nu }$ is the flat metric, which is used as usual to lower and raise the spacetime indices. From (25), it follows that the Lorentz transformation of the density matrix in (22) can be written as

Therefore, performing this Lorentz transformation of the density matrix is equivalent to replacing ${\beta }_{\mu }\to {\Lambda }_{\mu }^{\rho }{\beta }_{\rho }$ and ${\beta }_{\mu \nu }\to {\Lambda }_{\mu }^{\rho }{\Lambda }_{\mu }^{\sigma }{\beta }_{\rho \sigma }$. Lorentz covariance then tells us that where $f_1$, $f_2$ and $f_J$ are functions of the Lorentz-invariant quantities ${\beta }_a$, ${\beta }_{\mu }{\beta }^{\mu }$, ${\beta }_{\mu \nu }{\beta }^{\mu \nu }$, ${\beta }_{\mu \nu }{\beta }^{\nu }{\beta }_{\rho }{\beta }^{\mu \rho }$ only.

Taking now the generic equilibrium density matrix in (20), we have to set ${\beta }_{0i}=0$ that gives

In particular, for the components of the momentum, and we see that setting ${\beta }^i=0$ corresponds to having a system with zero (average) momentum. Moreover, the second equation in (29) indicates that zero (average) angular momentum corresponds to ${\beta }_{ij}=0$.

From now on, we set ${\beta }_{ij}=0$ and consider a system at equilibrium without angular momentum. In this case, we simply have

The convergence of the traces in (9) and (11) requires ${\beta }_0>0$ and ${\beta }_i{\beta }_i<{\beta }_0^2$. This can be easily seen by choosing a basis of momentum, energy, and charge eigenstates $|n\rangle$ to compute the traces such that where $p_n^{\mu }$ (and $C_n$) is the eigenvalue of $P^{\mu }$ (and, respectively, ${\beta }_a{Q}_a$) corresponding to the eigenstate $|n\rangle$. Focusing on the terms in (32) with zero spatial momentum and charge, one can easily see that the convergence of the sum requires ${\beta }_0>0$. Moreover, if we had ${\beta }_i{\beta }_i\ge {\beta }_0^2$, there would be directions in the space of states where the sum in (32) does not converge (for large momenta, one can neglect the masses, and the absolute value of the momentum can approach the total energy and also can feature a negative $p_n^i{\beta }_i$ of the form $-|{\overrightarrow{p}}_n||\overrightarrow{\beta }|$ even for vanishing charges). Therefore, there exists a proper Lorentz transformation that transforms ${\beta }_i\to 0$ and ${\beta }_0\to \beta$, where $\beta$ is some positive number. According to (30), this transformation corresponds to moving to the rest frame of the system. Therefore, in the rest frame, an equilibrium density matrix reads

The inverse of $\beta$ is the temperature T, while ${\mu }_a\equiv -{\beta }_a/\beta$ is the chemical potential associated with $Q_a$.

Let us start studying what happens in a field theory by considering the simple case of free fields. The generalization of the interacting case is provided in Section 4.