Field theory, but hot. · Donkey on the Edge

QFT at nonzero temperature: the Matsubara formalism, thermal propagators, and the machinery behind the quark-gluon plasma and the early universe.

Every QFT document so far has assumed zero temperature; computing vacuum expectation values, scattering amplitudes, Green’s functions in the vacuum state. But the universe isn’t at zero temperature. Nuclear collisions at RHIC and LHC produce matter at $T \sim 2\times 10^{12}$ K (200 MeV). The early universe went through electroweak symmetry breaking at $T \sim 100$ GeV and QCD confinement at $T \sim 170$ MeV. Inside neutron stars, matter reaches extreme temperatures and densities. Understanding any of these requires quantum field theory at finite temperature.

Fortunately, the framework extends naturally. The key insight, due to Matsubara (1955), is that Euclidean QFT with a compactified time direction of period $\beta = 1/T$ is thermal QFT. The formal equivalence between quantum fluctuations and thermal fluctuations that we’ve hinted at since document 9 becomes a computational tool.

This document develops thermal field theory from the Matsubara formalism through to applications. By the end, you’ll understand how modern particle physics calculations handle hot systems.

Prerequisites

Statistical mechanics reference (ensembles, partition functions)
QFT documents 1-12 (especially 9 for path integrals, 8 for RG)
Condensed matter document (for the analog of phase transitions)

Conventions

Mostly-minus metric
$\hbar = c = k_B = 1$
Temperature $T$ measured in energy units (so $\beta = 1/T$ , with $\beta$ in inverse energy)
At $T = 0$ , we recover vacuum QFT from previous documents

Why Thermal QFT?
The Fundamental Equivalence
The Matsubara Formalism
Imaginary-Time Propagators
Matsubara Frequency Sums
Thermal Feynman Rules
Worked Example: The Thermal Scalar Self-Energy
Real-Time Formalism: When Imaginary Time Isn’t Enough
Hard Thermal Loops and Effective Thermal Masses
Symmetry Restoration at High Temperature
Application: The Quark-Gluon Plasma
Application: The Early Universe
The Sign Problem, Revisited
Appendix: Thermal Formulas Reference

1. Why Thermal QFT?

The Physical Motivation

At $T = 0$ , a quantum field theory has a unique vacuum state $|0\rangle$ , and all correlation functions are vacuum expectation values. This is what we’ve computed so far.

At $T > 0$ , the system is in a thermal ensemble; a statistical mixture of energy eigenstates, weighted by $e^{-\beta E}$ . Correlation functions become thermal averages:

$\langle\mathcal{O}\rangle_\beta = \frac{1}{Z}\text{tr}[e^{-\beta\hat H}\mathcal{O}]$

where $Z = \text{tr}(e^{-\beta\hat H})$ is the partition function.

This is formally identical to the statistical mechanics of a quantum system. What makes it “field theory” is that the Hamiltonian $\hat H$ describes a relativistic quantum field (infinite degrees of freedom, Lorentz-invariant interactions).

Where Does Thermal QFT Matter?

Heavy-ion collisions. RHIC and LHC smash gold or lead nuclei together at ultra-relativistic energies, producing fireballs of deconfined quarks and gluons at $T \sim 200-500$ MeV. Understanding this “quark-gluon plasma” requires thermal QCD.

The early universe. In the first $\sim 10^{-10}$ seconds after the Big Bang, the universe was hotter than the electroweak scale. Electroweak symmetry was restored, the Higgs hadn’t “frozen” to its current vacuum, and baryogenesis (the generation of matter-antimatter asymmetry) may have occurred. All thermal QFT.

Neutron stars and compact objects. Dense matter at $T\sim 10$ MeV and huge chemical potential. The equation of state is determined by strong-coupling QCD at extreme conditions.

Condensed matter systems. Many condensed matter problems are naturally formulated as thermal QFTs; superconductors, superfluids, and strongly correlated systems at nonzero temperature. The same Matsubara formalism works there.

Cosmological phase transitions. The electroweak transition and the QCD transition in the early universe were quantum field theory phase transitions at finite temperature. Their character (first-order vs. crossover, latent heat, bubble nucleation) is thermal QFT.

The Big Mathematical Idea

The partition function $Z = \text{tr}(e^{-\beta\hat H})$ can be written as a path integral over fields that are periodic in imaginary time with period $\beta$ :

$Z = \int_{\phi(\tau = 0) = \phi(\tau = \beta)}\mathcal{D}\phi\, e^{-S_E[\phi]}$

This one equation; that thermal QFT equals Euclidean QFT on a circle of circumference $\beta$ ; is the foundation of everything that follows.

2. The Fundamental Equivalence

The Derivation

Start with the quantum thermal partition function:

$Z(\beta) = \text{tr}(e^{-\beta\hat H}) = \sum_n\langle n|e^{-\beta\hat H}|n\rangle$

where $|n\rangle$ runs over a complete set of states.

For a field theory, use position eigenstates (eigenstates of the field operator $\hat\phi(\vec x)$ at a fixed time):

$Z(\beta) = \int\mathcal{D}\phi_0\,\langle\phi_0|e^{-\beta\hat H}|\phi_0\rangle$

where $|\phi_0\rangle$ is an eigenstate of $\hat\phi(\vec x)$ with eigenvalue $\phi_0(\vec x)$ (a field configuration at one instant).

Connection to the Path Integral

The matrix element $\langle\phi_0|e^{-\beta\hat H}|\phi_0\rangle$ is precisely an amplitude to evolve from $\phi_0$ at imaginary time 0 to $\phi_0$ at imaginary time $\beta$ . Recall that at $T = 0$ :

$\langle\phi_f|e^{-i\hat H(t_f - t_i)}|\phi_i\rangle = \int_{\phi(t_i) = \phi_i}^{\phi(t_f) = \phi_f}\mathcal{D}\phi\, e^{iS[\phi]}$

Substituting $t_f - t_i = -i\beta$ (Wick rotation with imaginary time extending from 0 to $\beta$ ), and taking $\phi_i = \phi_f = \phi_0$ :

$\langle\phi_0|e^{-\beta\hat H}|\phi_0\rangle = \int_{\phi(0) = \phi_0}^{\phi(\beta) = \phi_0}\mathcal{D}\phi\, e^{-S_E[\phi]}$

where $S_E$ is the Euclidean action with $\tau \in [0, \beta]$ .

Putting It Together

Now integrate over $\phi_0$ (the trace):

