Thermal Field Theory in real-time formalism: concepts and applications for particle decays

This review represents a detailed and comprehensive discussion of the Thermal Field Theory (TFT) concepts and key results in Yukawa-type theories. We start with a general pedagogical introduction into the TFT in the imaginary- and real-time formulation. As phenomenologically relevant implications, we present a compendium of thermal decay rates for several typical reactions calculated within the framework of the real-time formalism and compared to the imaginary-time results found in the literature. Processes considered here are those of a neutral (pseudo)scalar decaying into two distinct (pseudo)scalars or into a fermion-antifermion pair. These processes are extended from earlier works to include chemical potentials and distinct species in the final state. In addition, a (pseudo)scalar emission off a fermion line is also discussed. These results demonstrate the importance of thermal effects in particle decay observables relevant in many phenomenological applications in systems at high temperatures and densities.


I. INTRODUCTION
The standard theoretical treatment of fundamental fields and their interactions primarily uses the language of quantum field theory (QFT), see for example [1][2][3]. This treatment has historically been extremely successful in describing the behaviour of all the known particles that constitute the Standard Model of particle physics. However, the theoretical framework of this zero-temperature theory does not naturally incorporate effects of the medium relevant at high temperatures and densities [4].
Currently, there is a limited understanding of particle interactions in hot and dense systems despite an enormous work done by several authors and intended to develop a universal formalism for treating the medium-induced phenomena both in and out of equilibrium. The problem of thermal dynamics has engaged physicists for a long time, and early attempts of analysis were made by, for example, Bloch [5]. A pioneer in the realm of thermal (quantum) field theory (TFT) was Matsubara [6] who developed the imaginary-time, or Matsubara, formalism that describes equilibrated systems. This treatment has great resemblance to zero-temperature QFT in the form of the propagators, the diagrammatic structure of the perturbative expansion, and the self-energies. It differs, however, in the way of how time is treated considering it as a purely imaginary quantity. Shortly after the developments by Matsubara, Kubo [7] and Martin and Schwinger [8] provided an important relation between propagators, the so-called Kubo-Martin-Schwinger (KMS) condition, which holds for thermal propagators in equilibrium. Further important developments of thermal theories with equilibrated media, such as an explicit consideration of real-time evolution, were made by Keldysh [9]. As a framework for a first principle analysis, thermo field dynamics was largely developed by Matsumoto et al. [10] while a detailed investigation of QFTs at finite temperatures for real times was made by Niemi and Semenoff [11]. Comprehensive outlines of the path-integral treatment of quantum fields in a thermal medium may be found in [11][12][13], and much of the theoretical concepts presented in this review draws upon that work. Three years before the work of Kobes and Kowalski [12], Weldon [14] provided self-energy calculations performed using the Matsubara formalism and presented a condensed and clear overview and interpretation of the thermal decay rates. Discontinuities of Green's functions, important for the decay rates have also been studied in detail by Kobes and Semenoff [15,16]. Further studies of thermal N -point function, important for higher order processes have been performed by e.g. Evans [17] This review aims to give a detailed and comprehensive overview of basic concepts and implications of finite-temperature effects, and thereby providing an insight on the importance of implementing the TFT framework for a consistent and universal description of high to low temperature and density phenomena. In particular, we present an explicit calculation of typical thermal effects on observable decay rates for several generic processes found in Yukawa-type field theories. In the following Sec. II, an overview of the underlying theory is presented starting from the statistical approach introduced in QFT by Wagner [18]. This general approach is valid in as well as out of equilibrium. Then, the equilibrium propagator is presented, followed by an overview of the Matsubara formalism, Sec. III, and the real-time formulation, Sec. IV. Having the real-time propagators and self-energies established, the concept and definition of the thermal decay rate is presented in Sec. VI. In Secs. VII-X, thermal decay rates are presented in the real-time formalism for several processes. Here, we present such examples as a (pseudo)scalar particle decaying into two distinct (pseudo)scalars, Sec. VII, a (pseudo)scalar decaying into a fermion-antifermion pair, Secs. VIII-IX, and finally in Sec. X, a (pseudo)scalar emission off a fermion line is discussed. In Sec. XI, a short summary is given.

II. ELEMENTS OF THERMAL QUANTUM FIELD THEORY
This section provides the theoretical formulation of TFT with the purpose of introducing the reader to the relevant methodology that will then further be used in explicit calculations of thermal decay rates. In order to obtain a diagrammatic formulation of TFT, the Green's function techniques are deployed in terms of a thermodynamic extension of the Wick's theorem. This technique is then used for analysis for two-body decay rates in the thermal medium that are related to the corresponding self-energies.

A. Statistical treatment of initial states
The underlying principle of TFT will be outlined in this subsection. The discussion is mainly based on the broad and comprehensive article by Wagner [18] and the initial formulas are valid in equilibrium as well as out of equilibrium for arbitrary initial distributions. This general treatment precedes a discussion of equilibrium theories, e.g. the Matsubara formalism or the real-time formalism.

Operator expectation values
Experimentally measurable observables in QFT are usually expressed as expectation values of certain operators evaluated at any given time. An arbitrary measurement is defined by two characteristic properties that can be introduced as follows.
The first property is the preparation of the initial state which fully specifies the system at some initial time t in . It is reasonable to assume that no initial state can be determined exactly; rather, one should adopt the view that the initial state can be prepared up to a probability distribution ρ over pure states so that the initial preparation of any system results in a mixed quantum state. Mixed states are statistical ensembles of pure states where the ensemble specifies some lack of knowledge of the system e.g. due to noise, entanglement with larger systems etc. Using the language of statistical quantum mechanics, the distribution ρ assigns a weight ρ(n) ∈ [0, 1] to the pure states |n such that the weight describes a fraction of the ensemble in each state. The weights are normalised according to n ρ(n) = 1 and, in the Fock space spanned by |n , they provide a definition of the density matrix :ρ (t in ) = n |n(t in ) ρ(n) n(t in )| . (II.1) This expression describes a classical probability distribution over the pure states. Note that any state, with pure or mixed initial preparation, may be described by this formalism 1 . Using statistical mechanics, the operator expectation value (the observable) is given by = n n(t in )|ρ(n)Ô(t in )|n(t in ) (II. 2) at the initial time t in for an observable O. Hence, the initial preparation of a system is equivalent to the determination ofρ(t in ).
With the initial state now dealt with, the second property to consider is the time evolution of the system. This process is specified by the HamiltonianĤ(t), and, in the Heisenberg picture, the system evolves according to the Heisenberg equation of motion: The states are time-independent in this picture and specified fully at t in . A consequence of the above equation of motion is that the time-dependent expectation value, the one-point function, is given by Ô H (t) = n n(t in )|ρ(n)Ô H (t)|n(t in ) . (II.4) The subscript H will be dropped from now on.

Expansion of the density operator
The initial density operator for an arbitrary experiment may be a multi-particle operator. The most general proof of Wick's theorem was provided by Danielewicz [19]. This work implemented the necessary condition that the expectation value must be taken over non-interacting operators with respect to a one-particle density operator in order for the theorem to hold exactly. However, Wagner [18] provided an appropriate expansion of correlation functions that allows the use of Wick's theorem. This expansion is two-fold and considers both an arbitrary (potentially multi-particle) initial density operator and the time evolution operator. Let us discuss such an expansion in more detail.
A generic density operator may always be expressed in an exponential form due to the physical requirements of hermiticity and positivity. This form is realised by the introduction of an operator B defined at the initial time t in . This is done in complete analogy to the formalism developed by Matsubara [6] so thatρ which is the starting point of a theory that describes thermal systems. Note that for the grand canonical ensembleB =Ĥ − µN and λ = β, with µ being the chemical potential,N -the particle number operator, and β = T −1 -inverse temperature, the formalism of Matsubara is recovered. The exponential form ofρ closely resembles the Boltzmann distribution which emerges in any system in thermal equilibrium. One may further define the one-particle density operatorρ 0 as followŝ Here, the general one-particle operatorB 0 is extracted fromB leaving the residual operatorB : B ≡B −B 0 . SinceB 0 is an analogue to the free one-particle Hamiltonian of the zero-temperature theory, a generalised interaction picture can be defined with respect to the time-independent operator in the Schrödinger picture usingB 0 so that for some parameter τ . Comparing the translation operators (exponentials) on the right-hand side of Eq. (II.7) to that of Eq. (II. 6), one arrives at a connection between τ and λ shown below. Making use of this generalised interaction picture to express the residual ofB after the extraction ofB 0 , Wagner [18] defined an expansion of the multi-particleρ in terms of the one-particleρ 0 . From Eq. (II.5), we haveρ Here, the parameter τ is purely imaginary for real λ and it defines an imaginary contour of evolution (or integration, in what follows) C ν in the complex plane, see Fig. 1. The new operatorŜ Cν , a generalisation of the time evolution operator of the zero-temperature theory, was introduced above purely by trivial extraction and takes the form The exponential on the right-hand side has been contour-ordered along the vertical contour segment that goes from 0 to −iλ according to T CνÔ1 (τ 1 )Ô 2 (τ 2 ) = Ô 1 (τ 1 )Ô 2 (τ 2 ) τ 1 ≥ τ 2 on C ν , ηÔ 2 (τ 2 )Ô 1 (τ 1 ) τ 1 < τ 2 on C ν . (II. 12) The sign expressed by η = +1 (−1) accounts for the case of commuting (anticommuting) operators. As mentioned above,B has been expressed in the generalised interaction picture with respect tô B 0 . Wagner [18] thereby provides an arbitrary expectation value at t in as follows Ô (t in ) = Z 0 Z n0 n 0 (t in )|ρ 0ŜCνÔ (t in )|n 0 (t in ) .
(II. 13) This has the same form as Eq. (II.4) if one regards Z 0 /Z Ŝ CνÔ (t in ) as a non-interacting operator (i.e. expressed in the generalised interaction picture) to be averaged with respect toρ 0 . Notably,ρ 0 is a one-particle density operator defined at the initial time t in . As a consequence, Wick's theorem holds for each term in the expansion ofŜ Cν following the criteria presented by Danielewicz [19]. The emerging contour of integration is the vertical line in Fig. 1. The formalism outlined above is valid for any system that is defined by a distribution in the exponential (thermal) form Eq. (II.5) at some initial point. Quite generally, the formalism has been shown [18] to hold for a wide family of systems since the multi-particle initial distribution may always, at least formally, be expressed in terms of quantities that obey the Wick's theorem.