$Z(\beta) = \int\mathcal{D}\phi_0\int_{\phi(0) = \phi_0}^{\phi(\beta) = \phi_0}\mathcal{D}\phi\, e^{-S_E[\phi]}$

The outer integral over $\phi_0$ combined with the constrained inner integral is equivalent to integrating over all fields that satisfy $\phi(0) = \phi(\beta)$ ; i.e., fields that are periodic in imaginary time with period $\beta$ :

$\boxed{Z(\beta) = \int_{\phi(\tau) = \phi(\tau + \beta)}\mathcal{D}\phi\, e^{-S_E[\phi]}}$

This is the fundamental formula of thermal field theory.

Geometric Picture

The spacetime on which the theory is defined has changed. At $T = 0$ , we work on 4D Euclidean space $\mathbb{R}^4$ . At finite $T$ , we work on $\mathbb{R}^3 \times S^1_\beta$ ; three-dimensional infinite space times a circle of circumference $\beta$ in the time direction.

The compactification of imaginary time into a circle is what encodes the thermal state. Fields must be periodic (for bosons) or antiperiodic (for fermions; see below) around this circle.

Fermions: A Twist

For bosons, fields are periodic in $\tau$ with period $\beta$ . For fermions, fields are antiperiodic:

$\psi(\tau + \beta) = -\psi(\tau)$

Why? The trace in $Z = \text{tr}(e^{-\beta\hat H})$ for fermions involves anticommuting fermion operators. When you derive the path integral representation, the spin-statistics structure flips the sign of the boundary condition. This antiperiodicity has dramatic consequences for fermionic Matsubara frequencies (section 3).

What Changed from Vacuum QFT

Compared to vacuum QFT:

Time is compactified: finite interval $[0, \beta]$ instead of $(-\infty, \infty)$
Fields are periodic (bosons) or antiperiodic (fermions) in $\tau$
Everything else (the Lagrangian, interactions, couplings) is the same

This geometric change is what makes thermal physics different from vacuum physics. The propagators, loop integrals, and physical predictions all differ because of it.

3. The Matsubara Formalism

Fourier Decomposition in Imaginary Time

Because $\tau \in [0, \beta]$ is compact, the Fourier modes in this direction are discrete. Instead of a continuous frequency $\omega$ , we have discrete Matsubara frequencies $\omega_n$ .

For bosons (periodic):

$\phi(\tau, \vec x) = \sum_{n = -\infty}^{\infty}\int\frac{d^3 k}{(2\pi)^3}e^{-i\omega_n\tau + i\vec k\cdot\vec x}\tilde\phi_n(\vec k)$

Periodicity $\phi(\tau + \beta) = \phi(\tau)$ requires $e^{-i\omega_n\beta} = 1$ , so:

$\omega_n = \frac{2\pi n}{\beta} = 2\pi nT, \quad n \in \mathbb{Z}$

These are bosonic Matsubara frequencies: $\omega_n \in \{0, \pm 2\pi T, \pm 4\pi T, \ldots\}$ .

For fermions (antiperiodic):

$\psi(\tau + \beta) = -\psi(\tau)$

Requires $e^{-i\omega_n\beta} = -1$ , so:

$\omega_n = \frac{(2n+1)\pi}{\beta} = (2n+1)\pi T, \quad n \in \mathbb{Z}$

These are fermionic Matsubara frequencies: $\omega_n \in \{\pm\pi T, \pm 3\pi T, \pm 5\pi T, \ldots\}$ .

Notice: no zero mode for fermions. This is crucial; it’s why fermionic effects are often exponentially suppressed at low temperature ( $\propto e^{-m/T}$ for the lowest Matsubara mode with $\omega_0 = \pi T$ ).

The Key Replacements

Going from vacuum to thermal QFT:

$\int\frac{dp_0}{2\pi} \to T\sum_{n}$

(Continuous frequency integral → discrete Matsubara sum, with factor of $T$ from the measure.)

And the Feynman rules change:

$p_0 \to i\omega_n$

The external momentum in the time direction becomes $i\omega_n$ (imaginary, because we’re in Euclidean signature). The " $i$ " arises from Wick rotation.

Temperature as Inverse Radius

A useful intuition: the temperature $T$ is the inverse circumference of the imaginary-time circle. High $T$ means a small circle; imaginary time is compactified to a short interval. Low $T$ (large $\beta$ ) means a big circle; imaginary time is nearly infinite, approaching $T = 0$ .

In the limit $T \to 0$ : $\beta \to \infty$ , the Matsubara sum becomes an integral (the frequencies become dense), and we recover vacuum QFT.

In the limit $T \to \infty$ : $\beta \to 0$ , only the $n = 0$ Matsubara mode dominates for bosons. The system becomes effectively 3-dimensional (time “shrinks away”). This is called dimensional reduction and underlies many high-temperature simplifications.

Why This Works

The Matsubara formalism is mathematically rigorous and computationally efficient. All the machinery of perturbation theory; Feynman diagrams, renormalization, RG; extends to the thermal case with:

Integrals over $p_0$ replaced by sums over $\omega_n$
Propagators modified to account for the discrete frequencies
Overall factors of $\beta$ from the time integration

Nothing is lost; we’re just computing in a slightly different geometry.

4. Imaginary-Time Propagators

The Free Scalar Propagator

The Euclidean scalar propagator at $T = 0$ :

$D_E(x - y) = \int\frac{d^4p_E}{(2\pi)^4}\frac{e^{ip_E\cdot(x - y)}}{p_E^2 + m^2}$

At $T > 0$ , the time integration becomes a Matsubara sum:

$D_T(\tau, \vec x) = T\sum_n\int\frac{d^3k}{(2\pi)^3}\frac{e^{-i\omega_n\tau + i\vec k\cdot\vec x}}{\omega_n^2 + \vec k^2 + m^2}$

with $\omega_n = 2\pi nT$ (bosonic).

The Matsubara Sum

The sum over $n$ can be evaluated using the identity:

$T\sum_{n=-\infty}^\infty\frac{1}{\omega_n^2 + E_k^2} = \frac{1}{2E_k}\coth\left(\frac{\beta E_k}{2}\right)$

where $E_k = \sqrt{\vec k^2 + m^2}$ and the $\coth$ encodes the thermal occupation.

Using $\coth(x/2) = \frac{e^{x/2} + e^{-x/2}}{e^{x/2} - e^{-x/2}} = 1 + \frac{2}{e^x - 1}$ :