Real-time evolution operator
The generalised time evolution operator, formally defined by Eq. (II.11), transforms an operator into its average along the contour of Fig. 1. In this section, the definition of this operator will be extended in order to incorporate evolution along the real-time axis thereby allowing for physical time dependence. This development is necessary in order to work out the so-called real-time formalism in which real-time quantities may be extracted directly.
Wagner [18] and Danielewicz [19] made an observation that the time evolution operator must be expressed as a non-interacting operator in order for the Wick's theorem to hold. However, the general time-dependent Hamiltonian contains multi-particle interactions. This issue was resolved in the previous section in Eq. (II.13) for the vertical contour in Fig. 1 since this expression provides an expansion in which the Wick's theorem holds for each term inŜ Cν . In general, in order to compute observables one is interested in the resulting n-point correlation functions. The latter can be generalised in the TFT as follows Ô (t 1 , t 2 , . . . , t n ) = Ô 1 (t 1 )Ô 2 (t 2 ) · · ·Ô n (t n ) , (II. 14) and then extended further to directly incorporate the real-time arguments as opposed to the purely imaginary temporal dependence along the C ν contour discussed in the previous section. Keldysh [9] developed an expansion of the time evolution operator that is valid on the real axis. As shown in Fig. 2, a contour segment C 1 that goes up to some largest time t fi on the real axis is added. A second segment C 2 goes back to t in , again along the real axis, before the piece C 3 goes down to t in − iλ. The segments on top of the real axis go through all temporal arguments of the n-point function of interest. Note that the apparent offset from the real axis in the figure is purely for display purposes: both contour pieces C 1 and C 2 lie exactly on top of the axis. The contour ordering T C may be defined along the extended Keldysh contour. In complete analogy to the previous section, the notion ofB may be extended by the introduction of The new operatorK(t) is defined as the real-time Hamiltonian of a given theory for real-time arguments and the generalised Hamiltonian discussed in the previous section on the vertical contour piece. The procedure may be followed through introducing the one-particleK 0 and residualK operators that are defined over the entire contour, similarly toK (see the analogous definition of B ,B 0 following Eq. (II.6)).
The two-point function may be written as The contour is extended along the real axis going through the time arguments of any n-point correlation function. The times t in and t fi are arbitrary and may be suitably chosen for a given system under the condition that all time arguments of the observable are incorporated. Note that the apparent offset from the real axis is superficial; C 1 and C 2 lie exactly on top of each other and the radii of the arcs connecting the contour segments are limited to zero.
The subscripts of the brackets here indicate that the trace has to be taken with respect to the one-particle density operatorρ 0 defined in terms ofK 0 . Besides, T C orders the operators along the entire contour C and the time evolution operator is expressed in the generalised interaction picture with respect to K 0 analogously to Sec. II A 2. Hence, a formulation of an arbitrary system, valid also out of equilibrium, exists in terms of the series expansion ofŜ C . For such series, the thermodynamic version of Wick's theorem holds [19]. Extracting the physical quantities from such an expansion is unfortunately anything but a straightforward procedure since it is clear that the series contains terms other than real-time correlation functions as well. Several authors have analysed the full expansion out of equilibrium, see for example [9,18,[20][21][22]. In particular, Wagner [18] provides a systematic series expansion of the contour-ordered two-point correlation function for fermion fields and sets up the framework in which one could attempt to find a solution to the generalised Dyson's equation. From this point on, this work will focus on approaches to thermal equilibrium theories, mainly developing their imaginary-and real-time formulations, with the main goal of presenting the decay rates for processes involving neutral scalars and Dirac fermions as practically relevant example calculations in the TFT.

B. Field operators in thermal theory
As was elaborated above in Sec. II A 2, the statistical approach to TFT produces a complex-time dependence through the formal equivalence between time and inverse temperature. The expectation values of operators that are needed for the evaluation of the corresponding observables are expressed above as traces over operators with a statistical weight implemented by the density matrix. This subsection provides relevant details on properties of the field operator and its modes in complete analogy to the zero-temperature theory.

KMS condition and ladder operators
The trace interpretation enforce a periodicity condition on the n-point functions of Eq. (II.14). Let us discuss this condition in more detail.
Consider, as an example, the simplest class of n-point functions -the statistical average with respect to the canonical density matrix,ρ = exp [−βĤ], of N scalar field operators expressed in the Heisenberg picture. In the case of a one-particle Hamiltonian, the average procedure results in Using the cyclic property of the trace, it is straightforward to show that It was noted already in 1932 by Bloch [5] that the sandwiched field operator that appears in the above calculation can formally be represented as a result of temporal evolution due to the Heisenberg equation of motion of Eq. (II.20) The condition of Eq. (II.19) is the well-known Kubo-Martin-Schwinger (KMS) condition [7,8] that holds for all thermal n-point functions. Now, any field operatorφ(x) of a free theory, as a function of the space-time variable x, may be expanded in the basis of creation and annihilation (ladder) operators. For a bosonic (fermionic) field, these operators (anti)commute. In the case of a bosonic field operator, the commutation relations of the ladder operators may be expressed according to for modes k, l and the commutator vanishes for any other combination of operators. Introducing the mode functions f k (x), f * k (x), the free field can be expanded aŝ The mode functions of the free massive field are the solutions to the Klein-Gordon equation (∂ 2 + m 2 )f k (x) = 0 which, for spinless bosons in a finite volume V = L 3 , read with ω k = √ k 2 + m 2 denoting the energy of the kth mode, and k = 2πn/L where n = n x n y n z T is the unit 3-vector along the particle momentum. As usual, the dynamics of the field operator above is governed by the free field Hamiltonian with the conjugate momentum π(x) = ∂ 0φ (x) and the free-field potential V = − 1 2 m 2φ (x). In a theory described by the free-field Hamiltonian, the thermal average of any combination of ladder operators may be evaluated. In analogy to the Heisenberg picture, one may define the shifted operator asâ † where the second equality is extracted from the commutator Ĥ 0 ,â † k = ω kâ † k . Making use of the KMS condition, one notices that Applying the commutation relation of Eq. (II.21) results in From this, one obtains where n(ω k ) is the thermal distribution function of Bose-Einstein statistics. It is a straightforward exercise to show, in a similar fashion, that Adding the case of anticommuting ladder operators, the thermal distribution generalises to for commuting (η = +1) and anti-communiting (η = −1) operators. The latter guarantees that the fields satisfy the periodicity condition imposed by the trace interpretation.

Thermodynamic Wick's theorem
The most general proof of the Wick's theorem can be found in Ref. [19]. Following the basic set up introduced above, one may deploy the same strategy as used for deriving Eq. (II.28) and examine the general expression It is a simple procedure of applying the commutation relation of Eq. (II.21) together with the KMS condition Eq. (II.26) to show that this expression vanishes if i = j. If the number of creation and annihilation operators is equal, this expression allows for a simple proof of the thermodynamic version of the Wick's theorem. Indeed, applying the KMS condition once one gets Then, moving the exponential factor to the left-hand side and commuting the creation operatorâ † k1 through the series of annihilation operatorsâ lj , one finds Here, the first term is the average of interest and it can be moved to the left-hand side where it combines to gives the thermal function. Collecting all terms proportional to the δ-functions results in Moving the thermal factor to the right-hand side, one identifies the Bose-Einstein distribution such that Eq. (II.33) becomes (II. 35) The rightmost factor in the sum is expanded in a similar fashion using {â mn } Hence, by induction, the full average â † k1 · · ·â † kiâ lj · · ·â l1 0 reduces to for j > i. Note, the rightmost factor in the nested sum must vanish since the average of a product of annihilation operators is zero. Similarly for j < i, Hence, it can be concluded that So, any average over an equal number of creation and annihilation operators may be expanded as a sum over the product of every possible combination of two-operator averages. All other combinations of creation/annihilation operators have a vanishing thermal average. This is the essence of the thermodynamical Wick's theorem.

Thermo field dynamics
In Secs. II B 1-II B 2, an operator formalism for treating dynamics of the fundamental fields has been introduced. As seen already in Sec. II A 3, the physical operators of TFT can be defined on the real-time contour shown in Fig. 2. The theory gets simplified significantly if one considers a thermal system in equilibrium. The study of thermal operators acting on the states of the equilibrated Hilbert space is referred to as thermo fields dynamics. This approach was investigated first by Umezawa et al (see e.g. Refs. [10,23,24]) and the idea results in a doubling of the Hilbert space. Thermo field dynamics has been further developed by Fujimoto et al [25] and incorporates, for example, a thorough structure of the thermal vacuum as well as provides a fundamental understanding of gauge theories in TFT [26].
The emergent property of the doubling of degrees of freedom is present also in the path-integral formulation of the real-time formulation. This makes the evaluation of dynamical thermal quantities, e.g. the full thermal propagator, cumbersome in comparison to the case of the zero-temperature theory. However, both thermo field dynamics and its path-integral version provide the means for evaluating the observables without considering the imaginary times, a property which aids the physical interpretation. We refer the interested reader to several books written on this topic, see for example [23,26].

C. Thermal two-point functions and propagators
Knowing the properties of the thermal trace from Sec. II B, the first non-trivial correlation, the two-point function, may be discussed in some detail. Its explicit form will be obtained and a comparison to the corresponding quantity from the zero-temperature theory is made straightforward. Besides, the thermal propagator of the Klein-Gordon theory is described, to be further generalised to charged anti-commuting multi-component fields in later sections.

Thermal two-point correlation functions
Owing to the discussion in previous sections, the tools have been outlined in order to calculate the two-point correlation function of a free neutral scalar field which will be used as an example for thermal calculations later on. First, let us define the two-point correlation functions ordered in the temporal arguments according to Now, the first one can be computed by expansion of the field in terms of the ladder operators as Only averages over an equal number of creation and annihilation operator will contribute to this expression. Hence, Substituting the expression for the mode functions f k from Eq. (II.23), we obtain In the thermodynamic limit 2 , an integral replaces the sum [13] so that In order to produce an integral in four-space over k, the following identity may be employed so that the correlation function becomes 2 Keeping the density of states constant, the replacement 1 Since ω k = ω −k by definition of ω k , the second integral could be recast by a change of variables according to k → −k. Hence, the exponential could be taken out for k = (k 0 , k). Now, by identifying and by defining the free spectral density the correlation function finally transforms to Note that the free spectral density ρ 0 (k) here should not be confused with the initial density matrix. Similarly, the Fourier transform of C < (x, x ) may be found as for which for following relation holdsC The latter can be straightforwardly proven by utilising the KMS condition.
Let us now consider the average of the commutator Making use of the relation between the two Fourier transforms in Eq. (II.51), one finds that utilising the fact that [1 + n(k 0 )] 1 − e −k0β ≡ 1. Such an average be used in analogy to the zero-temperature theory e.g. in order to define the thermal retarded and advanced propagators as follows where Θ(x) is the standard Heaviside step function.

Thermal propagators
In analogy to the correlation function Ô (t 1 , t 2 , . . . , t n ) of Eq. (II.14), one may define the thermal Green's functions as that are the time-ordered thermal averages of products of the field operators φ (x i ) as prescribed by Eq. (II. 16).
As an important non-trivial example, the two-point thermal Green's function for the free theory reads and in order for G 0 (x, x ) to be a Green's function of the theory, it must obey the Klein-Gordon equation of motion Using the definition of G 0 (x, x ), one identifies the thermal propagator of the free field which should itself be a Green's function of the Klein-Gordon equation, i.e.
In order to demonstrate this explicitly, let the derivative act on the variable x as will be denoted by ∂ 2 (x) while treat x as a constant shift. Then, Using the identity dΘ(t)/dt = δ(t) and rewriting the terms with second-order derivatives as and we notice that the terms proportional to δ(t − t ) disappear. Besides, one observes that due toD > (k) ∝ δ(k 2 − m 2 ), and the same holds true for the term containing (∂ 2 (x) + m 2 )D < (x, x ). Hence, the full expression is straightforwardly evaluated to thus, completing the proof.

Spectral representation
A special form of the thermal propagator D(x − x ) is known as its spectral representation in terms of the spectral density defined above in Eq. (II.48). The propagator introduced in the previous section may be re-written as an integral in momentum space as Note that e −k0β [1 + n(k 0 )] = n(k 0 ), recovering the spectral representation for a bosonic scalar field (see e.g. Ref. [27]) For practical applications in Feynman diagrams computation, it is instructive to turn into momentum space and to compare the Fourier-transform of the thermal propagator to its zero-temperature counterpart. To recover an explicit form of the Fourier transform of D(x − x ), one starts with Eq. (II.68) and writes Using the following representation of the Heaviside function, the integral in Eq. (II.69) may be re-written as where k 0 = k 0 − s 0 . This rather compact expression for the thermal propagator has a straightforward interpretation. The first term contains no reference to the thermal (medium) parameter β and resembles the propagator of the zero-temperature theory. The last term contains the thermal distribution function n(k 0 ) that characterises the medium in equilibrium through the temperature parameter T ≡ β −1 and for the Bose-Einstein (Fermi-Dirac) distribution governing the statistics of bosons (fermions). Note, this function vanishes when T → 0. The thermal propagator has therefore naturally been split into thermal and non-thermal terms.

D. Contour path-integral formulation
From the outline of the TFT given above, together with some illuminating examples, it has become clear that in order to perform explicit calculations, one must specify the initial thermal density. Focusing on the case of density operators that describe the thermal systems in equilibrium, this subsection overviews the path-integral formulation of such systems starting from the density distribution of the grand canonical ensemble. Then, the generating functional is discussed before the perturbative expansion of the thermal path-integral formalism is considered. The general propagator of TFT is discussed and the KMS condition is then restated in the context of path integral formulation.
One-particle quantities will be assumed initially since it has been shown above that the thermal functions of multi-particle systems may be recast into a series expansion over one-particle quantities and much of the theory from this point on is initially introduced for free particles. Interactions are thereafter introduced perturbatively. The material presented here is primarily based on an extensive discussion by Landsman and van Weert in Ref. [13].