$\coth(\beta E_k/2) = 1 + 2 n_B(E_k)$

where $n_B(E) = 1/(e^{\beta E} - 1)$ is the Bose-Einstein distribution.

So:

$D_T(\tau = 0, \vec x) = \int\frac{d^3k}{(2\pi)^3}\frac{e^{i\vec k\cdot\vec x}}{2E_k}[1 + 2n_B(E_k)]$

The “1” is the vacuum contribution; the “2 $n_B$ ” is the thermal contribution.

Decomposition: Vacuum + Thermal

This decomposition is general and important. Any thermal propagator can be written as:

$G_T(p) = G_0(p) + G_{\rm th}(p)$

where $G_0$ is the vacuum propagator and $G_{\rm th}$ vanishes at $T = 0$ .

For momentum-space propagators after Wick rotation back to real frequencies:

$G_T(p^0, \vec p) = \frac{i}{(p^0)^2 - E_p^2 + i\epsilon} + 2\pi n_B(|p^0|)\delta((p^0)^2 - E_p^2)$

The first term is the usual Feynman propagator. The second term is the thermal piece; it puts particles on-shell with weight given by the occupation number.

This decomposition has a physical interpretation:

Vacuum piece: quantum fluctuations (virtual particle pairs)
Thermal piece: real thermal particles that happen to be present in the ensemble

Loop integrals in thermal QFT involve both.

Fermion Thermal Propagator

For fermions with antiperiodic boundary conditions, the Matsubara frequencies are $\omega_n = (2n+1)\pi T$ . The analogous identity:

$T\sum_n\frac{1}{\omega_n^2 + E^2} = \frac{1}{2E}\tanh(\beta E/2) = \frac{1}{2E}[1 - 2n_F(E)]$

where $n_F(E) = 1/(e^{\beta E} + 1)$ is the Fermi-Dirac distribution. The ”-” sign flips because of the antisymmetric tanh instead of symmetric coth.

So the fermion thermal propagator decomposes as:

$S_T(p) = S_0(p) - 2\pi n_F(|p^0|)(\slashed p + m)\delta((p^0)^2 - E_p^2)$

The minus sign in front of the thermal piece is what encodes Fermi-Dirac statistics at the level of propagators.

Photon/Gluon Thermal Propagator

For gauge bosons, the thermal propagator has the same structure as the scalar, but multiplied by the appropriate tensor structure (transverse projector, with gauge choices affecting the longitudinal piece). In Feynman gauge:

$D^{\mu\nu}_T(k) = -\eta^{\mu\nu}\left[\frac{i}{k^2 + i\epsilon} + 2\pi n_B(|k^0|)\delta(k^2)\right]$

At finite temperature, gauge bosons have thermal distributions just like matter particles.

Useful Limits

Low temperature $T \ll m$ : The thermal piece is exponentially suppressed: $n_B(m) \sim e^{-m/T}$ . Thermal effects are negligible, and we recover vacuum physics.

High temperature $T \gg m$ : The thermal piece dominates. Both $n_B(E)$ and $n_F(E)$ become $O(1)$ for $E \lesssim T$ , and thermal particles fill the phase space up to momenta $\sim T$ .

5. Matsubara Frequency Sums

Evaluating thermal loop integrals requires summing over Matsubara frequencies. There’s a beautiful contour-integration technique that makes this systematic.

The Contour Trick

The key identity: a Matsubara sum can be rewritten as a contour integral using the poles of $\coth$ or $\tanh$ .

Bosonic case. The function $\coth(\beta p^0/2) = 1 + 2n_B(p^0)$ has simple poles at $p^0 = i\omega_n = 2\pi niT$ for all integers $n$ , with residue $2/\beta = 2T$ .

So for any function $f(p^0)$ analytic at these poles:

$T\sum_{n=-\infty}^\infty f(i\omega_n) = \frac{1}{2\pi i}\oint_C f(p^0)\cdot\frac{\beta}{2}\coth\left(\frac{\beta p^0}{2}\right) dp^0$

where $C$ is a contour encircling all the imaginary-axis poles.

Fermionic case. The function $\tanh(\beta p^0/2) = 1 - 2n_F(p^0)$ has simple poles at $p^0 = i\omega_n = (2n+1)\pi iT$ , with residue $2/\beta = 2T$ .

So:

$T\sum_{n=-\infty}^\infty f(i\omega_n) = \frac{1}{2\pi i}\oint_C f(p^0)\cdot\frac{\beta}{2}\tanh\left(\frac{\beta p^0}{2}\right) dp^0$

Deforming the Contour

The contour $C$ can be deformed away from the imaginary axis. Typically, you close it around the poles of $f(p^0)$ (which are the physical poles of the propagator). Using the residue theorem:

$T\sum_n f(i\omega_n) = -\sum_{\rm poles\ of\ }f\text{Res}\left[f(p^0)\cdot\frac{\beta}{2}\coth\left(\frac{\beta p^0}{2}\right)\right]$

(Or $\tanh$ for fermions.)

The sign comes from deforming the contour.

Example: The Scalar Bubble

Consider the Matsubara sum for a simple loop:

$I(\vec p) = T\sum_n\int\frac{d^3k}{(2\pi)^3}\frac{1}{(\omega_n^2 + E_k^2)((\omega_n)^2 + E_{p-k}^2)}$

Wait, this is more complex. Let me do a simpler one:

Simpler example. Sum over $n$ of $1/(\omega_n^2 + E^2)$ for a single energy $E$ .

The sum has poles at $p^0 = \pm iE$ (using $p^0 \to i\omega_n$ analytic continuation). Applying the contour formula:

$T\sum_n\frac{1}{\omega_n^2 + E^2} = -\text{Res}_{p^0 = iE}\left[\frac{1}{-(p^0)^2 + E^2}\cdot\frac{\beta}{2}\coth(\beta p^0/2)\right] - \text{Res}_{p^0 = -iE}[\ldots]$

Computing the residue at $p^0 = iE$ :

$\frac{1}{-(p^0)^2 + E^2} = \frac{1}{(E - p^0)(E + p^0)}$

Pole at $p^0 = iE$ : the factor $(E - p^0)$ vanishes. Near the pole: $1/((E - p^0)(E + iE)) = 1/(2iE(E - p^0))$ .

Residue of $[\beta/2\coth(\beta p^0/2)/(E - p^0)]$ at $p^0 = iE$ : $\beta\coth(i\beta E/2)/(2iE)\cdot(-1) = i\beta\cot(\beta E/2)/(2E)\cdot(-1)$

This is getting tangled. Let me just cite the final result from the textbook:

$T\sum_n\frac{1}{\omega_n^2 + E^2} = \frac{1}{2E}\coth(\beta E/2)$

(I did this in section 4.)

What the Sum Captures

The Matsubara sum with coth (bosons) or tanh (fermions) automatically includes:

Vacuum contribution (the “1” in $\coth = 1 + 2n_B$ )
Thermal contribution (the $n_B$ or $n_F$ pieces)
Correct signs for bosons vs. fermions

This is why the Matsubara formalism is so powerful: one sum captures all the thermal physics in one shot. Loop integrals in thermal QFT look just like vacuum ones, with the Matsubara sum replacing the frequency integral.

6. Thermal Feynman Rules

Putting everything together, the Feynman rules for thermal QFT are:

Propagators

Scalar boson: $1/(\omega_n^2 + \vec k^2 + m^2)$ where $\omega_n = 2\pi nT$

Fermion: $1/(-i\omega_n\gamma^0 - \vec\gamma\cdot\vec k + m)$ where $\omega_n = (2n+1)\pi T$

Gauge boson (Feynman gauge): $-\eta^{\mu\nu}/(\omega_n^2 + \vec k^2)$

Vertices

Same as vacuum; thermal effects don’t change the local interactions, only the propagator structure and the loop measure.

Loop Integrations

Replace:

$\int\frac{d^4p}{(2\pi)^4} \to T\sum_{n}\int\frac{d^3 k}{(2\pi)^3}$

Each loop contributes a factor of $T$ (from the conversion of the time-integral to a sum).

Overall Structure

A thermal Feynman diagram calculation goes:

Write down the standard diagram with appropriate external lines
Replace frequency integrals with Matsubara sums
Use propagator structure with $p_0 \to i\omega_n$
Evaluate the Matsubara sums using contour tricks or direct computation
Do the remaining 3-momentum integrals
Take the answer: usually has vacuum piece + thermal piece, with the thermal piece involving distribution functions

The vacuum piece may need renormalization (same UV divergences as at $T = 0$ ). The thermal piece is typically UV-finite (because $n_B$ , $n_F$ exponentially suppress high momenta).

This is one of the beautiful features of finite-temperature field theory: thermal corrections are automatically UV-convergent. Temperature acts as a natural cutoff.

Example: A Simple Scaling Argument

For a generic one-loop thermal integral at high temperature:

$T\sum_n\int\frac{d^3k}{(2\pi)^3}\frac{1}{\omega_n^2 + E_k^2}\cdot(\text{something})$

The thermal contribution from the $n_B$ piece scales as:

$\int\frac{d^3k}{(2\pi)^3}n_B(E_k) \sim \int\frac{d^3k}{(2\pi)^3}\frac{1}{e^{k/T} - 1} \sim T^3\int\frac{d^3x}{(2\pi)^3}\frac{1}{e^x - 1} \sim T^3$

So thermal loop corrections scale as $T^3$ (times whatever powers of coupling).

This is the origin of thermal masses: mass-squared contributions $\propto g^2 T^2$ from loop effects. We’ll see this explicitly in section 9.

7. Worked Example: The Thermal Scalar Self-Energy

Let’s compute the one-loop thermal mass for a scalar in $\phi^4$ theory.

Setup

The theory: $\mathcal{L} = \tfrac{1}{2}(\partial\phi)^2 - \tfrac{1}{2}m^2\phi^2 - \tfrac{\lambda}{4!}\phi^4$

The one-loop tadpole diagram contributes to the scalar self-energy:

$\Sigma_{\rm tadpole}(p) = \frac{\lambda}{2}T\sum_n\int\frac{d^3k}{(2\pi)^3}\frac{1}{\omega_n^2 + \vec k^2 + m^2}$

Note: the factor of $1/2$ is the symmetry factor for the tadpole.

Step 1: Evaluate the Matsubara Sum

Using the identity from section 4:

$T\sum_n\frac{1}{\omega_n^2 + E_k^2} = \frac{1}{2E_k}\coth(\beta E_k/2) = \frac{1}{2E_k}[1 + 2n_B(E_k)]$

So:

$\Sigma_{\rm tadpole} = \frac{\lambda}{2}\int\frac{d^3k}{(2\pi)^3}\frac{1}{2E_k}[1 + 2n_B(E_k)]$

$= \frac{\lambda}{4}\int\frac{d^3k}{(2\pi)^3}\frac{1}{E_k} + \frac{\lambda}{2}\int\frac{d^3k}{(2\pi)^3}\frac{n_B(E_k)}{E_k}$

Step 2: Split into Vacuum and Thermal

The first term is the usual vacuum contribution (divergent, renormalized away). The second term is the thermal correction:

$\Sigma_{\rm thermal}(T) = \frac{\lambda}{2}\int\frac{d^3k}{(2\pi)^3}\frac{1}{E_k(e^{\beta E_k} - 1)}$

Step 3: High-Temperature Limit

For $T \gg m$ (high temperature), the integrand is dominated by momenta $k \sim T \gg m$ , so $E_k \approx k$ . The integral becomes:

$\Sigma_{\rm thermal} \approx \frac{\lambda}{2}\int\frac{d^3k}{(2\pi)^3}\frac{1}{k(e^{k/T} - 1)}$

$= \frac{\lambda}{2}\cdot\frac{4\pi}{(2\pi)^3}\int_0^\infty dk\,\frac{k}{e^{k/T} - 1}$

$= \frac{\lambda}{4\pi^2}\int_0^\infty dk\,\frac{k}{e^{k/T} - 1}$

Substitute $x = k/T$ :

$= \frac{\lambda T^2}{4\pi^2}\int_0^\infty dx\,\frac{x}{e^x - 1}$

The integral $\int_0^\infty dx\, x/(e^x - 1) = \pi^2/6$ (a standard result).

So:

$\boxed{\Sigma_{\rm thermal}(T) \approx \frac{\lambda T^2}{24}}$

Step 4: The Thermal Mass

The physical mass of the scalar at temperature $T$ :

$m^2_{\rm eff}(T) = m^2 + \Sigma_{\rm thermal}(T) = m^2 + \frac{\lambda T^2}{24}$

The scalar acquires a thermal mass proportional to $\sqrt{\lambda}T$ , even if its vacuum mass $m$ is small or zero.