Grand canonical ensemble
In thermal equilibrium, the HamiltonianĤ is time-independent, and a system of charged particles is characterised by the density operator of the grand canonical ensemblê Here,Q a are the conserved charges with a ∈ {1, A}, and Φ is the thermodynamical potential related to the partition function through Φ = − log Z. The charge operatorsQ a and the Lagrange multipliers α a may be related to the number density operator and the chemical potential, respectively, asN where V is the volume of the system and µ a corresponds to the chemical potential related to each type of charge a. As seen in the previous section, the form ofρ determines the expectation value of any n-point function 3 .

Generating functional
Analogous to zero-temperature QFT and to the scalar-field theory example in Eq. (II.56), the n-point correlation functions (or thermal Green's functions) of arbitrary (commuting or anticommuting) field operatorsÎ ∞ can be defined as (II.74) An important detail brushed over in the example calculations above is the fact that the thermal Green's functions depend explicitly on the choice of integration contour C. This is evident mainly in Eq. (II.17). The contour of integration may be chosen quite arbitrarily but it becomes restricted for diagrammatic formalisms that require the thermodynamic version of Wick's theorem to hold. The thermal Green's functions can be generated by the functional expressed in terms of the c-number sources j(x) through functional differentiation as (II.76) The normalisation here is The generating functional of Eq. (II.75) contains all n-point functions G C (x 1 , . . . , x n ). In analogy to the zero-temperature theory, all disconnected diagrams can be factorised out by considering only the connected Green's functions G (n) The connected Green's functions satisfy the following the cluster property: which may be verified by direct evaluation. That is, the connected pieces fall off to zero when the distance between the field points increases. The expansion of W [j] is often referred to as the cumulative (or cumulant) expansion in statistical mechanics. The resulting series trivially generates the connected Green's functions, i.e. the sum of all diagrams that are fully connected. As a consequence of the logarithmic relation between Z[j] and W [j], the contribution from vacuum bubble diagrams, Z[0], factorises and may be absorbed by the normalisation of the Green's functions. An illustration of the order-by-order cancellation of disconnected bubble diagrams can be found in A.
Time-ordering procedure on the contour, T C , may be explicitly expressed by parametrising C according to t = z(τ ), with τ being real and monotonically increasing along the contour. Through the definition of the contour step function the two-point function, taken as an example, reads where the sign factor η = ±1 for commuting/anticommuting (bosonic/fermionic) field operators. Weldon [14] showed that this expression is well behaved only if the imaginary component of the contour C does not increase with τ . Feynman, Matthews and Salam [28,29] provided the first proof of this statement in the context of the path-integral formulation. It should be mentioned at this point that the terms scalar, boson and fermion will be used extensively throughout the following theoretical outline below, as well as commuting and anticommuting fields and operators. Some of those labels are used in a somewhat different context by various authors and hereby we would like to clarify on how such terms are used within this work.
The term 'scalar' intends to reference a scalar quantity both in Lorentz and in any internal space. Hence, a 'scalar field' associates a single value to each point in space-time but it may be either a commuting or an anticommuting field (the latter being a Grassmann variable) in this review. A 'boson' refers to an operator that obeys canonical commutation relations. Hence, a 'scalar boson' is the specification of a scalar (i.e. not a vector) field that obeys Bose-Einstein statistics. The use of the term 'fermionic scalar' in this review might induce some confusion for the reader; the label refers to a scalar field that anticommutes (which is not a vector in a given space). This distinction of scalars is introduced simply because it is an easier exercise and a common approach in the literature to present the initial theory in terms of either commuting or anticommuting scalar (neutral one-component) fields. The formalism is then later extended to include fields that carry Lorentz or spin components. Hence, proper care must be taken when constructing the Lagrangian so as to only include field types with the proper structure in Lorentz and spin space respecting the relevant symmetries of the theory. The term 'fermion' is applied along its most common usage to Dirac spinors that carry components in spin space.

Path-integral formulation
The functional that generates the n-point functions was determined in the previous section as a thermal expectation value over operators (Eq. (II.75)). In the path-integral formulation, the state vectors that span the Fock space introduced in Sec. II A can be used as a basis of coherent states that are eigenstates of the annihilation component of the field operator. The discrete modes are then exchanged for the continuous three-space variable 5 , x. This coherent state basis may be chosen to express the field states |φ(x); 0 at time t = 0 that are eigenstates of the field operatorφ(x): (II.82) The theory of neutral bosonic field considered in this subsection provides a simple example with vanishing chemical potential. The coherent states generalised to the continuum picture can be used to express the functional integration measure and further to express the thermal trace in the generating functional as Note that the bra and ket states here must be taken at equal times to uphold the trace interpretation. Given the time evolution of the field operator in Eq. (II.82), the states absorb the exponential through |φ(x); t = e iĤt |φ(x); 0 . (II.84) The action of the canonical density operator inside the trace of Eq. (II.83) can then analogously be interpreted as a complex shift in time, i.e.
This, so far, purely formal equivalence of the thermal density operator to the time-evolution operator was first noted by Bloch [5] already in 1932 and from this observation follows the introduction of the complex temporal contour discussed in Sec. II A. Note specifically that the eigenvalue of the field operator, ϕ, remains in order for the trace interpretation to hold. The functional measure [Dϕ] is such as it may only pick out the fields with this periodicity condition over iβ so that With the formal observation above, the shifted trace can be restated in terms of some arbitrary initial and final times t i , t f with the aid of the Feynman-Salam-Matthews (FSM) formula [28,29], if the path-integral on the right-hand side is taken over all c-number fields that satisfy the boundary conditions In fact, Eq. (II.86) remains valid for arbitrary initial and final times under the single restriction that the imaginary component of the integration contour C connecting t i and t f may not increase along the direction of the contour [13]. Hence, C must go downwards in the complex plane or extend in a parallel direction to the real axis. When specifying the theory to one with no derivative couplings, the normalisation N in Eq. (II.86) may absorb the Gaussian integration over the conjugate momentum so that N → N . The action is then and hence the generating functional becomes in terms of the field φ(x) and its source j(x). Here, t i = t in and t f = t in − iβ are connected via the contour C.
Expressed above as a path-integral in terms of fields rather than operators, the generating functional accounts for the time ordering on the contour. Further, the averaging procedure has been reinterpreted in the path-integral formulation to be taken with respect to the action as a statistical weight. The ill-defined 6 normalisation Z[0] cancels out in the thermal Green's functions of Eq. (II.76). This quantity contains all vacuum bubbles and is briefly discussed in terms of diagrams in Sec. A.
In order to apply the tools of perturbation theory, one first splits the Lagrangian into its free and interaction parts as L = L 0 + L I , such that the prefactor Z[0] in Eq. (II.89) factorises into free, Z 0 [0], and interaction, Z I [0] parts. In order to obtain the Feynman rules and a diagrammatic series, the Z I [0] then undergoes a perturbative expansion into a series over the interaction term i.e.
By observing that the expansion series can be obtained by functional differentiation w.r.t. the source j(x), the argument of S I may be replaced by the functional differential, and the exponentiated interaction term, no longer a functional of φ, can be taken out of the functional integral. The full Z[j] can thus be rewritten in analogy to that in the zero-temperature theory as Note, however, that the series expansion is formally presented in full generality w.r.t. a given contour C and the individual terms in this expansion may only have a clear interpretation once a particular choice for the contour C is specified.

E. Thermal propagators
In this subsection, the propagator of a free neutral boson will be connected to the two-point correlation function. The latter is of great importance for practical computations and it is found by the evaluation of expectation values in the form while making use of the time-ordering procedure on the contour, see Sec. II D 2. This section will further generalise the discussion of thermal propagators to include anticommuting charged multicomponent fields and we follow closely the notation and procedure provided by Landsman and van Weert [13].

Propagator of the free neutral scalar boson
In the special case of a free scalar particle with no derivative couplings 8 , here a neutral (µ = 0) boson with the free Klein-Gordon Lagrangian the free generating functional in Eq. (II.89) may be rewritten by means of a change of variables as follows The shift in φ(x) is chosen carefully in order to complete the square of the free Lagrangian. The emerging thermal bosonic propagator D C (x − x ) is defined in relation to the differential operator of the equation of motion and it satisfies (recall the definition t ≡ z(τ )). Therefore, the propagator is related to the free two-point Green's function as follows Generalising the discussion above in Sec. II C 3, it can also be rewritten in the spectral representation, where n(k 0 ) is the thermal distribution function defined in Eq. (II.32). This result can be found by a generalisation of the expression first derived by Mills in Ref. [30] for the non-relativistic case and then by Dolan and Jackiw in Ref. [27] for ordinary time ordering. Here, the spectral density for the free neutral bosonic field reads Note that this quantity is not related to the distribution found in Sec. II A and will not appear from this point on other than in its explicit exponential form.

Generating functional for an arbitrary free field
Consider a general multi-component covariant complex field (bosonic or fermionic) described by the operatorÎ ∞ i α (x) [13]. Here, α is the index of the Lorentz representation of the field, e.g. a spinor or vector index and i is an index of a representation of an internal symmetry or a field generation index, i.e. not related to space-time transformations. Let the operator transform under some representation of the Lorentz group J αβ [Λ]. The field is considered to be charged under a given set of internal symmetries and the conserved charges are denoted as q ij a so that The chemical potential for the fields charged under multiple symmetries is generally given by the summed quantity µ = a µ a q a ; see the grand canonical ensemble in Sec. II D 1 applied to the eigenstates ofQ a , such that Generalising the plain Klein-Gordon case discussed above, it is reasonable to start with a free Lagrangian in the following quadratic form Here, the conjugated or adjoint field is defined as 9 is the differential operator of the free field in the medium which generally depends on the chemical potential. The free thermal Green's functions that arise from this Lagrangian can be obtained from the generating functional. In the path-integral formulation, this functional may be written in terms of two source fields where j(x) is the source of I ∞ while(x) is the source ofĪ ∞ , and the fields are fully contracted over all discrete indices, e.g. I . The µ-dependence has been shifted away from K to the source terms by means of a change of variables introduced by Weldon [31], under which the fields and the differential operator are transformed as thus treating the field and its adjoint counterpart as independent variables.

Thermal propagator and KMS condition: a generic treatment
The functional integral of Eq. (II.103) is defined in order to generate the thermal Green's functions (n-point functions) under functional differentiation and therefore it must represent the thermal trace of Eq. (II.4): In order to recover the trace interpretation of the generating functional, the functional measure should be taken only over c-number fields I ∞ ,Ī ∞ such that the following periodicity condition is satisfied over a period iβ such that 10 The periodicity leads to the Kubo-Martin-Schwinger condition (KMS) on the propagator which was introduced by Kubo [7], implemented by Martin and Schwinger [8] and is presented in what follows.
The free generating functional of Eq. (II.103) can be rewritten in terms of the general free thermal propagator by a shift of the field and its adjoint analogously to the case of the scalar field in Sec. II E 1. Contracted over all indices, the result is In order to complete the square and extract the propagator D ij C αβ ≡ K −1 ij αβ , the thermal propagator must satisfy the differential equation The generating functional must be invariant under the introduced shifts of the fields. This condition fixes the solution of the differential equation so that the propagator can be related to the two-point function as This is a special case of the general relation Eq. (II.16) with the final expectation value expressed in the path-integral formalism over the fields rather than the operators, analogously to Eq. (II.89). Guided by the definition of the time ordering along the contour T C , the free propagator may be split using the contour step-function of Eq. (II.80) so that Let us compare this expression to Eq. (II.81) for the general n-point function. The imposed boundary condition on the field of Eq. (II.106) enforces a boundary condition on the propagator that relates its advanced and retarded components. As a consequence, invariance of the path-integral over the time shift by iβ holds. This way, we arrive at the KMS condition in the following generic form, which is a proper generalisation of the condition Eq. (II. 19) previously introduced in the case of a neutral scalar field. This condition is important for thermal theories as relates several thermal quantities in the real-time formalism (see Sec. IV). Note that any relativistic, multi-component field I ∞ i α (x) should satisfy the Klein-Gordon equation. Therefore, assuming a multi-mass Klein-Gordon divisor d ij αβ (i∂) to exist according to Ref. [32] and further references therein, we have The m l 's represent the mass spectrum and can be inferred from the structure of K [32,33]. For an extended and comprehensive discussion regarding the mass spectrum, see the lecture notes by Wightman [33]. The solution to Eq. (II.109) simply becomes in terms of the scalar propagator encountered earlier. A slight reservation should be made regarding the form of D C (x − x ). The solution presented in Eq. (II.98) of Sec. II E 1 is that of a neutral scalar boson. When incorporating anticommuting as well as charged fields (µ = 0) this propagator is modified slightly to become with the distribution function The ± corresponds to charged particles/antiparticles through the definition ω ± = β(k 0 ∓ µ). Generally, the spectral density of a free field can be found as (II.117) Here, the discontinuity introduced by Ref. [32] and further discussed in Ref. [13] is defined over the real axis as Disc f Note, the results of this section rely on the fact that the contour step-function commutes with the Klein-Gordon divisor. Further, the charge operators are assumed to commute with L 0 so that the matrices q a commute with both K(i∂) and d(i∂) [13].
Let us now summarise the most important points from the above discussion as follows: • Time integration is to be taken along a contour C that begins at an arbitrary initial time t in and goes down to the final time t in −iβ with the restriction that the imaginary component of the contour can not increase 11 .
• The propagator D C (x − x ) depends explicitly on the choice of contour.
These points are valid for all thermal propagators as a consequence of the very general expression of Eq. (II.114). Since the quantities used in practical calculations depend explicitly on the choice of the contour C, it is imperative to choose a formalism where the physical quantities defined at real times can straightforwardly be extracted and in such a way that the final result does not depend on this choice of C. Two commonly used formalisms that provide reliable methods to extract such quantities will be presented in the following.

III. IMAGINARY-TIME FORMALISM
The thermal formulation of QFT presented in the sections above formally resembles to a high degree the well-known QFT at zero temperature. The main distinguishing property is that the explicit form of the thermal n-point function (Green's function) depends on the choice of integration contour C in the complex temporal plane. This contour connects some arbitrary initial time t in with the final time t in −iβ. The latter time point emerges as a consequence of the formal interpretation of the thermal equilibrium distribution operator as an evolution operator. In the equilibrium case, the contour that connects the initial and final time points can be arbitrarily chosen with the single restriction that the imaginary component of the contour must not increase when integrating along the contour. In this section, a specific choice of contour will be considered in order to make explicit calculations possible. This discussion will be followed by a diagrammatic expansion of the generating functional in order to derive the Feynman rules on this contour.

A. Matsubara contour
The simplest possible choice of the integration contour has been proposed by Matsubara in Ref. [6]. It constitutes a straight vertical line that connects t in and t in − iβ, cf. Fig. 1 with λ = β. The time variable may then be parametrised on a simple and purely imaginary form as t = iτ with τ ∈ [0, β] being a real variable. This particular choice has the advantage that it generates a perturbative expansion that contains the diagrams recognised from the corresponding zero-temperature theory (vacuum bubbles, connected diagrams, 1PI, etc.) [26]. This section will mainly be concerned with a single scalar field with µ = 0 as an illuminating example, but the results may be generalised to a theory with charged fields in a non-trivial representation of both Lorentz and internal symmetries, see Ref. [13,26].
Through a suitable parametrisation of time, the action may be written as an integral in Euclidean space: where the Euclidean Lagrangian reads The n-point Green's functions may be obtained by functional differentiation w.r.t. a source j(τ, x) of the generating functional Since the complex shift of the temporal variable is always present in theories with an equilibrated medium, TFT is in a sense intrinsically Euclidean. Let us discuss the basic ingredients of such a Euclidean formulation of Matsubara.
B. Imaginary-time propagator As stated in Eq. (II.95), the free generating functional Z 0 [j] (with V(x) = 0) may be written in terms of the free propagator by completion of the square. In the Matsubara formalism, the result takes the following form The inverse of the differential operator must satisfy the boundary condition of Eq. (II.96) in order to fix a proper solution of the differential equation. The solution is the thermal propagator of which is given by Eq. (II.115) over the chosen contour 12 . The propagator is periodic in τ over a period β as imposed by the KMS condition (Eq. (II.112)), which can be straightforwardly verified.
Completion of the square of the free Lagrangian renders the generating functional in the form Comparing with Eq. (II.95), the Euclidean propagator of a neutral scalar may be found as Note the purely imaginary temporal coordinate of the Euclidean propagator here. Thus, Z 0 [j] shown above generates all the Green's functions with imaginary time arguments. The above discussion is analogous to the zero-temperature theory, with the following two differences: • the temporal integration is reduced to the periodic interval [0, β]: i dt → β 0 dτ ; • the appropriate propagator of a free scalar field is no longer the Feynman propagator but the thermal propagator of Matsubara: Due to the periodicity condition in τ applied to the bosonic fields (antiperiodicity for fermionic fields), and therefore on the thermal propagator, over β, the Matsubara propagator may be written in Fourier space as a sum over discrete frequencies ω n . Transformation of Eq. (II.98) using the Matsubara contour leads to with ω n = 2nπ/β (bosonic), ω n = (2n + 1)π/β (fermionic), and ω 2 k = |k| 2 + m 2 . Note that the mode energy here reads ω k ≡ ω(|k|), meaning that k denotes an argument of the function rather than an index to be summed over. The ω n 's are the Matsubara frequencies [6,13,26,34] and they parametrise an infinite number of poles of the momentum-space propagator located on the imaginary axis provided that the Fourier transform is A more general multi-mass version of Matsubara propagator including charged fields was presented by Ref. [13]. In order to obtain an expression analogous to the zero-temperature theory, define k 0 = iω n in Eq. (III.7). Then where the energy k 0 is purely imaginary.

C. Imaginary-time Feynman rules
In the path-integral formulation, the generating functional may be expanded to yield the Feynman rules of the imaginary-time formalism through a perturbative series. Following the reasoning of Sec. II E 1, the Euclidean action may be split into a free part that is quadratic in the field and an interaction part so that S E = S 0 + S I . The generating functional of Eq. (II.89) may then be expanded as Each term in this series may be reinterpreted as a statistical average with respect to the weight exp[−S 0 ] so that The Feynman rules are directly extracted from this series. For the specific case of φ l (l = 3, 4, . . . ) interacting theory with a coupling constant κ, the rules are given as follows [13,35]: • Draw diagrams using Wick contractions and determine symmetry factors.
• Define k 0 = iω n and assign a propagator∆(k) to each line and a factor of −κ to each vertex i.
• Sum and integrate over all internal energies and momenta k i : 1 β n Again, one should note that the diagrammatic structure is identical to that of the zero-temperature theory. Further, the definition of the cumulant in Eq. (II.79) allows for the extraction of all connected diagrams with at least one external leg so that

D. Matsubara frequency summation
This subsection contains a short discussion on practical calculations in Matsubara theory. Every term in the perturbation theory series of log Z[j] may be expanded as a product of Matsubara propagators due to the thermodynamic Wick's theorem. Such propagator products must generally be summed up over Matsubara frequencies. An example of Matusbara summation from φ 3 -theory is the propagator product of the scalar loop diagram with one propagator for each internal line in the loop, with ω m−n ≡ ω m − ω n . The summation corresponds to the integration over time in the zerotemperature theory and may be arbitrarily complicated. A general scheme should therefore be provided in order to perform the summation over any number of propagators in the product. Generically, summation over Matsubara frequencies ω n comes in the following form with the function g(iω n ) being a product of arbitrarily many Matsubara propagators, and η indicates bosonic or fermionic frequencies i.e.
FIG. 3: Poles of the product g(z)h η (z). The crosses on the imaginary axis (red) are the infinite number of poles of the weighting function h η (z). The cross at complex z (green) symbolises the poles of g(z). The contour in (a) may be deformed as shown in (b) [35,36]. The poles of the weight function indicated in red may be used to evaluate the sum of the residues indicated in green.
In the case of a single propagator∆(iω, k), the function g reads Note, S η converges if g(z = iω n ) approaches 0 when |z| → ∞ faster that z −1 .
In order to evaluate S η , one may deploy calculus of residues in a procedure described by, for instance, Refs. [35,36]. The Matsubara sum Eq. (III.14) may be treated as the sum of residues of simple poles at iω n and, thus, may be turned into a closed integral. A suitable weighting function h η (z), different for the bosonic and fermionic cases, that has simple poles at every instance of iω n should be introduced to accomplish this. Then, The poles of h η (z) are illustrated in Fig. 3. A number of crosses on the imaginary axis in the figure indicate the poles of the appropriately constructed h η (z) at iω n . The contour picks up the residues of all these poles giving rise to S η . To perform the Matsubara summation, the contour lines are deformed to enclose the poles of g(z), see Fig. 4. The sum S η is then the sum of the residues of g(z)h η (z) over all poles of g(z): The leading negative sign is due to the clockwise orientation of the contour. This procedure may be carried through given that g(z) has no poles on the imaginary axis [35].

Choice of hη(z)
The choice of a suitable Matsubara weighting function h η (z) is a somewhat delicate matter. In order to guarantee simple poles at z = iω n , for example, for bosonic frequencies ω n = 2πn/β, any of the following weighting functions may be utilised: + (z) may be used to control convergence in the half-plane with z < 0, while h (2) + (z) controls convergence in the right half-plane where z > 0.
Adding the two functions to provide convergence in the entire complex plane yields and this relation will (by means of the contour deformation) further reduce down to In the case of the tadpole diagram in φ 4 -theory This function clearly has simple poles at z = −sω k and the two residues (one for each term in the sum) are Res g(z 0 ) = s/2ω k due to the structure of the Laurent series of g(z). Hence, Employing the following identities, the above equation simplifies to (III.26)

Example: eye diagram
In the case of the eye-diagram of φ 3 -theory discussed further down in the context of thermal decay rates, the sum S η appears to be somewhat more complicated. In this case, The poles of g(z) are simply the four poles z 0 = −sω q and z 0 = iω n + s ω k−q . Hence, accounting for the residues over these poles, the Matsubara frequency sum takes the form Noticing that and taking out signs of the coth-functions in order to cancel the sign in front of the sum, we arrive at . (III.30) Then, expanding the sum over s, s and using Eq. (III.25), one gets the final result for the Matsubara summation: that can be exploited in practical calculations.

IV. THE REAL-TIME FORMALISM
Matsubara found that, by the procedure of analytic continuation of the Euclidean Green's functions, the retarded correlation functions with real time-arguments may be obtained in the imaginarytime formalism presented above. However, the evaluation of three-momentum integrals over Matsubara propagators summed over the frequencies is by no means straightforward for any but the simplest cases. By evaluating the Green's functions with respect to a contour that contains the real axis, in the so-called real-time formalism, an analytic continuation is avoided completely. Such a real-time procedure is originally due to Schwinger [37] and Keldysh [9]. The formalism makes use of a special contour of temporal integration, the real-time contour presented in Fig. 2, which has the effect of doubling the degrees of freedom. This property, in particular, results in a 2×2-matrix structure of the propagator for a given field to be discussed below.
Other contours have been investigated as well. Particularly, one could mention the suggestion by Evans [38] to use an asymmetric contour.

A. Keldysh contour
The real-time contour in Fig. 2 is a special case of a family of real-time Keldysh-esque contours [10,24] parametrised by σ ∈ [0, β] and presented in Fig. 5. The Keldysh contour (or sometimes Kadanoff-Baym-Keldysh contour, see e.g. Ref. [39]) is recovered for σ = 0. An even more general family of contours was presented in Ref. [24] with pieces going back and forth N times between t in and t fi in parallel to the real axis. However, it can be shown that such additional pieces cancel out for N ≥ 2 and the fundamental reasons for discarding this extended class of contours are laid out in Ref. [24]. Hence, the contour presented in Fig. 5 is the most general real-time contour needed to include all n-point functions with real arguments, with σ being a freedom of choice. Commonly, this parameter is set to 0, β or β/2 [12,13,21,40].
It was argued previously, using Eqs. (II.68) and (II.114), that the contour propagator depends explicitly on the choice of contour C. It is clear from these equations that the contour dependence of the propagator enters only through the contour step-function Θ C (t − t ). In the case of the real-time contour shown in Fig. 5, the time coordinates of x, x may be distributed arbitrarily over C and several unique two-point functions (propagators) arise. Let us label them as t r , t s so that t r ∈ C r , t s ∈ C s for r, s = {1, 2, 3, 4}. Then, the different propagators may be expressed as (IV.1) Note e.g. that t 2 , t 3 always are later times than t 1 , t 4 .
The contour is extended along the real axis in order to pick out all the real time-arguments of the n-point function of interest. It then goes down by iσ before going back in parallel to the real axis. Finally, the contour goes down to the final time t in −iβ. The points t in and t fi are arbitrary and may be suitably chosen for any given thermal system. In the equilibrium limit, these end-points are taken to ±∞.
The appearance of several distinct propagators is a very general consequence of the statistical formulation of the two-point function in Eq. (II.16). The real-time contour of Fig. 2 is treated in great detail by Wagner [18]. In an approach valid both in as well as out of equilibrium, Wagner organised the two-point functions (Green's functions) in a matrix. In the special case of σ = 0, this matrix takes the following form Adding the piece C 4 (σ = 0) is a trivial extension. For any system in equilibrium, Landsman and van Weert [13], Wagner [18] and Das [26] argued that the contour pieces C 3 , C 4 decouple entirely from the parallel pieces due to in the limit t in → −∞ and t fi → +∞. This may be proven by the application of the Riemann-Lebesgue lemma [41] if the source field is constrained to This condition is compatible with the KMS condition [13]. Due to the decoupling of the vertical contour segments, the generating functional factorises as Here, C rs = C r ∪ C s . Therefore, the main focus, when considering the real-time Green's functions, is the upper-left 2×2-corner of the matrix in Eq. (IV.2) while the contribution from the vertical pieces is treated as a multiplicative factor. This factor only enters the connected Green's functions as an additive term in the cumulant expansion of Eq. (II.77).

B. Real-time propagator
An explicit form of the real-time thermal propagator, in the case of a neutral scalar particle, can be found by examining Eq. (IV.6) in the light of Eq. (IV.5). Having factorised out the vacuum bubbles Z[0] and the vertical contour pieces Z 34 [j], only the real-time part Z 12 [j] will be considered from here on. The two exponentials in Z 12 [j] can be further factorised with respect to C 1 ∪ C 2 by conveniently organising the propagator components in a matrix. From now on, the spatial variables will be suppressed. At the same time, a new notation will be introduced. The contour pieces C 1 , C 2 in Fig. 5 are relabelled as C + , C − , respectively, to indicate chronological and antichronological time evolution on the different parts of the contour.
The factorisation of the first exponential becomes . (IV.7) The relative minus sign is a result of the antichronological evolution of the fields on C − . The two functions here j + (t), j − (t) are introduced as sources of two separate field components φ + , φ − that exist only on C + , C − , respectively. Together, the components form the full real-time field The physical field component φ + evolves chronologically and its integration limits are recognised from the zero-temperature field theory. The component φ − is an unphysical degree of freedom since it evolves antichronologically. This field component can not appear on external lines of diagrams generated by the above exponential. Nevertheless, its presence is an unavoidable consequence of the Keldysh contour, and contributions from virtual φ − -fields in loops must be taken into account. Such contributions are introduced by the second term L I [δ/iδj − ]. Summing over r, s = ±, the second exponential factor can be rewritten as The matrix components of the propagator must fulfil (IV.10) Here, the internal indices i, j have been suppressed, and is the Feynman propagator of the zero-temperature theory. Recall that the spectral density ρ 0 is the discontinuity over the real axis of this propagator (Eq. (II.99)) and may therefore also be written in terms of i∆ F , i∆ * F as Hence, the propagator matrix may be diagonalised in a basis of zero-temperature functions by means of a Bogolyubov transformation. The result is with the thermal (Bogolyubov) matrix defined as follows (IV.14) Here the thermal angle θ k was introduced 13 through the functions Here, the thermal function reads N (ω) = Θ(k 0 )n(ω + ) + Θ(−k 0 )n(−ω + ). The square-parenthesis notation conveys the point that for fermions the goniometric functions are considered while the hyperbolic functions apply to the scalar case. Note that the propagator is symmetric under simultaneous transposition and reversal of momentum, i.e. iD rs (k) = iD sr (−k), for neutral bosons 14 .
In the limit of T → 0, the Bogolyubov matrix becomes the unit matrix. In such a case, the field components completely decouple such that Z +− → Z + Z − and the contribution from φ − may be fully absorbed into the normalisation constant. The diagonalisation of the propagator is particularly useful since any thermal dependence may be absorbed into vertices giving rise to a basis where the propagator itself is defined in terms of zero-temperature functions. Hence, the real-time formalism and the propagator are intrinsically connected to the imaginary-time formalism since they meet when the formalism of Matsubara is analytically continued to real time-arguments. This statement is formulated as follows where ∆ is the analytic continuation of the Matsubara propagator of Eq. (III.7) down to the real axis.

Real-time scalar propagator(s)
The thermal decays of interest here involve interactions between fermions and scalar bosons. In this subsection and the following one, the real-time free propagators for such fields will be explicitly presented for further references.
For the choice of σ = 0, the components of Eq. (IV.10) may be written as Here, one recognises the retarded and advanced propagators known from the zero-temperature theory. For practical calculations, it is common to rewrite the four components in terms of the Feynman propagator (and its conjugate), consequently rewriting also the thermal terms. Through the relation 19) and by insertion of the explicit spectral density Eq. (II.99), the components can be found as Here, the new thermal distributions introduced for compact notation read Note especially that, in the case of a neutral boson, one might define and, for further convenience (see Secs. VII-X), in the case of charged fermions we have Here, the lower sign reflects the opposite charge of antiparticles.

C. Real-time Feynman rules
It was seen in the previous section that the Keldysh contour in Fig. 2 necessarily gives rise to a doubling of the degrees of freedom in the real-time thermal theory. This is achieved by separating the field into two components, each appearing on the respective segment of the real-time contour. One component is physical while the other should be regarded as a thermal ghost arising from the choice of the specific contour. For notational convenience, a neutral scalar field φ will be considered in the derivation of the Feynman rules below. For illustration, the field is assumed to have an interaction term proportional to κφ l . The generality of the method, however, allows for the treatment of the multicomponent vector field I ∞ i α introduced above. The derivation of the Feynman rules for the case of I ∞ i α is analogous to the case of a neutral scalar theory, but with the presence of contraction of Lorentz and internal indices in addition to the summation over the thermal propagator components.
It was seen in Sec. IV B that, for a field with no derivative couplings, the relevant factor in the generating functional, Z 12 [j], can be written as The series expansion of this function generates all the diagrams of the theory in analogy to the zerotemperature case. Two distinct vertices appear due to the two interaction terms in the leftmost exponential: one interaction vertex involves only the physical field component and one vertex involves only the ghost field. See Fig. 6 for the vertices of the φ 4 -theory as an instructive example. Hence, the field components do not mix in the vertices. However, due to non-zero off-diagonal components of the propagator on the Keldysh contour, the fields may propagate into each other. The connected Green's functions of Eq. (II.77) are generated by functional differentiation according to (IV.27) It is clear that, for any value of σ = 0, the real-time Green's functions may only be recovered by functional differentiation with respect to the source j + since the second contour segment includes no real times. It may be proven [13] that the real-time Green's functions are independent of the parameter σ and, hence, they are given by ∀σ . (IV.28) The thermal ghost field φ − may only propagate in internal lines that connect the vertices. These lines are generated by functional differentiation with respect to j − inside Z[j + , j − ], see Eq. (IV.7). The full set of Feynman rules can easily be extracted through the generating functional in its path-integral form if the Lagrangian is split into its free and interaction parts. Then, in terms of the action, the generating functional may be written as Summation over r, s = ± is assumed. Expanding the second exponential, each term may be interpreted as an expectation value in the path-integral formalism with respect to a weight e iS0 so that By means of the cumulant expansion, it is possible to factor out not only the free Z 0 [0, 0] but all vacuum bubbles Z[0, 0] so that the cumulant becomes The expansion of the final term gives rise to all connected graphs with at least one external leg of the physical field φ + . The real-time Feynman rules are extracted from this series as follows: • Draw diagrams using φ + -, φ − -vertices: n external φ + -lines and arbitrarily many internal φ +and φ − -lines. Determine symmetry factors.
• Integrate over each internal momentum k i : 4 and sum over all internal distributions of r, s.
Comparing the real-time formalism to the Matsubara theory, the propagator has gained a matrix structure and, further, two separate types of vertices have appeared.

V. THERMAL SELF-ENERGIES
This section discusses the expression of the thermal decay rate in terms of the self-energy. The topic has been treated by e.g. Weldon [14], Ho and Scherrer [34] and discussions on the thermal self-energy has also been advanced by [13,15,16,25,26,40,42,43]. Weldon [14] related the self-energy of the Matsubara formalism to a quantity interpreted as the thermal decay rate. In the sections below, we will examine the self-energy of Matsubara, its analytic continuation and finally the connection of real-time quantities to it. This connection has been discussed by several authors, see for example Refs. [13,25,40].

A. Matsubara self-energies
The connected two-point function of the Matsubara formalism is obtained by functional differentiation of the cumulant in Eq. (III.12). In order to obtain a perturbative series, and hence the Feynman rules, in momentum space, the full generating functional of Eq. (III.11) may be Fourier transformed. It may be further manipulated in order to recover a form similar to that of Eq. (II.91). The full propagator appearing in Feynman diagrams is introduced as staring from its zero-temperature counterpart. Since global momentum conservation is imposed for all diagrams, an overall factor of β(2π) 3 δ n0 δ k i may be extracted, leaving the interacting manybody propagator labelledG(k) by Landsman and van Weert [13]. Note here that the propagator carries the discrete energies in terms of the Matsubara frequencies. Real-time quantities must be obtained through analytic continuation. Analogous to Eq. (III.8), this propagator may be defined in terms of the full Matsubara propagator that has been extended away from the discrete frequencies by means of analytic continuation:∆(iω n , k) →∆ (z, k). Hence, For the full propagator, the free spectral density ρ 0 has been replaced by the corresponding full quantity defined as iρ(k) = Disc∆ (k). The procedure of analytic continuation allows the propagator to be defined for real energies. This continuation is not unique and that is usually resolved by letting lim |z|→∞∆ (z, k) = 0 and by taking∆ (z, k) to be analytic away from the real axis [13,27] which provides the above expression. Thus, we have uniquelỹ provided the property k 0 ρ(k) ≥ 0. This function can be shown to have neither zeroes nor poles off the real axis provided that this inequality property holds true [13]. Given that the extended propagator is analytic off the real axis, its analytic inverse exists and fulfils a Dyson-like equatioñ which defines the analytic self-energy. Note the free analytic propagator on the right-hand side. Guided by the above consideration, it is useful to define the thermal Feynman propagator as This propagator may be inverted by writing k 2 0 − |k| 2 = m 2 so that the thermal self-energy for the imaginary-time formalism can be extracted: As in the zero-temperature theory, one may observe that the real part of the self-energy may be absorbed as a correction to the mass while its imaginary component will be interpreted as a decay rate [14].

B. Real-time self-energies
The full real-time propagator may be assumed to satisfy the Schwinger-Dyson equation (see Ref. [13] and references therein) similar to the zero-temperature case: Like in the zero-temperature theory, the real-time self-energy is most often 15 the sum of all oneparticle irreducible (1PI) insertions on the propagator line as visualised in Fig. 7. This self-energy coincides with the self-energy of the many-body propagator provided by Ref. [13] in connection to their discussion on the Matsubara theory. The full real-time components are not independent as a consequence of the KMS condition and due to translational invariance [13] and they satisfy similar relations as the components of Eq. (IV.10).

FIG. 7:
The graphical interpretation of the Schwinger-Dyson equation. All 1PI insertions are parametrised by the self-energy. By recursive insertion of the full propagator, the 1PI series expansion may be obtained.

C. Self-energy relations
Viewing the free-field propagator of Eq. (IV.10) as a special case of the full propagator of Eq. (V.7), the matrix of Eq. (IV.14) also diagonalises the full propagator for some functionD Fαβ (k) so thatD The new function is assumed to satisfy the Schwinger-Dyson equation that defines the self-energȳ Π F :D Fαβ (k) =D Fαβ (k) + D FΠFDF αβ . As an example, the relation betweenD Fαβ (k) and the componentD (++) 16 A 1PI diagram is defined to be any diagram that cannot be split into two parts by removal of a single propagator line. The 1PI blob shown in Fig. 7 is defined as the sum of all such diagrams with the eye-diagram (the single loop) being the lowest-order contribution. See Ref. [1] for a discussion of the 1PI insertion procedure. One must note, however, that the properly defined self-energy includes one-particle reducible diagrams as well, e.g. tadpole diagrams that can be regulated separately.

The insertion of M η and its inverse in the Schwinger-Dyson equation allows for the identificatioñ
The derivation of this is rather lengthy but can be found in Ref. [13] where, most importantly, it is argued thatΠ F is indeed the analytically continued self-energy of the re-summed propagator∆ F in Eq. (V.6). For notational convenience, possible spinor and Lorentz indices have been suppressed but are generally present in the structure of the Klein-Gordon divisor d ij αβ . By inserting the explicit expression for the thermal matrix M η , the following relations may be deduced:Π It is similarly straightforward to relateΠ (−+) toΠ (++) . Hence, a number of components has reduced from the initial four to a single independent one. The two relations above hold for any σ, also for σ = 0. Importantly for the forthcoming calculations, the real and imaginary parts of the self-energies are related as 14) The self-energy on the left-hand side will be calculated to the first order for several types of interaction vertices in Secs. VII-X.

VI. THERMAL DECAY RATES
In this section, the self-energies of the previous section will be related to the thermal decay rates. Initially, a relation between the thermal quantity that corresponds to the zero-temperature decay rate and the self-energy of the imaginary-time formalism due to Weldon [14] will be presented. This relation will then be restated in terms of the self-energy components of the real-time formalism by making use of the self-energy relations in Sec. V C.

A. Thermal rates in Matsubara formalism
The decay rate γ D for a given process in the zero-temperature theory may be related to the discontinuity of the self-energy (i.e. its imaginary part) by means of the optical theorem as where E 0 is the energy of the decaying particle. Weldon [14] defined a similar quantity Γ(p 0 ) through whereΠ F is the self-energy of the imaginary-time formalism. It is important to reflect upon the physics that this quantity describes. One may assume that the distribution of a particle Φ at some time t in is described by a nonequilibrium function which Weldon labelled as f (p 0 , t in ) [14]. This function will approach the thermal distribution of the equilibrium (Bose-Einstein or Fermi-Dirac) in a simple manner if a deviation from equilibrium is small (∂f /∂t 1). The rate of this approach is parametrised by Γ(p 0 ) to the first order. The evolution of this distribution may be formulated according to in terms of the quantities Γ D , Γ I . Here, the first term expresses the loss of Φ-particles through decay modes while Γ I takes into account particle production by the medium. The thermal medium is filled with particles that couple to Φ and, hence, the second parameter Γ I , the inverse decay rate (production rate), adds the contribution from processes in the medium that produce Φ particles. As an example, production may occur through reactions φφ → Φ. Weldon comprehensively presents an analysis of possible production and decay channels for Φ through interactions with particles in the thermal medium. The solution to Eq. (VI. 3) is for some function c(p 0 ) which is constant in time. This solution requires that the temperature T is constant over time which is true for equilibrated media. Since T characterises the background medium, the deviations of the Φ distribution from its equilibrium limit are required to be small and one may assume that the distribution of the medium particles corresponds to the thermal equilibrium distribution. The net decay rate Γ of the distribution of Φ in the medium is the rate at which the distribution of Φ approaches the equilibrium regime and amounts to The evaluation of the amplitude for the thermal forward decay Γ D is of interest in this work. Weldon provides a thermal relation between the forward and inverse decay rates due to unitarity so that one therefore may write the decay rate of Φ as The expression provided by Weldon has hereby been extended in order to explicitly take into account the role of the chemical potential µ.

B. Thermal rates in the real-time formalism
The relations of Sec. V C connect the self-energy of the Matsubara formalism to the components of the real-time self-energy. Hence, Eq. (VI.6) may be rewritten as In the following sections VII-X, the (++)-component of the real-time self-energy for several types of interactions has been evaluated. The results are used to extract the decay rates for fields in an equilibrated thermal medium. Note that, ifΠ (++) comes with internal or Lorentz indices, one may follow the procedure advised by Weldon [14] and define the scalar function which represents a contraction with asymptotic states. This allows for a probabilistic interpretation of Σ as a decay rate, similar to a Breit-Wigner resonance, for any particle species. For future convenience in the following sections, the ratio of the thermal decay rate above to the zero-temperature limit, as used in the analysis by Ref. [34], may be defined as This quantity is a parametrisation for the deviation of thermal theory predictions from those in the zero-temperature limit. If this ratio is equal to unity for all temperatures, the medium has no effect on decay rates.

VII. (PSEUDO)SCALAR DECAY INTO TWO (PSEUDO)SCALARS
We are now equipped with the formalism and the relations needed for computation of thermal decay rates. This section provides a detailed analysis of the thermal decay rate at which a neutral (pseudo)scalar particle Φ of mass M decays into a pair of neutral (pseudo)scalar particles φ 1 φ 2 of masses m 1 , m 2 . The scalars in the final state may not be identical in the following expressions. Note, however, that all scalars considered here are neutral particles. The interaction term responsible for such a decay process reads (VII.1) The resulting scalar boson self-energy discussed in this section has been previously published in Ref. [40] and represents an important example of the usage of the real-time TFT framework.

A. Real-time self-energy of the scalar-scalar eye-diagram
Using the interaction term of Eq. (VII.1), one may draw the following eye-diagram: .

(VII.2)
This is the matrix eye-diagram that appears at the first order in the Dyson series representing the full interacting propagator of Φ. The crosses denote either external legs or connections to vertices.
The Φ-lines in the loop diagram are amputated in the following self-energy calculation. Note that this diagram represents all allowed propagator combinations 17 . The intermediate particles are not restricted to the physical field component φ + while the external legs are. As a consequence, all valid combinations of propagators in the loop must be considered when calculating the total amplitude for this diagram. The (+)-and (−)-field components do not mix in the vertices, cf. Fig. 6, but propagate into each other due to the off-diagonal propagator elements. Besides, each vertex contributes with a factor of ∓iκ, respectively. In order to obtain the thermal decay rate, the self-energy of Eq. (VII.2) will be evaluated below and it is defined by the loop diagram   Note here the pre-factor; δ 12 is defined as unity if the two particles in the loop are identical and zero otherwise. Hence, the correct symmetry factor will be taken into account when considering either identical or distinct loop particle species. The two-component structure of the scalar propagator splits the bubble into four terms indicated above: a non-thermal (temperature independent) integral I The evaluation of the types of integrals above is discussed in some detail in Sec. B.

Non-thermal self-energy term
The non-thermal contribution to the self-energy is well-known from zero-temperature QFT. It may be evaluated by dimensional regularisation techniques in D = 4 − 2 dimensions at the renormalisation scaleμ. The resulting expression for the non-thermal integral is Here, 1/˜ parametrises the divergence emerging from the dimensional regularisation scheme and, denoting the Euler-Mascheroni constant by γ E , it is defined as 1 = 1 − γ E + ln 4π. The last term I loop is the term of interest to the decay rate since it is complex. (VII.9) The cases are defined as the momentum regions (VII.10) The imaginary part of Eq. (VII.8) reads as follows (VII.11)

Mixed self-energy term
The mixed self-energy contribution consists of two complex propagator cross terms. The evaluation of the imaginary component of those integrals results in (VII.12) Here, Two further cases have been defined as The real part of the mixed term is (VII.14) Here ω 2 2,k = |k| 2 + m 2 2 .

Thermal self-energy term
The thermal contribution to the self-energy has no real component and is given by the purely imaginary expression The further two cases were defined as The leftmost inequality corresponds to the upper sign. The rightmost inequality corresponds to the lower sign.
The leftmost curly bracket corresponds to the upper sign. The rightmost curly bracket corresponds to the lower sign.
It is important to notice that in the above expression for a purely thermal contribution, the mass parameters m 1 , m 2 in ω 2,p± should be replaced by expressions min{m 1 , m 2 } and max{m 1 , m 2 }, respectively. The combined (++)-self-energy component of the eye-diagram given in Sec. VII A 1-VII A 3 was provided by Nishikawa et al. [40] as part of the evaluation of the sunset diagram and has been included here for completeness.

B. Decay rate of Φ → φφ
This section presents the thermal decay rate of a scalar particle into two scalar particles with the rate given by Eq. (VI.7), in Sec. VI B. In the case of identical loop masses, the expression for the self-energy simplifies significantly due to vanishing terms in the mass function √ C and ω 2,p± . The ratio of Eq. (VI.9) may be obtained by considering the limiting case of p 2 = M 2 (m 1 + m 2 ) 2 , rendering the particles in the loop effectively massless. The plotted ratio can be seen in Fig. 8a-c for a non-relativistic, relativistic and ultra-relativistic incoming particle Φ. The figures reproduce the findings of Ho and Scherrer [34] who considered the case of identical loop masses m 1 = m 2 = 0 evaluated in the imaginary-time formalism and their result has hereby been verified in the real-time formalism. The limit as T → 0 for the ratio is R Φ→φφ → 1 as required. Ref. [34] that were evaluated in the Matsubara formalism. Note that R → 1 when T vanishes.
The deviation from Ref. [34] when loop masses m 1 , m 2 are not identical was investigated. Such deviations are clearly relevant only in the case of M ∼ m 1 , m 2 , close to the equality of M 2 = (m 1 + m 2 ) 2 . An analysis of the self-energy presented earlier in this section shows that the behaviour is qualitatively similar to Fig. 8a-c. Furthermore, in the mass region 0 ≤ M 2 < (m 1 − m 2 ) 2 and for virtual Φ with p 2 < 0, the behaviour of Γ Φ→φφ is also quadratically growing with temperature. Hence, these regions exhibit a qualitatively similar behaviour compared to the plotted cases.

VIII. SCALAR DECAY INTO A FERMION-ANTIFERMION PAIR
This section presents the thermal rate at which a neutral scalar particle Φ decays into a fermionantifermion pair ψ 2ψ1 . One model that gives rise to such decay process is (VIII.1) The particles are associated with masses M for Φ and m i for ψ i .

A. Real-time self-energy of the fermion-antifermion eye-diagram
In analogy to the case of the scalar-scalar loop considered in the previous section (Sec. VII), the self-energy of the eye-diagram is related to the thermal decay rate of the process Φ → ψ 2ψ1 . Also here, all valid combinations of field components must be taken into account but only one self-energy component is independent, see Eq. (V.10), and the evaluation of the self-energy (++)-component is sufficient to extract the thermal decay rate. This component reads Note, the overall negative sign appears together the trace over spinor indices in the case of a fermion loop. The bracketed arrow represents the second mixed integral that arises from propagator cross terms. As in the case of the scalar-scalar loop, the bubble has split into four terms and the explicit integrals are The trace is common to all integrals and evaluates to The non-thermal contribution is again known from the zero-temperature theory and the imaginary part of interest to the decay rate calculation may be extracted. The imaginary part of iI where C is identical to the definition in Sec. VII. In the case of identical loop masses, the expression reduces to This matches the corresponding result in the literature, see e.g. Ref. [44].

Mixed self-energy term
The mixed self-energy contribution consists of two complex cross terms. The evaluation of the imaginary component of those terms results in (VIII.9) Here, ω 1,p± has been defined analogously to ω 2,p± in Sec. VII with interchanged masses.
In compact notation, the real part may be written as where

Thermal self-energy term
The purely thermal contribution to the self-energy has no real component and may be found as (VIII. 16) Note that in this final expression for the thermal contribution, the mass parameters m 1 , m 2 should again be replaced by min{m 1 , m 2 } and max{m 1 , m 2 }, respectively; cf. Sec. VII A 3. The explicit (++)-component of the real-time self-energy of the fermionic eye-diagram presented above in Secs. VIII A 1-VIII A 3 has not been found in the literature to the best of our knowledge.

B. Decay rate of Φ → ψψ
The thermal decay rate of a neutral scalar is given by Eq. (VI.7) in Sec. VI B. Similar to the decay rate treated in Sec. VII, the particles in the loop were considered to have the same mass m and the limit of p 2 = M 2 4m 2 has been plotted in Fig. 9. In this mass region, the ratio of Eq. (VI.9) may be obtained since the zero-temperature decay rate does not vanish here. The plotted ratio can be seen for non-relativistic, relativistic and ultra-relativistic Φ. The figure agrees with the findings of Ho and Scherrer [34] who evaluated the loop in the Matsubara formalism. The limit as T → 0 is R Φ→ψψ → 1 as expected.

IX. PSEUDOSCALAR DECAY INTO A FERMION-ANTIFERMION PAIR
This section presents the thermal rate at which a neutral pseudoscalar particle Φ 5 decays into a fermion-antifermion pair ψ 2ψ1 as calculated in this work in the real-time formalism. This case was explicitly considered in order to contrast the pseudoscalar result with the result for a scalar particle provided in Sec. VIII. The corresponding interaction term is Here, the vertex factor of γ 5 requires some care when evaluating the loop diagrams in comparison to the previous section. In what follows, masses have been assigned according to M 5 for Φ 5 and m i for ψ i .

A. Real-time self-energy of the fermion-antifermion eye-diagram
In analogy to the case of the scalar-scalar loop in Sec. VII, the self-energy of the eye-diagram is related to the thermal decay rate of the process Φ 5 → ψ 1ψ2 . Also here, all valid combinations of field components must be taken into account but only one self-energy component is independent, see Eq. (V.10), and the evaluation of the self-energy (++)-component is sufficient to extract the thermal decay rate. This component is found as follows  denotes the second mixed cross term. The emerging trace may be taken in the naïve dimensional regularisation scheme since the diagram appears at one-loop level 18 . The trace evaluated in the naïve dimensional regularisation scheme reads The sign differences compared to the case of a decaying scalar particle are an overall positive sign together with the minus sign in front of m 1 m 2 . The overall sign change may be reverted by imposing hermiticity of the Lagrangian and thereby force the pseudoscalar coupling toψ 1 γ 5 ψ 2 to become purely imaginary, i.e. b = i|b|. The sign difference of the mass term is more intriguing and signifies a fundamental difference between the scalar and pseudoscalar decay rates. Apart from the mentioned differences in relative signs, the resulting self-energy integrals for the pseudoscalar case are identical to those for the scalar case. To obtain the correct non-thermal, mixed and purely thermal self-energy integrals it is sufficient to replace m 1 → −m 1 in the integrals from the previous section.

Non-thermal self-energy term
The non-thermal contribution in the case of an external pseudoscalar particle is related to the case of an external scalar particle as (IX.4) The right-hand side is the result found in Sec. VIII A 1.

Mixed self-energy term
The mixed thermal contribution in the case of an external pseudoscalar particle is related to the case of an external scalar particle as (IX.5) The right-hand side is the result found in Sec. VIII A 2.

Thermal self-energy term
The purely thermal contribution in the case of an external pseudoscalar particle is related to the case of an external scalar particle as  (IX.6) The right-hand side is the result found in Sec. VIII A 3.

B. Decay rate of Φ5 → ψψ
The thermal decay rate of a neutral pseudoscalar is given by Eq. (VI.7) in Sec. VI B. Similar to the decay rate treated in Secs. VII and VIII, the particles in the loop were considered to have the same mass m and the pseudoscalar was assumed to fulfill the on-shell condition p 2 = M 2 5 4m 2 . The decay rate behaves qualitatively similar to that of a scalar decay (see Fig. 9) as stated by Ho and Scherrer [34] while one should note a discrepancy between the two cases close to the equality M 2 5 , M 2 = 4m 2 . This is a non-thermal effect due to the mentioned sign differences in the above subsection and it is presented in Fig. 10 for clarity, at T = 0, with normalised masses M , M 5 given in units of m. The ratio Γ Φ5→ψψ /Γ Φ→ψψ is plotted as a function of the mass of the external (pseudo)scalar. Close to the equality M 2 = 4m 2 , the expression diverges which implies that the scalar decay rate vanishes quicker than the pseudoscalar decay rate. Normalising the fermion mass to m = 1, the ratio diverges at threshold M = 2m.

X. EMISSION OF A (PSEUDO)SCALAR OFF A FERMION PAIR
The interaction term of Eq. (VIII.1) in Sec. VIII allows for the emission of a scalar off a fermion line. The amplitude for such emission may be evaluated by considering (pseudo)scalar corrections to the fermion line. In this section, the thermal self-energy of the scalar-fermion one-loop diagram is presented.
The four integrals extracted above correspond to the non-thermal contribution, the mixed thermal contribution and the purely thermal contribution to the self-energy. The second mixed term, expressed through the bracketed arrow, is similar to its preceding term but with the thermal term of the fermion propagator and the non-thermal term of the boson propagator interchanged. The explicit forms are Relevant to the emission is the imaginary part of the non-thermal contribution to the (++)component of the self-energy. It may be extracted by a similar procedure performed in Sec. VIII A 1. Evaluated by the dimensional regularisation techniques, the imaginary part of Eq. (X.4) is found to be Here, (X.9)

Mixed self-energy term
The mixed self-energy contribution consists of two complex cross terms from the product of the propagators: the integrals in Eqs. (X.5) and (X.6). The resulting imaginary component of the sum of the integrals is a rather long expression when written down explicitly. For clarity, the result is therefore presented for each momentum region separately.
a. If where Li 2 (z) is the analytic continuation of Spence's function, the dilogarithm which cancels exactly the imaginary part that arises due to the negative sign of 1−e βω M,p ± in the two analytically continued logarithmic terms in Eq. (X.10). Here Here, the upper and lower signs correspond to the modified case (d2) defined in Sec. VII A 3. The modification is analogous to the modification of the case (d1) in the section above for 0 ≤ p 2 < (m 1 − M ) 2 and, as a consequence, the result depends on the mass hierarchy of M , m 1 .

The purely thermal self-energy term
The thermal self-energy contribution in Eq. (X.7) is purely imaginary, cf. previous self-energy evaluations in Secs. VII A 3, VIII A 3 and IX A 3. The resulting expression for this term, after momentum integration, is again rather long. Therefore, different regions of p 2 are presented separately below for clarity.
a. If Again, the upper and lower signs correspond to the modified case (d1) discussed in Sec. X A 2 in order to take into account the hierarchy between M , m 1 .
d. If p 2 < 0: Also here, the sign cases correspond to a modified (d2). Note that in the final expression for the purely thermal contribution, Eqs. (X.16)-(X.19), the mass parameters m 1 , M should be replaced by min{m 1 , M } and max{m 1 , M }, respectively; cf. Sec. VII A 3.
The explicit (++)-component of the real-time self-energy of the fermionic eye-diagram presented in Secs. X A 1-X A 3 above is not present in the literature to the best of our knowledge.

B. Pseudoscalar emission off the fermion line
Analogously to Sec. IX, the pseudoscalar interaction termψγ 5 ψ provides an additional vertex factor of γ 5 compared to the case of scalar emission considered in the previous three sections. In the naïve dimensional regularisation scheme, γ 5 anticommutes with the γ-matrix in the numerator of the fermion propagator and thereby provides an overall sign change together with a sign change of the numerator mass m 1 . These sign changes are stated relative to the obtained results of Secs. X A 1-X A 3 and a conclusion completely analogous to Sec. IX is reached.
C. Scalar emission rate of the process ψ 2 → ψ 1 Φ The thermal rate for the scalar emission provided in this section is given by Eq. (VI.7) in Sec. VI B. In order to extract an observable from the self-energy presented in Secs. X A 1-X A 3, that carries spin structure, it is advisable to follow the strategy provided by Weldon [14] to define This is the self-energy contracted with incoming and outgoing states. These states satisfy the Dirac equation (p /−m)u(p) = 0 as well asū(p)u(p) = 2 p 2 . In terms of Σ, the discontinuity over the real axis of this new scalar function may be related to the thermal decay rate. More specifically, the contraction of the three different pre-factors that appear in the self-energy and which contain explicit γ-matrices, is of interest. Assuming the external ψ 2 to be on-shell (p 2 = m 2 2 ), the contraction results inū The ratio of Eq. (VI.9) may be obtained for example in the mass region of p 2 = m 2 2 (m 1 + M ) 2 effectively rendering the particles in the loop massless. The plotted ratio can be seen in Fig. 11 for a non-relativistic, relativistic and ultra-relativistic incoming particle ψ 2 . In the limit T → 0 the decay ratio R ψ 2 →ψ 1 Φ → 1.
The linear dependence on T of the self-energy seen in Fig. 11 is initially surprising in comparison to the quadratic temperature dependence of the self-energy of the scalar-scalar loop seen Fig. 8. This difference in behaviour may be explained in the high-temperature limit. The dependence on temperature of the free propagator comes from the term proportional to the thermal distribution function n(|k 0 |), see Eq. (IV.20). Expanding this function in the high-temperature limit (or small β) results in different powers of T of the leading term for bosons and fermions, respectively. For bosons, the temperature-dependence of the propagator is proportional to T while the leading term for the fermion propagator is constant in T . Hence, it is seen that for T → ∞, the scalar-scalar loop of Sec. VII contains a term that is quadratic in T (the purely thermal part) while the scalar-fermion loop of Sec. X contains terms proportional at most only to T . (m 1 + M ) 2 . Curves display ratios for varying = |p|/m 2 parameter for a non-relativistic ∈ {0.001, 1}, relativistic = 10 and ultra-relativistic = 100 incoming ψ 2 . Note that the ratio R → 1 when T → 0.

XI. SUMMARY AND OUTLOOK ON THERMAL DECAYS
Thermal effects on key observables in particle physics have caught a strong interest and attention of both theorists and experimentalists over past few decades. The key results of the thermal field theory (TFT) appear to be scattered over many different sources in the literature, and a complete and consistent compendium of these results is yet missing. This review provides a pedagogical overview of the most important elements of the equilibrium TFT in both imaginary-and real-time formulations, together with the most immediate and relevant applications for all possible particle decays in thermal medium in a generalised Yukawa theory. A special attention is given to the real-time formulation and its predictions for particle two-body decay rates that are compared and verified against the existing results in the imaginary-time (Matsubara) framework, widely used in the literature. This represents an important step towards a proper generalisation of TFT towards in-equilibrium particle dynamics that can be consistently achieved in the real-time framework.
Starting from a theoretical overview of the very general and basic principles of TFT formulated by Wagner [18], we present a collection of possible decay rates for a theory containing generic (pseudo)scalar bosons and fermions. Several thermal one-loop calculations were performed within the scope of this work in the real-time framework. Components of the self-energy were related to a quantity interpreted as a thermal decay rate of a field in an equilibrated medium which is considered in various possible kinematical domains.
A (pseudo)scalar state decaying into a pair of (pseudo)scalar particles was considered as well as a (pseudo)scalar decaying into a fermion-antifermion pair. The thermal rate of a (pseudo)scalar emission off a fermion line is also provided. For all these processes, except for the (pseudo)scalar decay into a fermion-antifermion pair, an enhancement of the thermal decay rate was quantified relative to the corresponding zero-temperature decay rate. In the case of a (pseudo)scalar decaying into two distinct (pseudo)scalars, this enhancement is quadratically growing with temperature, see Fig. 8a-c. Considering instead a scalar boson emission off a fermion line, the T -dependence is linear at high temperatures, see Fig. 11. As a practically important example, we notice a relative thermal suppression present in the case of a fermion-antifermion final state emerging from an incoming (pseudo)scalar. The rate of this process was shown to be 0.25 times the zero-temperature quantity at high temperatures. The real-time selfenergy is proportional to a factor of (1 − n F/F ) for each outgoing fermion/antifermion in the same way as in the imaginary-time formalism [14]. Hence, the outgoing particles are accompanied by Pauli suppression through the distribution function n F ; at high temperatures, this distribution approaches 1/2 and the outgoing particles are Pauli blocked by medium particles with this probability as is manifested in Fig. 9.
The one-loop self-energy for (anti)commuting fields had previously been written down in the imaginary-time formalism by Weldon [14]. Ho and Scherrer [34] extracted decay rates for the processes of Secs. VII-VIII using the results of Weldon. This work confirms the thermal decay rate of Ref. [34] in the referenced sections within the real-time formalism. The decay rates have been further extended beyond the results of Ref. [34] to include the chemical potentials, and, more importantly, the masses of particles in the loop have been allowed to take on different values. Hence, a new momentum region 0 ≤ p 2 < (m 1 − m 2 ) 2 has been investigated for a (pseudo)scalar decaying into two (pseudo)scalars, a region in which decays are kinematically forbidden in the zero-temperature theory.
The explicit real-time self-energy of the scalar-scalar loop has been published previously in Ref. [40] and was reproduced in this work. The fermion-fermion self-energy has been discussed in Ref. [14] in the imaginary-time formalism. The calculated decay rate of Fig. 9 reproduces the imaginary-time result of Ref. [34] if the masses of virtual states are taken to be identical. The explicit real-time self-energy of the scalar-fermion self-energy does not exist in the literature to the best of our knowledge. The thermal decay rate has been extracted from such a self-energy, see Fig. 11.
Here, κ has explicitly been extracted from the interaction-term S I of the action to visualise the order of expansion in the perturbation series. To each order in κ it is then possible to calculate the contribution to Z[0] from the vacuum bubbles by computing

Cumulant expansion
Recall the generating functional for the connected Green's functions, W [j] = log Z[j] and the cumulant expansion from Sec. II D 2. Using these expressions, one arrives at an alternative expansion for the generating functional at j = 0 as This is the defining relation for the connected averages (S I ) n con . By combination of Eq. (A.1) and Eq. (A.3), the following identification can be made This can be considered as a perturbative series expansion in κ of the logarithm on the left hand side. The expression on the right-hand side contains, apart from powers of κ, the newly defined connected averages, the cumulants (S I ) n con . These are defined explicitly order by order in κ through the full expansion of the sums in the above expression. In order to obtain the explicit expansion of Eq. (A.4) below, one may exponentiate the left-and right-hand sides and expand the resulting exponential and the two sums to identify the cumulants in each order in κ: (−κ) n+m n!m! (S I ) n con (S I ) m con + . . .
The cumulants can now be identified order by order in κ as: The zeroth-order relation provides nothing but a trivial identity. By making use of the relation at order κ, the second-order relation can be solved for the square-cumulant (S I ) 2 con in terms of ordinary averages: This relation together with the definition at order κ further gives the third order cubed connected term as Collecting Order Cumulant definition κ 0 1 ≡ 1 κ 1 S I con = S I 0 κ 2 (S I ) 2 con = (S I ) 2 0 − S I With this scheme, (S I ) n con to any order n may be defined iteratively. By explicitly calculation of the cumulant at the first four orders (up to κ 3 ), their diagrammatical interpretation is clear.
b. Order κ 1 This calculation may straightforwardly be performed.
c. Order κ 2 At second order, calculations are here shown explicitly.
d. Order κ 3 At third order in κ, the calculation is more cumbersome but not any less straight forward.

Diagrammatically this expression is
It is clear from the diagrammatic calculations above that the cumulants precisely correspond to the fully connected diagrams at each order in κ. This can be shown in general to any order, see for example Refs. [1,13].

Appendix B: An example of a thermal loop calculation in the real-time formalism
In this section, the evaluation of a thermal loop integral is provided as a reference problem. The scenario considered is that of Sec. VII and the diagram is generated by the interaction term The explicit integrals of interest are

SS -integral
The first term in Eq. (B.2) is the cross-term of the non-thermal part of the two propagators and for the loop of interest it is This integral is recognised from the zero-temperature theory and includes no reference to the thermal parameter β. One may in this instance appreciate the separation of the thermal part from the zero-temperature contribution since it is immediately realised that the calculation reduces to the zero-temperature result if T → 0. In this limit, the above integral is the only remaining term since the others contain the thermal distribution function which vanishes when T → 0. The evaluation of this integral is a standard procedure from zero-temperature QFT and the result may be found in any textbook on loop calculations (e.g. [1]). Therefore, the calculation of this integral will not be repeated in any detail here. It is enough to say that the procedure involves the Feynman parametrisation of the two fractions followed by a shift of the integration momentum k. By a subsequent Wick rotation of the energy component according to the integral may be recast into a Euclidean integral on polar form, generally in D dimensions as Generally, integrals on this form may be UV divergent at D = 4 and a scheme for regulating such divergences must be deployed. A common approach is dimensional regularisation techniques in which the dimensionality of the above momentum integral is taken to be D = 4 − 2δ with δ ∈ (0, 1). The innermost integral may be solved in the general case by shifting to polar coördinates and solving the angular and radial contributions separately. The well known result is where the Γ-function parameterises the divergence. Poles arise for D → 2m, m ∈ N such that arguments of the Γ-functions becomes negative whole numbers. The divergence may be explicitly extracted as γ E is the Euler-Mascheroni constant. The above momentum integral is an important relation since many loop integrals may be expressed on this form. Hence, there is seldom any need for explicit momentum integration 19 .
In the case of interest to the evaluation of the integral iI It is important to note a consequence of modifying the dimension D = 4 → 4 − 2δ when regularising the loop integral. The action have dimension 0 when the Lagrangian density is integrated over space-time. When modifying the dimensionality s.t. the space-time integration over the Lagrangian d 4 x → d D x, this requires the missing dimensionality to be accounted for elsewhere. In the interaction term, it is common to extend the dimension of the coupling constant in order to absorb extra dimensions. This is usually done s.t. κ → κμ p whereμ has dimension 1 and is defined as the renormalisation point to some power p. The exponent x is determined through bookkeeping of dimensionality.
Introducing the renormalised coupling constant κμ p squared as a pre-factor to the above integral, one may move the resulting μ 2 into the logarithm. The final step is to evaluate the resulting integral in x over the logarithm which is proportional to the following form In the case of interest, the constants in x are The integral may be evaluated through factorisation of the polynomial as (B.14) All but the last four terms of the evaluation simplify trivially s.t.
The final functional form of this expression depends on the sign of b 2 − 4ac which is determined by the sign and magnitude of p 2 . Disregarding corrections to order , this expression is positive when: and negative when In the negative case when p 2 lies between the difference squared and the sum squared of the two masses, a further simplification can be made.
where the three quantities x ± and C was defined as In the positive case, p 2 is either above or below a certain threshold value: Each of the four logarithmic terms may acquirer a negative argument thus making it undefined. A convenient approach is the identification, and further expansion, of, the following square roots: An identity for the principal value of complex logarithms ln [x ± i ] = ln |x| ± iπΘ(−x) + O( ) may be used in order to expand the logarithms as: Collecting the above result, the full non-thermal part of the loop integral of interest to this subsection may be written down: The last integral is (B.28)

Isolating the imaginary part
In the case of decay rates, the main point of interest is the imaginary part of the loop diagram. In the above case of the non-thermal part, it proved to be a decent amount of work for essentially a contribution of 2πi and it should be noted that a quicker calculation is viable in order to extract the imaginary part of the loop.
After dimensional regularisation in the previous subsection, the logarithm to be integrated over the Feynman parametrisation may be written on explicit complex form using the identity preceding Eq. (B.23). Thus, the imaginary part of the loop may be extracted at an early stage: Recall the definition of a, b (Eq. (B.13)) and note that after expansion, the i contribution has been taken into account. Therefore, c now reduces to c with The only remaining task is to identify for what regions of p 2 the roots of the polynomial in the Θ-function fall inside the integration domain and for which regions the entire interval is to be considered. Factorise To exhaust the possibilities of sign combinations, the values of the roots has been tabulated, see Tab. I.
In the case 1 in the table, one may show that In the case 2 in the table, one may show that Λ + ≤ Λ − and Hence, only for p 2 ≥ (m 1 + m 2 ) 2 is the polynomial negative inside the region of integration giving a finite imaginary component. Both roots Λ − , Λ + ∈ [0, 1], ordered as Λ + ≤ Λ − . The result of the previous subsection for the imaginary part is the swiftly recovered: This section provides an example calculation of the so-called mixed first-order contribution to the thermal self-energy. The mixed integrals are due to the cross-terms of the two propagators that mix the thermal and non-thermal propagator terms.
Since the two integrals differ in structure only by a substitution of variables (q = p − k) the second integral has simply been denoted as (1 ↔ 2). This notation implies a naïve interchange of indices 1 and 2. The full mixed contribution is The integrals are evaluated through the Sochocki-Plemelj theorem 20 which requires a smooth closed simple curve C in the complex plane. Splitting the Feynman propagator and the δ-functions results is Thus, the real and imaginary parts of the loop has been separated and are written out below for clarity. The real part: and the imaginary part of the mixed contribution to the loop: a. The imaginary part The limits k , , + , k , , − are given by restrictions from the δ-functions on k listed above. Now, the integral of interest may be recast on this polar form: where the function f (rs) (x) was defined as with r, s = ±1. Now, simplify the δ-function using where the roots x root are points x 0 so that f (rs) (x 0 ) = 0. Solving for x 0 gives The end points of the interval are of primary interest as to determine the limits of integration and one may consider specifically the values of |k| that fulfill the equalities, labelled |k| 0 . Separately, the two equalities may then be squared and stated as These resulting limits on the polar integral for the imaginary part of iF (2) (p; m 1 , m 2 ) + (1 ↔ 2), indicated above as k , , + , k , , − . The limits are p 2 -dependent and contributions may come from any of the three δ-functions in the integral over cos θ. In Tab. II, the different regions are presented.
TABLE II: Integration limits on the radial variable |k| are displayed below. The total integral of interest contains three terms, each proportional to a δ-function of different arguments. The δ-functions give rise to integration limits on the radial coordinate and labelled accordingly as k , , ± to distinguish the three contributions. The integration limits for 0 ≤ p 2 ≤ (m 2 − m 1 ) 2 as m 2 < m 1 and the limits in parenthesis for p 2 < 0 as m 2 > m 1 correspond to the case when 1 + < 0. Empty entries for the limits imply no contribution from the δ-function in this region.

(B.71)
Note that the right-hand side is purely imaginary. The fourth term above proportional to the integrand δ(p 0 + ω 2,k + ω 1,p−k ) may be discarded due to the physical condition p 0 ≥ 0.
If m 2 > m 1 , the contribution come from δ p 0 − ω 2,k + ω 1,p−k while for m 2 < m 1 the contribution come from δ p 0 + ω 2,k − ω 1,p−k . The integral becomes Secondly, consider the case of r = +1. In this case, one should take the relevant δ-function as δ p 0 + ω 2,k − ω 1,p−k , the condition imposed by this function is βp 0 + y > 0. The integration proceed very similar as to the first case in this momentum region using: For the momentum region where p 2 ≤ 0, contributions come from both δ p 0 − ω 2,k + ω 1,p−k and δ p 0 + ω 2,k − ω 1,p−k . The imposed limits of integration on |k| are different in the two terms.
First consider a mass hierarchy of m 2 > m 1 . The limits on |k| are different in the cases of 1 + > 0 and hence, the integration result is identical.

Summary of the results
Collecting the results of the thermal loop integration presented in Secs. B 1-B 4 for all different momentum regions are presented in his subsection.
The non-thermal contribution comes from the term iI Here, the last integral term is (B.85) with The mixed thermal-non-thermal integration gives the contributions   This section presented, in detail, the type of integrals encountered in the evaluation of thermal loop diagrams by making use of the simplest possible loop of two scalar particles. It is clear from the result that the final expressions simplify further if the two loop masses are identical. Additional complications for the procedure of integration arise when the propagators carry Lorentz structure. It should be noted, however, that both such complications and loop divergences are introduced mainly through the non-thermal integration and, therefore, the thermal terms does provide any major increase in complexity.