This is a general feature: at finite temperature, fields acquire thermal masses from loop corrections involving thermal particles. For a coupling $g^2$ and temperature $T$ , thermal masses typically scale as $gT$ .

Physical Interpretation

Why does the scalar get heavier at high temperature? Because it’s propagating through a medium of thermal fluctuations. Scattering off these thermal excitations effectively gives the particle an additional self-energy; a thermal mass.

This is analogous to a photon acquiring an effective mass in a plasma (plasmon), or a quark acquiring thermal mass in the quark-gluon plasma.

8. Real-Time Formalism: When Imaginary Time Isn’t Enough

The Matsubara formalism is elegant for equilibrium thermal quantities; pressure, susceptibilities, static correlators. But for dynamical quantities (transport coefficients, nonequilibrium processes, time-dependent correlators), it struggles because imaginary time isn’t physical time.

The Problem

Imaginary-time correlators $\langle T_\tau\{\phi(\tau_1)\phi(\tau_2)\}\rangle$ are defined on the Euclidean torus $\tau \in [0, \beta]$ . They don’t directly give physical (real-time) observables.

To get a physical retarded or advanced Green’s function, you need to analytically continue $\omega_n \to \omega \pm i\epsilon$ . This works for simple quantities but becomes cumbersome for complex nonequilibrium problems.

Schwinger-Keldysh Formalism

The real-time (Schwinger-Keldysh, or closed-time-path) formalism works directly with real time. The key idea: to compute thermal averages of time-dependent quantities, integrate over a contour in complex time that goes forward from $-\infty$ to $+\infty$ , then back from $+\infty$ to $-\infty$ , and finally down to $-\infty - i\beta$ to close the contour.

This produces a $2\times 2$ matrix structure for propagators, with components:

$G^{++}$ : both times on the forward branch (standard Feynman propagator at $T$ )
$G^{--}$ : both times on the backward branch (anti-Feynman propagator)
$G^{+-}, G^{-+}$ : one time forward, one backward (Wightman functions)

When You Need Real-Time

Real-time formalism is essential for:

Transport coefficients (viscosity, conductivity, diffusion)
Spectral functions ( $\rho(\omega) = -2\text{Im}\, G^R(\omega)$ )
Nonequilibrium dynamics (quenches, thermalization)
Particle production rates from time-dependent backgrounds

For equilibrium static quantities (pressure, condensates), Matsubara is simpler and more efficient.

Kubo Formulas

In linear response theory, transport coefficients are given by Kubo formulas; thermal averages of retarded commutators, which are naturally computed in real-time formalism. For example, shear viscosity:

$\eta = \frac{1}{20}\lim_{\omega\to 0}\frac{1}{\omega}\int dt\, e^{i\omega t}\int d^3x\,\langle[T_{ij}(t, \vec x), T^{ij}(0)]\rangle_\beta$

Computing $\eta$ for the quark-gluon plasma requires the full real-time machinery. The result; that the QGP is a “nearly perfect fluid” with $\eta/s \sim 1/(4\pi)$ ; is one of the most striking results of thermal field theory.

9. Hard Thermal Loops and Effective Thermal Masses

At high temperatures, thermal loops generate new effective interactions and modify particle dispersion relations significantly. The organizing framework is the hard thermal loop (HTL) resummation.

The Scale Hierarchy

In a hot gauge theory like QED or QCD at $T \gg m$ , there are three scales:

Hard scale: $T$ (typical momentum of thermal particles)
Soft scale: $gT$ (where $g$ is the coupling)
Ultrasoft scale: $g^2 T$ (smaller corrections)

At weak coupling ( $g \ll 1$ ): $T \gg gT \gg g^2 T$ . At strong coupling (QCD at $T \sim \Lambda_{\rm QCD}$ ): all scales are similar, and perturbation theory breaks down.

The HTL Effective Theory

For momenta much smaller than $T$ (soft scale), the dominant loop contributions come from hard momenta $\sim T$ flowing in the loop. These can be computed systematically and resummed into an effective action for soft modes.

The leading HTL contributions modify gauge boson propagators:

$D^{\mu\nu}(k) \to D^{\mu\nu}_{\rm HTL}(k; T)$

where the HTL propagator includes thermal corrections from integrating over hard thermal particles.

Thermal Gluon Mass

For QCD at high temperature, the thermal gluon mass squared is:

$m_g^2 = \frac{g^2 T^2}{6}(C_A + \frac{1}{2}n_f) = \frac{g^2 T^2}{6}(3 + \frac{1}{2}n_f)$

for $SU(3)$ . With $n_f = 3$ active light flavors: $m_g^2 \approx 0.75 g^2 T^2$ .

At the LHC (QGP temperatures $T \sim 400$ MeV with $\alpha_s \sim 0.3$ ): $m_g \sim 0.7$ GeV. The gluons in the QGP are far from massless!

Thermal Fermion Mass

Similarly, fermions acquire thermal masses from their interactions with gauge bosons:

$m_f^2 = \frac{g^2 T^2}{4}\cdot C_F$

For quarks in QCD with $C_F = 4/3$ : $m_q^2 \approx g^2 T^2/3$ .

Consequences

The thermal masses modify:

Dispersion relations: $E^2(k) = k^2 + m_{\rm thermal}^2$
Debye screening: static electric fields are screened on length scale $1/m_D$ where $m_D$ is the Debye mass
Magnetic screening: magnetic fields at ultra-soft scales need non-perturbative treatment
Transport coefficients: computed using HTL-resummed propagators

The HTL framework makes thermal QFT calculations at high temperature systematic, analogous to how the standard perturbation theory works at $T = 0$ .

The Infrared Catastrophe

One crucial feature: magnetic gluons in high- $T$ QCD are not screened at any perturbative order. Their propagators remain infrared-divergent. This is the magnetic mass problem; at scale $g^2 T$ , perturbation theory breaks down entirely, and lattice gauge theory must take over.

This is one of the deepest unsolved technical problems in thermal QCD: there’s no known systematic expansion for magnetic sector quantities at weak coupling.

10. Symmetry Restoration at High Temperature

One of the most dramatic predictions of thermal QFT: broken symmetries can be restored at high temperature.

The Physical Picture

At $T = 0$ , a theory might have a spontaneously broken symmetry; e.g., the Higgs vacuum expectation value $\langle\phi\rangle = v \neq 0$ breaks $SU(2)_L \times U(1)_Y \to U(1)_{\rm EM}$ .

At $T > 0$ , thermal fluctuations can destabilize this broken vacuum. If $T$ is high enough, the effective potential is minimized at $\langle\phi\rangle = 0$ , and the symmetry is restored.

This is analogous to a ferromagnet: at low $T$ , it has a magnetization that breaks rotational symmetry. At $T$ above the Curie temperature, thermal fluctuations randomize the spins and restore symmetry.

The Effective Potential at Finite T

For a scalar with potential $V(\phi) = -\tfrac{1}{2}m^2\phi^2 + \tfrac{\lambda}{4!}\phi^4$ (Mexican hat, $m^2 > 0$ ), the one-loop thermal effective potential is:

$V_{\rm eff}(\phi, T) = V(\phi) + V_{\rm thermal}(\phi, T)$

Using the result from section 7 extended to a background field:

$V_{\rm thermal}(\phi, T) = \frac{T^4}{2\pi^2}\int_0^\infty dk\, k^2\ln[1 - e^{-\sqrt{k^2 + m^2_{\rm eff}(\phi)}/T}]$

where $m^2_{\rm eff}(\phi) = -m^2 + \frac{\lambda}{2}\phi^2$ is the effective mass in the presence of the background.

High-Temperature Expansion

For $T \gg m_{\rm eff}$ :

$V_{\rm thermal}(\phi, T) \approx -\frac{\pi^2 T^4}{90} + \frac{T^2 m^2_{\rm eff}(\phi)}{24} - \frac{T m^3_{\rm eff}(\phi)}{12\pi} + \ldots$

Substituting $m^2_{\rm eff}(\phi) = -m^2 + (\lambda/2)\phi^2$ :

$V_{\rm eff}(\phi, T) \supset -\tfrac{1}{2}m^2\phi^2 + \tfrac{\lambda}{4!}\phi^4 + \frac{T^2}{24}\left(-m^2 + \frac{\lambda}{2}\phi^2\right)$

$= -\tfrac{1}{2}\left(m^2 - \frac{\lambda T^2}{24}\right)... \text{wait let me redo}$

Let me combine the $\phi^2$ terms:

$-\tfrac{1}{2}m^2\phi^2 + \frac{\lambda T^2}{48}\phi^2 = \tfrac{1}{2}\left(\frac{\lambda T^2}{24} - m^2\right)\phi^2$

So the effective mass-squared of the scalar at temperature $T$ :

$m^2_{\rm eff}(T) = m^2 - \frac{\lambda T^2}{12}$

Wait, I need to be careful about signs. The original $V$ has $-\tfrac{1}{2}m^2\phi^2$ with $m^2 > 0$ (symmetry broken). The thermal correction adds $+\tfrac{\lambda T^2}{48}\phi^2$ … which is positive.

So the coefficient of $\phi^2/2$ becomes $-m^2 + \lambda T^2/24$ . At $T = 0$ : negative (broken). At high $T$ : positive. The crossover (symmetry restoration) happens at:

$T_c^2 = \frac{24 m^2}{\lambda}$

Above $T_c$ : symmetric phase. Below $T_c$ : broken phase. This is a finite-temperature phase transition.

Example: The Electroweak Transition

In the Standard Model, the Higgs has $|\mu|^2 \sim (80 \text{ GeV})^2$ (from $\mu^2 = \lambda v^2$ with $v = 246$ GeV, $\lambda \approx 0.13$ ). The electroweak transition temperature is:

$T_{\rm EW} \sim 100 - 160 \text{ GeV}$

depending on the details (which include Yukawa contributions, gauge boson contributions, and higher-loop corrections).

Above $T_{\rm EW}$ : $SU(2)_L \times U(1)_Y$ is unbroken. All gauge bosons are massless. All fermions are massless. The universe is radiation-dominated.

Below $T_{\rm EW}$ : the Higgs gets its VEV. W and Z acquire masses. Fermions acquire masses. The electroweak interaction “freezes” into its current form.

This transition happened in the early universe at $t \sim 10^{-10}$ seconds after the Big Bang, when $T$ dropped through the electroweak scale.

First-Order vs. Crossover

First-order transition: two coexisting phases at $T_c$ , separated by a latent heat. Bubbles of the new phase nucleate in the old phase and grow.

Second-order (continuous) transition: smooth change in the order parameter, diverging correlation length.

Crossover: no real phase transition, just a smooth change in thermodynamic quantities.

In the minimal Standard Model, the electroweak transition is a crossover (confirmed by lattice simulations); not a real phase transition. This is important because many baryogenesis scenarios require a first-order electroweak transition for out-of-equilibrium conditions. The observed crossover rules out simple baryogenesis in the SM.

Extensions of the SM (with additional Higgs fields, supersymmetry, or other new physics) can make the transition first-order. This is one motivation for new physics searches.

11. Application: The Quark-Gluon Plasma

What Is the QGP?

At temperatures $T > \Lambda_{\rm QCD} \sim 200$ MeV, QCD undergoes a “deconfinement transition.” Quarks and gluons, normally confined into hadrons, become effective degrees of freedom in a hot, dense plasma; the quark-gluon plasma (QGP).

The QGP existed in the early universe for the first $\sim 10^{-5}$ seconds, between the electroweak transition and the QCD confinement transition. It’s recreated experimentally at RHIC (2000-present) and LHC (2010-present) in heavy-ion collisions.

Thermodynamic Predictions

At $T \gg \Lambda_{\rm QCD}$ , the QGP is approximately a free gas of quarks and gluons. The pressure:

$P_{\rm free}(T) = \left[g_g + \frac{7}{8}g_q\right]\frac{\pi^2 T^4}{90}$

where $g_g = 2 \cdot (N^2 - 1) = 16$ (gluon spin × color) and $g_q = 2 \cdot 2 \cdot N \cdot n_f = 12n_f$ (quark + antiquark, spin, color, flavors).

For $N = 3$ , $n_f = 3$ (active light flavors at QGP temperatures):

$P_{\rm free} = [16 + 21/8\cdot 12/3]\cdot\ldots$

Let me redo this more carefully. The Stefan-Boltzmann pressure for QCD:

$P_{SB}/T^4 = \frac{\pi^2}{90}\left[2(N^2 - 1) + \frac{7}{8}\cdot 2\cdot 2N n_f\right] = \frac{\pi^2}{90}\left[2\cdot 8 + \frac{7}{8}\cdot 12 n_f\right]$

For $n_f = 3$ : $P/T^4 = (\pi^2/90)[16 + 31.5] \approx 5.2$ .

Lattice QCD calculations show that the QGP pressure at $T \sim 2-3 T_c$ is about 80% of this free-gas value. The remaining 20% comes from strong-coupling effects; non-perturbative contributions from residual interactions.

The “Nearly Perfect Fluid”

One of the most striking results from RHIC and LHC: the QGP has very small viscosity. The ratio of shear viscosity to entropy density:

$\eta/s \sim 1/(4\pi) \approx 0.08$

This is close to the conjectured lower bound $\eta/s \geq 1/(4\pi)$ from holographic (AdS/CFT) arguments. Real fluids typically have $\eta/s \sim 10^{-1}$ to $10^3$ , so the QGP is remarkably close to a “perfect fluid.”

Jet Quenching

High-energy partons (quarks or gluons) traversing the QGP lose energy through interactions with the medium. This jet quenching shows up as suppressed jet yields in heavy-ion collisions compared to proton-proton.

The energy loss is calculable in thermal QCD using HTL-resummed propagators. The coefficient $\hat q$ (transverse momentum broadening per unit length) has been measured: $\hat q \sim 1-2 \text{ GeV}^2/\text{fm}$ in LHC conditions.

Transport Coefficients

Calculating transport coefficients in the QGP requires:

Real-time formalism (Kubo formulas)
HTL resummation for soft scales
Non-perturbative input (lattice QCD) for strong-coupling corrections

Modern work combines all three. The state of the art is:

$\eta/s \approx (1-2)/(4\pi)$ (near-perfect fluid)
Bulk viscosity small
Electrical conductivity $\sigma \sim 0.2 T$

The QGP Phase Diagram

The full QCD phase diagram in the $(T, \mu_B)$ plane (temperature vs. baryon chemical potential):

Low $T$ , low $\mu_B$ : hadron gas
High $T$ , low $\mu_B$ : QGP (crossover from hadrons)
Low $T$ , high $\mu_B$ : dense nuclear matter, neutron stars, possibly color superconductors
First-order critical line: somewhere at moderate $T$ and high $\mu_B$ , not yet experimentally established

Finding the critical endpoint where this line terminates in a second-order phase transition is a major goal of ongoing heavy-ion experiments.

12. Application: The Early Universe

The Thermal History

In the first moments after the Big Bang, the universe was hot and dense; a thermal QFT state. As it expanded and cooled, it went through various transitions:

$T > 10^{19}$ GeV: Planck epoch. Unknown physics (quantum gravity).

$T \sim 10^{16}$ GeV: GUT scale. If GUTs are correct, $SU(5)$ or larger symmetry broken here. Phase transitions could produce topological defects.

$T \sim 100$ GeV: Electroweak transition. Higgs VEV develops. $SU(2)_L \times U(1)_Y \to U(1)_{\rm EM}$ .

$T \sim 200$ MeV: QCD confinement. Free quarks and gluons bind into hadrons. Nucleons form.

$T \sim 1$ MeV: Nucleosynthesis. Protons and neutrons fuse into light nuclei.

$T \sim 0.3$ eV: Recombination. Electrons and protons combine into atoms. Photons decouple; this is the CMB.

Every one of these is a thermal QFT process, calculable with the techniques in this document.

Electroweak Baryogenesis (Schematic)

One proposed mechanism for the matter-antimatter asymmetry: during a first-order electroweak phase transition, bubbles of the broken phase nucleate in the symmetric phase and sweep through space. At the bubble walls, CP-violating processes combined with baryon-number-violating sphaleron transitions (non-perturbative in the SM) can generate a net baryon asymmetry.

The trouble: in the Standard Model, the electroweak transition is a crossover (not first-order), and CP violation is too small. So SM baryogenesis fails. This motivates extensions with additional scalar fields or SUSY that provide a first-order transition and more CP violation.

BBN: Big Bang Nucleosynthesis

At $T \sim 1$ MeV, when the universe was about 1 second old, the Hubble expansion rate became comparable to the weak interaction rate. Neutrons and protons stopped interconverting via weak processes like $n \leftrightarrow p + e^- + \bar\nu_e$ . The neutron-to-proton ratio “froze out.”

Subsequently, as $T$ dropped to $\sim 100$ keV, protons and neutrons combined into light nuclei: deuterium, helium-3, helium-4, lithium-7. The observed abundances of these elements (compared to proton) provide one of the most stringent tests of Big Bang cosmology.

BBN calculations require:

Thermal weak interaction rates (computed in thermal QFT at $T \sim 1$ MeV)
Nuclear reaction cross sections
Cosmological expansion (GR Friedmann equation)

All match observations to exquisite precision; one of the cornerstones of standard cosmology.

The CMB

The cosmic microwave background is the thermal radiation left over from decoupling at $T \sim 3000$ K (redshifted today to $T \sim 2.7$ K). It’s nearly perfectly thermal, with small fluctuations ( $\Delta T/T \sim 10^{-5}$ ) that tell us about the early universe.

The CMB provides:

The temperature $T_{\rm CMB} = 2.725$ K
The baryon-to-photon ratio $\eta_b \sim 6\times 10^{-10}$
Evidence for inflation (scale-invariant perturbations)
Constraints on dark matter, dark energy, neutrino masses

Understanding the CMB requires the whole machinery of thermal QFT plus general relativity plus particle physics.

Relic Abundances

Many dark matter candidates are predicted to be thermal relics: particles that were in equilibrium with the SM plasma in the early universe, then decoupled when their annihilation rate fell below the Hubble rate.

The relic abundance is given by a Boltzmann equation:

$\dot n + 3Hn = -\langle\sigma v\rangle(n^2 - n_{\rm eq}^2)$

where $n$ is the number density, $H$ is the Hubble rate, and $\langle\sigma v\rangle$ is the thermally-averaged annihilation cross section.

For a WIMP (weakly-interacting massive particle), the observed dark matter abundance requires $\langle\sigma v\rangle \sim 3\times 10^{-26}$ cm³/s; remarkably close to the weak-scale cross section. This “WIMP miracle” motivates much of the dark matter search.

13. The Sign Problem, Revisited

Thermal QFT brings back the sign problem we encountered in document 10.

At Finite Chemical Potential

Consider QCD with a chemical potential $\mu_B$ for baryon number (relevant for neutron stars, heavy-ion collisions at low energies). The partition function:

$Z = \text{tr}\, e^{-\beta(\hat H - \mu_B \hat N_B)}$

After integrating out fermions, the resulting action for gauge fields has a complex (not real) weight $\det M(\mu_B)$ . Monte Carlo methods fail.

This is the infamous sign problem at finite density. It blocks lattice QCD from systematically studying:

Neutron star equations of state
The QCD critical endpoint at finite $\mu_B$
Quark matter at low temperatures and high densities

Possible Escapes

Various approaches are being pursued:

Complex Langevin dynamics: complexify the fields and evolve stochastically
Lefschetz thimbles: deform the contour of integration
Taylor expansion in $\mu_B/T$ : works for small $\mu_B$
Tensor network methods: alternative to Monte Carlo
Quantum simulation: use quantum computers

None have fully solved the problem. The sign problem remains one of the major technical obstacles in thermal QFT.

Real-Time Sign Problem

A similar issue arises in real-time nonequilibrium calculations; the Schwinger-Keldysh measure isn’t positive-definite, making Monte Carlo sampling difficult.

This limits our ability to compute real-time transport coefficients from first principles. State-of-the-art work combines lattice methods for equilibrium quantities with analytical continuation and uncertainty quantification to extract spectral functions.

14. Appendix: Thermal Formulas Reference

Key Formulas

Matsubara frequencies:

Bosons: $\omega_n = 2\pi n T$ , $n \in \mathbb{Z}$
Fermions: $\omega_n = (2n + 1)\pi T$ , $n \in \mathbb{Z}$

Bose-Einstein distribution: $n_B(E) = \frac{1}{e^{\beta E} - 1}$

Fermi-Dirac distribution: $n_F(E) = \frac{1}{e^{\beta E} + 1}$

Key sum identities: $T\sum_n\frac{1}{\omega_n^2 + E^2} = \frac{1}{2E}\coth(\beta E/2) \text{ (bosons)}$

$T\sum_n\frac{1}{\omega_n^2 + E^2} = \frac{1}{2E}\tanh(\beta E/2) \text{ (fermions)}$

Decomposition: $\coth(\beta E/2) = 1 + 2n_B(E)$ $\tanh(\beta E/2) = 1 - 2n_F(E)$

Useful Integrals

$\int_0^\infty\frac{x}{e^x - 1}\, dx = \frac{\pi^2}{6}$

$\int_0^\infty\frac{x^3}{e^x - 1}\, dx = \frac{\pi^4}{15}$

$\int_0^\infty\frac{x}{e^x + 1}\, dx = \frac{\pi^2}{12}$

$\int_0^\infty\frac{x^3}{e^x + 1}\, dx = \frac{7\pi^4}{120}$

Stefan-Boltzmann Pressure

For free relativistic bosons of mass 0:

$P_{\rm boson} = g_B\frac{\pi^2 T^4}{90}$

For free relativistic fermions of mass 0:

$P_{\rm fermion} = g_F\frac{7\pi^2 T^4}{720} = \frac{7}{8}g_F\frac{\pi^2 T^4}{90}$

Total energy density $\rho = 3P$ (for relativistic degrees of freedom).

Thermal Masses (Order $gT$ )

Scalar in $\phi^4$ : $m^2_{\rm th} = \lambda T^2/24$

Fermion in gauge theory: $m^2_{\rm th} = g^2 T^2 C_F/4$

Gauge boson (Debye, SU(N) with $n_f$ fermions): $m_D^2 = g^2 T^2(N/3 + n_f/6)$

Characteristic Scales

Hard scale: $T$
Debye (electric) screening: $m_D \sim gT$
Magnetic screening: $g^2 T$ (non-perturbative)
Critical endpoint (QCD, if it exists): $T_{\rm cep} \sim 150$ MeV, $\mu_B \sim 600$ MeV (estimates vary)

Problems

Derive the Bose-Einstein and Fermi-Dirac distributions from the appropriate Matsubara sum identities.
Compute the thermal scalar mass from the tadpole diagram at one loop, and verify $m^2_{\rm th} = \lambda T^2/24$ .
Using the high-temperature expansion of the effective potential, derive the critical temperature $T_c^2 = 24m^2/\lambda$ for $\phi^4$ theory with $m^2 > 0$ (symmetry restoration).
For QED at high temperature, compute the photon Debye mass. Show that longitudinal photons acquire a thermal mass while transverse ones do not (at leading order).
Compute the free energy of a gas of massless scalar particles at temperature $T$ by evaluating the logarithm of the partition function.
Estimate the fraction of the universe’s energy density in electroweak gauge bosons and fermions at $T = 100$ GeV.

Closing Note

Finite-temperature QFT extends the framework we’ve built to the hot systems where most of interesting physics actually happens. The key ideas:

Matsubara formalism: imaginary time is compact with period $\beta = 1/T$ , frequencies are discrete
Propagators acquire thermal pieces proportional to $n_B$ or $n_F$
Everything factorizes into vacuum + thermal contributions
Loop integrals are Matsubara sums plus 3-momentum integrals
Thermal masses emerge naturally and scale as $gT$
Phase transitions can be studied by computing the effective potential at finite $T$
HTL resummation handles the hierarchy of scales ( $T \gg gT \gg g^2 T$ )
Real-time formalism is needed for dynamical/transport quantities

Applications

Quark-gluon plasma physics: thermodynamics, transport, jet quenching
Early universe: phase transitions, baryogenesis, nucleosynthesis, CMB
Condensed matter: superconductors, superfluids, quantum critical points
Neutron stars: equation of state at high density and moderate temperature

Connections to Other Topics

This document connects to:

Statistical mechanics: thermal partition function as Euclidean path integral
Condensed matter: same formalism for electron systems, phonons, etc.
Cosmology: thermal history of the universe
Nuclear physics: QGP, neutron stars
Lattice QCD: numerical methods for thermal QFT

The Matsubara formalism is one of the great unifying concepts in modern physics; a single mathematical framework that describes everything from cuprate superconductors to the quark-gluon plasma to the first second of the universe.

Where to Go Next

Related directions for further exploration:

Option B (EFTs in depth): HTL is an EFT, and the formalism you’ve learned here generalizes
Option C (anomalies in depth): anomalies have finite-temperature avatars (anomaly-induced transport, chiral magnetic effect)
Option D (non-perturbative QFT): instantons, solitons, and finite-T critical phenomena
Back to the main sequence: the QFT reference documents are the foundation for all of this

You now have the core framework of thermal field theory. Everything from the phase diagram of QCD to baryogenesis is built on these tools. The physics frontier is wide open.