Two-loop thermal spectral functions with general kinematics

Spectral functions at finite temperature and two-loop order are investigated, for a medium consisting of massless particles. We consider them in the timelike and spacelike domains, allowing the propagating particles to be any valid combination of bosons and fermions. Divergences (if present) are analytically derived and set aside for the remaining finite part to be calculated numerically. To illustrate the utility of these 'master' functions, we consider transverse and longitudinal parts of the QCD vector channel spectral function.


INTRODUCTION
In a relativistic plasma, the rates of processes like particle production and damping are derivable from the imaginary part of a particle's self-energy [1,2]. That quantity, also called the spectral function, depends on the energy k 0 and momentum k which can occur only in the combination K 2 ≡ k 2 0 −k 2 at zero-temperature. This is not so for thermal systems, where the medium's rest frame is distinguished and the temperature T = 0 joins k 0 and k = |k| as an important scale in the problem. Introducing another scale can dramatically alter the naive weak coupling expansion: New infrared singularities foreshadow that next-to-leading order (NLO) corrections are large, or even that resummation is obligatory.
One such instance is the photon spectral function in hot QCD [3][4][5]. Truncating the perturbative result for the self-energy to order e 2 in the electromagnetic interactions, we denote by Π (l) the ensuing contribution from g 2l to the strong coupling expansion. The series then takes the form with a supposed ordering by powers of g 2 . For a strict loop expansion, the 'coefficients' of g 2l are themselves functions of k 0 and k but independent of g. However their dependence on the external momentum K can (and does) spoil this power counting, e. g. when |K 2 | ∼ < g 2 T 2 .
For many observables only leading-order (LO) or partial NLO results are known, making it unclear where (1.1) actually breaks down. To obtain an approximation that is justified for all k 0 , the fixed order expansion can be 'matched' with the resummed approach (which works near the light cone). That was the idea put forward in Ref. [11], where it was tested for Im Π µ µ with k 0 > k . Here we also consider energies below the light cone and separately the polarisation state Im Π 00 (1) , as inspired by Ref. [12]. Our goal is to assist in the effort of quantifying another order in perturbation theory by cataloging a general class of two-loop spectral functions. (One of the earliest attempts in this spirit provided the first correction to the gluon plasma frequency [13].) What follows is rather technical, but lays out a generic approach to evaluate those integrals frequently needed in NLO computations. All code used for determining the finite thermal parts (defined as specified below) is supplied in Ref. [14]. The primary task of that code is a phase space integration of amplitudes squared, with thermal weightings appropriate to each process.
To compute loop integrals at finite temperature, we apply the imaginary time formalism for massless particles. Free scalar propagators, carrying either bosonic (s = +1) or fermionic (s = −1) momentum, are denoted by ∆ s (P ) = 1 (1. 2) The integer n specifies the Matsubara frequencies and Θ is the Heaviside step function.
We regularise the spatial momentum in d = 3 − 2 dimensions with the modified minimal subtraction (MS) scheme and renormalisation scaleμ. The trace over momentum P = (p 0 , p) at finite temperature is defined by where γ is Euler's constant. We follow [15] to carry out the sums over p 0 , defined in (1.2).
This paper is organised as follows. In section 2 a general class of master sum-integrals is introduced and those considered here are specified. They are then evaluated, one by one, in sections 3, 4, 5 and 6. (For completeness, and as an important cross-check on our results, the K 2 T 2 behaviour of each sum-integral is derived analytically in Appendix D.) Finally, we 'sum up' in section 7 and mention some potential applications.

LIST OF INTEGRALS
Let us define, for generic sum-integrals I as functions of the external four-momentum K = (k 0 , k) , a uniform notation (for m = n = 0, cf. [16]) where ∆ i ≡ ∆ s i (P i ) . The case n = 0 is abbreviated by I (m) abcde . Together with P and Q, the integration variables, K determines all the propagating momenta (as depicted in Fig. 1), P 1 ≡ P , P 2 ≡ Q , P 3 ≡ R = K − P − Q , Thus we summarise the statistical content of (2.1) by (s 0 , s 1 , s 2 ) = (±, ±, ±) , not including it explicitly on the notation. At T = 0 , the statistics play no role and these integrals can be evaluated using well-known methods [17]. The vacuum contributions will dominate over the thermal ones for large K 2 .
Relative corrections are suppressed by powers of K 2 that can be formally organised with an operator product expansion (OPE) [18]. Thermal effects are important for K 2 being of similar order to T 2 , the regime we consider here to calculate the imaginary part of (2.1) 1 .
The remaining limit, K 2 T 2 , is frequently associated with the need to resum all orders in perturbation theory to cure a diverging spectral function. With that in mind, we shall often also discuss the master integrals for K 2 ≈ 0 . 1 By evaluating it at an energy k 0 + i0 + .

Strategy: m, n ≥ 1
To take care of the powers of the energy in the numerator of (2.1), namely p m 0 and q n 0 with positive integers m and n , we employ the following strategy. In special cases, the corresponding graph has a symmetry in momenta P and Q (i. e. if a = b and d = e), and one can take advantage of the same symmetry for the powers of p 0 and q 0 . But in general, one should make use of the Fourier representation of the (massive) scalar propagator where E p = λ 2 + p 2 and λ is the particle mass. Here we also introduced f − s = sn s , f + s = 1 + sn s and the distribution function n s (E) = [exp(E/T ) − s] −1 .
Beginning with the case m = 1 and differentiating with respect to τ under the Fourier transformation, one effectively multiplies 2 by the conjugate variable, This is inserted into (2.1) before carrying out the frequency sum. It is straightforward to differentiate (2.2) with respect to τ . Hence (2.3) provides an extra factor of (vE p ), counting the minus sign from integrating by parts.
The case m = 2 is also elementary from the relation p 2 0 ∆ s (P ) = 1 + E 2 p ∆ s (P ) . (2.4) By applying (2.4) and (2.3) in sequence one can reduce p m 0 for m ≥ 2 in (2.1) to a sum of powers of (vE p ) with simpler masters. (And the same strategy works for q n 0 .) The benefit of all this, is that the frequency sums are relatable to cases with m = n = 0 ; any complications will move to the integration over the three momenta p and q that follows. 2 Here we generalise (1.2) to have a mass: ∆ s (P ) = (p 2 0 − p 2 − λ 2 ) −1 . This will help later on, as an infrared regulator.

Example: A QCD spectral function
Integrals of the form (2.1) [with statistics (+, −, −)] can be used to express the NLO photon self-energy in an equilibrated QCD plasma at zero chemical potential. The emission rate is derived from the (contracted) spectral function Im Π µ µ [19,20], but here we also study Π 00 . Due to the Ward identity at non-zero temperature, the polarisation tensor has two independent components. identified with the longitudinal and transverse polarisations: The difference between Π L and Π T is purely thermal, at zero temperature there is none [12].
Accordingly, Π µ µ and Π 00 are enough to completely specify Π µν at finite temperature. Denoting the number of colours by N and the group factor by C F ≡ (N 2 − 1)/(2N ) , they 11020 − I (2) 10120 As part of the procedure to reduce Π 00 and Π µ µ to a minimal set of integrals, we removed angular variables in the numerator thanks to relations like p · k = p 0 k 0 : : (1,1) : : : We have grouped the master integrals into classes (designated by the numeral on the graph), according to the associated topology. The topology of the class does not necessarily follow directly from the assignment of loop momenta in Fig. 1. A change of integration variables is sometimes required to relate them. The first three classes are all presented together in Sec. 3 because they consist of simpler one-loop subgraphs that factorise. Classes IV, V and VI can be considered genuinely two-loop and will receive the most attention, being discussed in Secs. 4, 5 and 6 respectively.
Before moving on, a brief comment on one-loop diagrams is in order. They have been studied extensively in the literature and are usually considered in the hard thermal loop (HTL) approximation. (Higher order HTL results have also been investigated, cf. Ref. [21].) This is common in the high temperature limit [15] because it affords analytic expressions for the self energy, and supplies results that are automatically gauge invariant. If one relaxes the HTL assumption that the external momentum K 2 is much smaller than T 2 , the complete self-energies must be evaluated numerically [22]. Effective field theory methods have also been developed recently to compute associated power corrections [23]. The imaginary parts involve phase space integrals for 1 ↔ 2 'decays' which are needed for our master diagrams II and III as well as certain terms arising in V and VI. Appendix A gives details on the integration measure, where we also discuss the 2 ↔ 2 and 1 ↔ 3 processes to be utilised when more intermediate states can go on-shell.

DIAGRAMS I-III (FACTORISABLE TOPOLOGIES)
Those diagrams we assigned to classes I, II and III are reducible to products of simpler one-loop integrals. They can all be recast as those having c = 0 in (2.1) 3 , which implies that P and Q dependence of the integrand does not mix. It is thus useful to recap a general one-loop function, defined by The frequency sum over p 0 is well known [15], and we organise the subsequent integration over spatial momentum p according to Appendix A .
The cases where b = 0 are local contributions. For the one with a = 1 we abbreviate the integral by where ζ is the Riemann zeta function [24]. Type I self-energies are then constant and we need not discuss them because they have no imaginary part. Moreover, since I n is zero in vacuum (for m ≥ 0), the type II integrals are entirely thermal corrections.
For a = b = 1 , integrals of the form (3.1) are usually considered in the limit k 0 ∼ k for which the HTL functions can be used. But in general, the emerging integral expressions must be evaluated numerically [22]. Only when taking the imaginary part, thus putting internal momenta on-shell, is the integral doable analytically.
Let us introduce three useful functions F, G and H that make the statistics explicit, (The same spectral functions in Ref. [25] were labelled by a 'd' and 'g' respectively.) Since I n was given above, we now turn to the K-dependence of F m , for the particular cases needed. As derived in Appendix B (and applicable for both k 0 > k and k 0 < k) where the distribution function n s was defined below (2.2) and is evaluated at light cone momenta k ± = (k 0 ± k)/2 . We abbreviated the quantity s 0 n s 0 (k 0 ) by n 0 .
where the summation extends over v 1,2 = ±1 [a definition of f v s is given below Eq. (2. 2)]. The strategy discussed in Sec. 2 was applied to cover the cases m ≥ 0 . Dimensional regularisation is adopted because the related function G m will include a customary ultraviolet divergence: Take G 0 and H 0 for instance, which are real and imaginary parts (respectively) of the same function. Their zero temperature limits can be read from Branches of the logarithm are made explicit; we write log X to mean log |X| . Thus G 0 bears an ultraviolet divergence and so O( ) terms must be kept in H 0 when multiplying them together. For that reason we write to make the dependence on the scaleμ explicit. Setting d = 3 in (3.5) allows us to find H Thus H where the last term is from a symmetry as specified in (3.3). The first line above includes all divergences and yields the entire result for T = 0 . The second line is purely a finite thermal function which is not present in vacuum. Note that the divergent first line is also a function of the temperature and we omit it in Fig. 3, where the master integral is displayed.
On the light cone, a logarithmic singularity in this thermal part may arise from ψ . This divergence is softened by the extra weight K 2 /T 2 that was used in Fig. 3, and the plotted function is zero at k 0 = k . We note that if s 0 = 1 (implying s 1 = s 4 and s 2 = s 5 ), some simplifying relations hold

DIAGRAM IV (SETTING SUN)
We now consider the first genuine two-loop structure, specifically the integral I 11100 , which gives e. g. the first non-zero contribution to the imaginary part of the self energy in a scalar ϕ 4 -theory [26]. Another master of the same class, I (0) 01110 = 0 , is identically zero due to integration by parts identities [24]. The 'setting sun' graph is given in vacuum by where κ was introduced in Eq. (3.6). This vacuum result has an imaginary part for K 2 > 0, associated with the threshold for massless particle production. In a thermal medium, the Landau-damping mechanism explains why this imaginary part also builds up below the light cone [27]. Explicitly, Im where f v s was defined just below (2.2) and takes the arguments at energies E 1 = |p|, E 2 = |q| and E 3 = |r| which are on-shell. This integral is labelled 'f' in Ref. [25].
Let us clarify the physical content of Eq. (4.2). The sum over the signs v i = ±1 enumerates eight distinct physical interactions, with external momentum K = (k 0 , k) . We denote the corresponding fields φ i for argument's sake with s i = +1. As an example, the term with v 1 = v 2 = v 3 = +1 represents the probability for decay φ 0 → φ 1 φ 2 φ 3 , with a statistical weight of (1 + n B )(1 + n B )(1 + n B ) for spontaneous emission, minus the probability for creation φ 1 φ 2 φ 3 → φ 0 , with a weight n B n B n B for absorption. There are many other processes, such as φ 0 φ 2 φ 3 → φ 1 minus φ 1 → φ 0 φ 2 φ 3 and so on [1]. Equation (4.2) may be simplified into a two-dimensional integral (now for general s i ) Im I (0) 11100 = n −1 0 (4π) 3 dp dq W IV (p, q) n 1 n 2 n 3 , where we abbreviated n i = s i n s i and agree that the arguments 4 of the distribution functions may be negative. The 'kernel' W IV (defined below) also depends on k 0 and k , but not on the temperature. The momentum moduli p and q have been generalised to negative values, which implicitly incorporates the sum over the signs {v i }. And the statistical weight has accordingly been re-expressed using An explanation that starts with Eq. (4.2) is given in appendix C, where we also show how to calculate W IV from kinematic constraints. Here we simply state the result: where k ± = (k 0 ± k)/2 are the light cone momenta.
Of note is that W IV (p, q) p/k for p → 0, which suppresses the log divergence from n B (p) . Furthermore, W IV = 0 in regions that are kinematically forbidden, providing limits on the p and q integrals. Continuity of Im I 11100 at k 0 = k follows from the very same property in (4.5). We note that the limit k → 0 is also well-defined and leads to W IV = sgn( pq(k 0 − p − q) ) where it has non-zero support.
The statistical factor (4.4) includes the vacuum contribution for v 1 i.e. where p and q are positive and p + q < k 0 . It is the leading term in Eq. (4.4), after expanding in combinations of the distribution functions In Fig. 4, the energy dependence of Im I (0) 11100 is shown for k = {0, 1 10 , 1, 10} × T . Here the vacuum result (4.1) was subtracted, i. e. we actually plot The exception, for massless particles, occurs when k = 0 so that there is an essential singularity at k 0 = 0 [28]. Hence the zero momentum curve in Fig. 4 is finite for k 0 → 0 and not equal to the same limit at fixed |k| > 0 .

Kinematics
The region(s) where Eq. will reveal why it holds in general for the real corrections. To illustrate, we first consider the simpler case k → 0 and then explain what happens when k < k 0 and k > k 0 separately.
The function W IV is not Lorentz invariant -if it were, we could perform the whole calculation in the rest frame. Nevertheless, specialising to k = 0 will be useful as a starting point. For example, in the channel with v 1 = v 2 = v 3 = +1 we require vectors p and q that satisfy Clearly p + q ≤ k 0 in general, with equality if and only if p = −q and |p| = 1 2 k 0 . Moreover, by the triangle inequality we have |p − q| < |p + q| < p + q (here p, q > 0). The upper limit thus gives 1 2 k 0 < p + q, while the lower limit leads to p < 1 2 k 0 if p > q and q < 1 2 k 0 if p < q. That defines the relevant domain in the (p, q)-plane; see '1' in Fig. 5, where W IV = +1.
This time p + q ≥ k 0 and the very same triangle inequalities give p > 1 2 k 0 if q > p and q > 1 2 k 0 if q < p. Similar results hold if v 3 = +1 and if exactly one of v 1 , v 2 is equal to −1. However if two or more of {v 1 , v 2 , v 3 } are negative, the equivalent of (4.7) cannot be satisfied (for k 0 > 0). Hence in the sum over {v i }, only half of the summands contribute. This is summarised by the wedge-shaped regions '2', '3' and '4' in Fig. 5. In each of these three regions, W IV = −1. Although they include arbitrary large momenta (in absolute value), those much larger than the temperature are cut off by the thermal distribution functions.
We now consider k > 0, but still less than k 0 so that k − = 1 2 (k 0 − k) is positive. The external vector k now plays a role, i. e. (4.7) is supplanted by (4.9) Of course p + q ≤ k 0 still holds, but now equality can occur for all q ∈ [k − , k + ]. The triangle inequality gives k + |p + q| > |k − p − q|, and therefore Hence the lower bound on p + q is diminished to k − . Similarly, the other side of the triangle inequality gives |k − p − q| > |p + q| − k. That produces k + ≥ max(p, q), which is higher than the upper bound for k = 0. Generalising to other channels is trivial; see As k is increased beyond k 0 , the virtuality K 2 = 4k + k − becomes negative. The channel with v 1 = v 2 = v 3 = +1 ceases to be accessible; conservation of energy (4.9) cannot be  New wedges open up in the (p, q)-plane -they correspond to having exactly two of Fig. 7 they are labelled '5', '6' and '7'. Regions '2', '3' and '4' are carried over from the case k < k 0 , with modified boundaries: For example in '2', the same triangle inequality as before gives |q| < −k − . Region '5' has v 1 = v 2 = −1 and v 3 = +1 , so in some sense it is the intersection of '2' and '3'. In '5' we have p and q negative such that p + q < k − . (And there are similar constraints for '6' and '7'.) We note that adjacent regions are separated by lines where p, q or r is zero. These boundaries are important because they mark the location of potential singularities coming from bosonic distribution functions.
The expression in (4.5) can be simplified in each of the demarcated regions discussed for k 0 > k and k 0 < k . In the former case, we recover Eqs. (33)(34)(35)(36)(37)(38)(39)(40)(41)(42)(43) As before, now with k > k 0 so that k − < 0. (The blue hatched region is forbidden by kinematics. It is, incidently, also where the vacuum contribution comes from.) Ref. [25] (with an adjustment in variables, namely p → k 0 − p). For the more complicated master integrals still to be studied, the same regions must be considered, though the kernel functions will be different. Hence, the foregoing analysis serves to spell out what must be included for the task of numerical integration.

DIAGRAM V (THE SQUINT)
In the previous section there is only one 'contribution' to the imaginary part: all of the three internal energies are set to their on-shell value giving (4.2). But more generally, the discontinuity may always be written in terms of products of two amplitudes that are separated by on-shell 'cut' propagators [1]. For classes V and VI, there is more than one way to do this; see Fig. 8. These contributions, which may be readily identified after carrying out the Matsubara sums, are separately infrared divergent. They differ by one extra loop momentum being on-shell in the real case, yielding a tree-level decay that may be treated in a similar manner to that of the previous section. The virtual correction has an internal loop due to one fewer final state particle than before and includes a two-body phase space integration. The latter is also ultraviolet divergent, seen in the vacuum result which, for K 2 positive, contains an imaginary part ∝ 1/ . This divergence will then acquire a temperature dependence in the medium. When assembled together, these infinite parts cancel in actual observables. The two terms shown in Fig. 8 also have separate collinear divergences which exactly compensate only in their sum. This cancellation is somewhat intricate at finite temperature and affirms the Kinoshita-Lee-Nauenberg (KLN) theorem, in this case applied to individual graphs [29,30]. The amplitude (a) in Fig. 8 may contain a large 'eikonal factor' if the denominator of an internal line is zero, where θ kp is the angle between p and k. Such configurations exhibit two collinear outgoing A fictitious mass λ (for the momentum R) is convenient to use as a calculational tool [17] 5 . It prevents L 2 → 0 in (5.2) and proves that the whole expression is rendered finite as λ → 0 . [See ahead (6.3) for the abbreviations n i .] The first term above represents the real correction, and involves the thermal weight (4.4) from before. Virtual corrections lead to the second term and depend on the mass λ so as to compensate for the contingent singularity in the first. Equation (5.3) is valid for m, n ≤ 1, otherwise there could be more terms.
We may obtain W V (p, q ; λ) by the techniques laid out in Appendix C. A compact formula for it can be given with the help of some funny notation: To also indicate whether p is in the interval Ω = [k − , k + ], we define Thus the complementary set Ω c = R\Ω contains p if the value ofō p = 1 − o p is equal to unity. Using o p andō p will signal terms that are included or not, depending on the relative ordering of k , k 0 and p . With these definitions, and regulating through λ → 0, the weight function can be expressed as 6 where we have introduced the ratios This implies that it is too soon to set λ = 0 in those domains of the (p, q)-plane and brings us to incorporate the missing virtual pieces.
Again deferring details to Appendix C, we write, for λ → 0, This reveals how the anticipated ultraviolet divergence in (5.1) emerges from the loop in (b) of Fig. 8. At the same time, it can be seen that the λ-dependence in (5.3) cancels for λ = 0, with the logarithmic mass singularities compensating perfectly, and in exactly the domains where (5.2) is satisfied. They coalesce in this way both above and below the light cone.
For the purpose of plotting, we subtract a piece that is ultraviolet divergent from (5.3) and coincides with its vacuum result for T = 0. What remains is thus finite and proportional to T 2 for large photon virtualities, see Fig. 9. The function is continuous across the light cone unless it diverges there, in which instance the singularity is the same if k 0 → k is approached from either above or below. That is why, in Fig. 9, we multiply the whole function by |K 2 | which is enough to render the blow-up finite on the light cone. (It also gives the whole master integral a dimension of T 2 .) To be clear, and following [25] (in which this master is labelled 'h'), the part subtracted is unaffected and hence the λ-dependence (and ultimate lack thereof) is the same as before.
For plotting in Fig. 9, we subtract from the result. The shape is similar to Im I 11110 , but ratios of the large-k 0 limit to its value on the light cone are different.
For m = 0 and n = 1 , the virtual parts need to be reconsidered. The result, with details available in Appendix C, can be written, for λ → 0 , Figure 10 displays the result as a function of positive invariant mass M = k 2 0 − k 2 . We have followed Refs. [16,25] by using M to define as a proxy for the average three momentum squared. (K ν are modified Bessel functions.) This quantity assumes Boltzmann distribution functions to average k 2 for a fixed mass.
The curves shown in Fig. 10 use k → k 2 ave (M ) as a rather crude substitute for 'typical' momenta. As M → 0 , we note the clear log-divergent behaviour of the master integral for the statistics shown in the figure.
The divergent piece that was subtracted is given by [note the ordering of s 4 and s 1 on We have also checked numerically that Im I (0,1) 11110 provided s 2 = s 3 . This relation follows by a shift of integration variables, but it seems difficult to discern this rule by simply looking at the explicit form of the integrand.
In Sec. 2 we listed among the type V integrals, one with a propagator of negative power for K 2 T 2 [3][4][5]. This is the singularity alluded to in the Introduction, which mandates screening effects to be incorporated through resummation. The log-singularity on the light cone can be traced back to the integral we are about to discuss.
Let us consider the particular combination of master integrals defined by (This one is labelled by 'h ' in Ref. [25].) After carrying out the Matsubara sum, one arrives at Eq. (C.2) in the Appendix. Along the same lines as (5.3), this can be expressed by where the two terms are the real and virtual parts respectively. The first weight includes the same manner of λ-dependence as that previously defined in (5.4), which allows the previous argument to be partially recycled here. To explicitly define W , let Accordingly, in the first summand of (5.12), we have The arguments of g are defined in Appendix C; see Eq. (C.9). Virtual corrections in (5.12) require the function We see that the auxiliary mass λ enters only in the previously defined functions, W V and U V , so that the pattern of cancellation is unchanged and allows us to set λ = 0 . The divergent piece that we choose to subtract is We have checked numerically that Im I 11110 = 1 2 K 2 Im I  This makes the square brackets in (5.14) a difference between two thermal moments; see (D.20), for the particular statistics s 2 and s 3 . If they are the same, it gives zero. A noteworthy feature is the behaviour for K 2 → 0 , where the functions W and U simplify.
There is a discontinuity due to the latter, cf. Ref. [31], defined by the function's limit k 0 → k + 0 + minus k 0 → k − 0 + . For a general statistical configuration, it is given by where the integration is meant in the principal valued sense.

DIAGRAM VI (CAT'S EYE)
The most complicated form of (2.1) that we consider here is a = b = c = d = e = 1 , which requires a careful cancellation of real and virtual diagrams. We can deploy the same strategy as in Sec. 5, albeit now with two real and two virtual amplitudes. (Each is related by the symmetry s 1 ↔ s 4 and s 2 ↔ s 5 .) An auxiliary mass λ is again attached to E 3 = √ r 2 + λ 2 so that collinear singularities can be regulated.
In Ref. [16], the function Im I is finite and has no imaginary part. The leading contribution to Im I The weight function W VI that is needed in (6.2), carrying over some notation from (5.13), reads for k 0 > k and λ → 0 , W VI (p, q; λ) = 1 4kK 2 r log W VI +ō p log W VI +ō q log W VI (6.4) Below the light cone (k 0 < k), the function should instead be In Eq. (6.4) the following ratios were defined: We note the appearance of a log-divergence, just where expected and signaled by the coefficients o p and o q . The formulas for W VI are symmetric in arguments p and q , as is (then) the other real correction, which comes from p → k 0 − p and q → k 0 − q .
The case m = 2 and n = 0 can also be written in the form of (6.2). One way of seeing this, is to follow (2.4) and rewrite where the explicit integrands for the two terms can be found in Appendix C. (After s 1 is interchanged with s 5 in the first term above.) Equation (6.2) can be recovered after some manipulations of integration variables.
The virtual corrections are triangle diagrams, one of which is given by the second line of (C.3) in the Appendices. Their calculation is similar to those studied in the previous section and is given explicitly in Appendix C. Taking up the first term in the second line of Eq. (6.2) (the other term follows analogously), for λ → 0 it can be expressed by Together, the last two terms in (6.2) are seen to combine with the real corrections [see (6.4)] so that the complete expression is λ-independent.
One needs to be careful for m = 2 because a genuine ultraviolet divergence must be subtracted. That divergence is temperature dependent and originates from the virtual contribution U VI (q) . It is equal to which we subtract for plotting purposes. (For the related master integral with m = 0 and n = 2 , one replaces s 2 → s 1 and s 5 → s 4 in the above.) Figure 14 depicts Im I An ultraviolet divergent part was subtracted, see Eq. (6.7). and the energy k 0 = M 2 + k 2 ave .

RESULTS AND CONCLUSIONS
We have computed a list of spectral functions that originate from self-energy diagrams with two loops, generalising the results of Refs. [16,25] to below the light cone and considering a larger set of master integrals with m, n > 0 . In so doing, we have separated the ultraviolet divergence where appropriate and shown how to determine the finite remainder numerically. Validity of the KLN theorem was explicitly demonstrated in Secs. 5 and 6, by careful analysis of the collinear phase space. We considered any arrangement of propagating bosons or fermions allowed by the diagram's topology. The code used for our numerical evaluation is publicly available at Ref. [14].
Returning at last to the QCD corrections for the photon spectral function, which was our original motivation, the imaginary parts of Eqs. (2.6) and (2.7) can now be evaluated. Firstly, all the temperature dependent divergences that were individually isolated end up cancelling and the vacuum NLO result is recovered. Since the thermal parts carry no ultraviolet difficulties, zero-temperature counterterms suffice for renormalisation. Some integrals give zero because they have no k 0 -dependence prior to taking the imaginary part, specifically Im I Other terms in (2.6) and (2.7) also do not contribute, in particular those proportional to without a compensating divergence in the master integral. It is important to keep some of these terms so that the vacuum result is recovered, but Im I At finite temperature the spectral function for the current-current correlator is specified by two scalar functions, ρ T,L = Im[Π T,L ] according to (2.5). Figure 15 shows the energy dependence of the NLO spectral functions for several momenta. The behaviour in (5.10) prevails near the light cone for the transverse polarisation, while the factor of K 2 is enough to ensure that ρ L is zero there. Both functions approaches zero for k 0 → 0 , as they should.
The large K 2 behaviour (of the NLO parts) can be found in Eq. (D.21) of the Appendix.
In heavy-ion collisions [19,20], the observable photon and dilepton rates are proportional to −Im Π µ µ = 2ρ T + ρ L for K 2 = 0 and K 2 ≥ 4m 2 l respectively (m l is the lepton's mass), where perturbative studies have hitherto focused 7 . The complementary region, K 2 < 0 , is not merely academic: It provides an opportunity to cross-check the weak coupling framework, e. g. with non-perturbative Euclidean correlators (of either polarisation) provided by lattice QCD [12].
As noted in the Introduction, our results are relevant for deploying (fixed-order) perturbation theory in contexts where simplifying kinematic assumptions are not justified. It is nevertheless worthwhile, e. g. if analytic expressions are available, to check that these numerical results reproduce the correct behaviour in those limits. We have done this using the 7 Including the special case m l → 0 . OPEs in Appendix D for K 2 T 2 in each integral studied. Our results are also compatible with recent HTL self-energies at NLO, in the limit k 0 , k T [32].
It is worth recalling that the individual masters are not (usually) themselves physical, rather they supply a convenient 'basis' from which observables can be built. In this appendix, we discuss two important phase space integrals which were used in the main text. The following results are valid both above and below the lightcone: k 0 > < |k| .
Without loss of generality, we assume k 0 is positive. Energy conservation is imposed by cutting the diagram, handled with the replacement A.1. Two-particle decay The production rate due to binary encounters, e.g. qq → γ * (the asterisk indicates a virtual photon), is given by a phase-space integration over the momenta of interacting partons. For a given external momentum K = (k 0 , k), we simplify the integration measure with the (on-shell) energies E 1 = |p|, E 2 = |q| . The v i = ±1 are summed over, for i = {1, 2}, so that in (A.1) each channel is represented. Here we also used dimensional regularisation for the terms O( ) needed later. Integrating over d d q is trivial; momentum conservation fixes q = v 2 (k − v 1 p) . The on-shell energy is therefore where θ kp is the angle between k and p .
The external vector k distinguishes an orientation with which to organise the remaining d d p-integration. We choose to align the p z axis with k, so that the azimuthal integration is also trivial. And instead of a polar integration over θ kp , we integrate over the magnitude q = E 2 given by (A.2). The angular limits accordingly translate into a kinematic restriction; |k − v 1 p| < q < |k + v 1 p| (see Fig. 16). Let us therefore express the angular integration bȳ All the necessary scalar products can be formed in terms of these variables. When the internal lines are on-shell, energy conservation implies that k 0 = p + q allowing the q integration in (A.3) to be carried out. Therefore the final result is supported only along dashed line shown in Fig. 16, and reads [1,2] = 1 8π k dp 1 + logμ The p-integral was flagged with a prime as it depends on whether K 2 is positive or negative.
Precisely one of v 1 and v 2 equals −1 in the latter case.
Altogether, one has The discontinuity at k 0 = k is thus given by an integral over all real p , in the principal valued sense. In the rest frame, k = 0 , (A.4) has the net effect of setting p = k 0 /2 .
As an application of the foregoing discussion, define and consider it for ν = {0, 1, 2} in particular. This function, which is needed for the one-loop discontinuity, was also used in the virtual parts of our two-loop calculations.
Letting the exponentials be abbreviated by we find that These quantities are discontinuous on the light cone, explicitly with .

A.2. Three-particle decay
It is sufficient to work in d = 3 dimensions for phase space integrations of the real two-loop corrections (those where three propagating particles are put on-shell). Extending (A.1), we define [1,2,3] ≡ v p,q,r with the (on-shell) energies E 1 = |p|, E 2 = |q| and E 3 = √ r 2 + λ 2 . (The regulator λ is necessary to observe the cancelling divergence between real and virtual corrections.) The v i = ±1 are summed over, for i = {1, 2, 3} . One of the integrals may be simplified using energy and momentum conservation; we choose the r-integral, and complete it by writing

Hence, by fixing
Let us now specify a coordinate system in order to proceed. We choose, following Ref. [33], to align the z-axis along the direction of k and orient the zy-plane to contain p, viz.
The integral over the azimuthal angle φ can then be performed in (A.11). To do so, express the argument of the δ-function by Here we have abbreviated x i = cos θ i . Consider then the integral of a function g(φ) , where h = B 2 − A 2 and the angles φ ± = π ± acos(A/B). For our purposes, the function g depends on φ via cos φ and hence ± g(φ ± ) → 2g(φ ± ) .
We have accomplished five of the nine integrals in (A.10), and are left with which cannot be simplified further in general. The combination v 1 p and v 2 q in (A.13) once again suggests that we formally extend the magnitudes p and q to negative values. The factor of Θ(r 0 ) in (A.11) has been dropped because r 0 = v 3 (k 0 − v 1 p − v 3 q) and so exactly one term in the sum over v 3 = ±1 will contribute. (The relevant value for v 3 is determined for a given p and q -extended to take negative values.) And then, as before, we can just forget the sum over {v i }.
The h-function of (A.14) is quadratic in each of its arguments p, x 1 , q, and x 2 . Requiring h ≥ 0 (equivalent to taking the real part of √ h) summarises the allowable phase space. The formula for F 0 was already given in the main text: see (3.4). Here we derive that result after discussing some simple properties. The expansion of the function F 0 about the light cone energy, k 0 = k , is Whether the leading term has a double-pole depends upon s 2 = +1 , implying the propagator that appears twice is bosonic. If rather s 2 = −1 , then the pole is only simple. Equation A mass regulator λ helps to evaluate F 0 , by writing the repeated propagator as This calls for the one-loop results, provided in Appendix A, to be endowed with a mass on one of the propagators.
The integration methods are only slightly modified by the presence of a mass term above.
On-shell energies are then E 1 = |p| and E 2 = | λ | , and the integration limits are determined by energy conservation. We find, by generalising definition (A.1) to massive particles, The signed energies introduced are e 1 = sgn(p)E 1 and e 2 = sgn(q)E 2 . It is necessary to take the derivative of (B.2) with respect to λ 2 (and evaluate it at λ = 0). For this purpose the integration variables are changed to e 1 and e 2 so that λ-dependence is swept into the limits of integration. The δ-function then imposes the following limits on the p-integration: We only need to evaluate the integrand at the appropriate boundary values (with λ = 0), to obtain F . The result is stated in Eq. (3.4), in a way that is valid for either sign of K 2 .
Using the approach stated in Sec. 2.1 to include powers of p 0 into the above derivation, one obtains expressions for F 1 and F 2 . Casting them altogether gives (3.4). The resulting functions are plotted in Fig. 18, omitting the large-k 0 behaviour of (3.4), namely [see also Eq. (D.9)] We note that in this limit, the coefficient of the T 2 thermal contribution is zero. (Which is also the case for the statistical combinations not plotted.) Some relations can be derived for these functions. For the special case that s 1 = s 2 , one has F 1 = k 0 F 0 . In general, F 2 is actually expressible by F 1 and F 0 due to the identity: which can be used in  Fig. 17.

B.2. H m and its O( ) contribution
In this section, we explain the functions H To first establish H m from (3.5), it can be written as for m = 0, 1, 2 with the special notation of Appendix A.1. In the special case that s 1 = s 2 , one may show that H 1 = 1 2 k 0 H 0 by a change of integration variables. But in general, one must adhere to (A.4) for the integration. According to that formulation we can set v 1 = v 2 = +1 enabling the distribution functions to be written where the arguments were omitted.
The results are immediate: [The primed integral is defined in (A.5).]

B.3. Finite part of G m
We now consider G In the sum above, where v 1 = v 2 the ' 1 2 '-terms combine into what becomes the vacuum result. Anything proportional to the distribution functions gives thermal effects, and we use a symmetry to write n s 2 with the argument E 1 . Making use of (A.3), without imposing the energy constraint, it is necessary to also add the contribution coming from k 0 → −k 0 . We can complete all but one of the integrations to obtain For small p relative to the external momentum, Γ 0 −8p k/K 2 and therefore the integral is finite even if one of the distribution functions is bosonic. We have decided to work with all logarithms taking the absolute values of their arguments, but if we were careful to keep the correct sign one could use the imaginary part of (B.5) to determine H 0 .

Appendix C: Two-loop Kernels
In the diagrams labelled IV, V and VI, index c is nonzero and therefore the momenta running in the loop do not decouple. After carrying out the sums over p 0 and q 0 in (2.1), following e. g. Ref. [15], the terms can be sensibly collected according to which energies are on-shell. The main virtue of our strategy for m, n ≥ 0 (discussed in Sec. 2) is that it means we can perform the frequency sums assuming m = n = 0 .
For IV, the outcome is (4.2) with all intermediate particles on-shell. That is generally the form the 'real' corrections will take. Having more propagators (via d, e = 0) will facilitate other permutations of cuts. The type-V master integral has more, which are represented by the two terms in the following expression for the imaginary part: (The notation of the outermost integrals was defined in Appendix A.) Three momenta are put on-shell in the third line, now with an internal propagator that was not needed before.
The virtual correction [equal to the part of (C.1)] has factored into a binary decay amplitude multiplying a one-loop vertex amplitude. Note that the 'mass' λ 2 has been introduced as a regulator, so that R has on-shell energy E 3 = √ r 2 + λ 2 . Later we will take the limit λ → 0 and find that the real and virtual pieces dovetail together, leaving a result that is both finite and λ-independent. The other energies are as before: E 1 = |p| and E 2 = |q|, and we also The imaginary part of the special integral in (5.11) is given explicitly by The most intricate two-loop topology, yet benefiting from the most symmetry, as four different cuts contribute to the discontinuity, given by the following expression: The meaning of the third line should be clear, the second line is one of the virtual corrections.
It possess three terms which ought to be clarified. The momenta P and L = v 4 (K − v 1 P ) are on-shell with E 1 = |p| and E 4 = | | . Inside the curly braces we have defined Here we derive the expression given in (4.5), by explicitly carrying out the angular integrals from (4.2). The result of Appendix A, and specifically Eq. (A.14), shows that we can identify The meaning of h and the boundaries on the integrals were given there. One may safely set λ = 0 (it does not generate the collinear logarithms), which simplifies the kinematical constraints from h ≥ 0. Starting with the angular integration over x 1 = cos θ kp in (A.14), we write the function h = h(x 1 ) as a quadratic: ax 2 1 + bx 1 + c [33]. The coefficients, which can be calculated from (A.13), are Because a < 0 , the Θ-function in (C.5) dictates the upper and lower limits on the x 1 integration. We can parametrise x 1 in terms of the 'angle' ξ ∈ (0, π) with where ∆ = b 2 − 4ac is the discriminant. Changing the integration variable from x 1 to ξ thus removes the Θ-function and yields It turns out that ∆ ≥ 0 summarises the allowable phase space, which we have elaborated previously (in the momenta p and q). We already assumed a permissible (p, q) configuration by writing out the x 1 -integral.
Returning to (C.5), this gives The limits above follow from requiring ∆ ≥ 0 in (C.7), so that the integral has non-zero support. They are, for sanctioned p and q values, One may check that this reproduces the formula in (4.5).
The same approach works for calculating W V and W VI , however λ = 0 must be kept Let us note that the angular limits in (C.9) apply to all the real two-loop contributions.
Harking back to Figs. 6 and 7, these limits are specifiable in each region of the (p, q)plane. Considering them individually reveals how for k 0 k , the prima facie incompatible regions interchange when moving from above to below the light cone: If one takes the limit k 0 → k + 0 + (see Fig. 19), the x 2 -limits coincide thus giving zero precisely where p and q are kinematically forbidden in the spacelike region. A similar argument works out for k 0 → k − 0 + . Moreover, the angular limits are well-defined for k 0 = k in the nontrivial regions that remain (i.e. '2', '3' and '4' from Fig 5). That implies 9 the functions W (p, q) that we calculate are continuous at k 0 = k .
This does not preclude the final master integrals from being infinite on the light cone.
However, for those that are, the singularity must be the same whether approached from above or below. Some of the master integrals are finite at k 0 = k , others diverge either logarithmically or due to a pole.
The two terms in the second line of (C.1) may be calculated by a standard Feynman parametrisation. Let us manipulate these terms assuming v 1 = v 4 = +1 (it does not lose 9 We skipped over the technicality of setting λ = 0 . However the conclusion is still true: W is continuous for finite λ , as are the virtual corrections. Therefore once combined, and the limit λ → 0 is taken, they remain continuous at k 0 = k . any generality). We find where L = ( 0 , ) = K − P and R = L − v 2 Q . (Here, as before, we absorb the sum over v 2 into the sign of q.) Similarly, with a signed energy e 3 = sgn(r) √ λ 2 + r 2 . In this term, we change integration variables from r to q = k 0 − p − e 3 so that (C.11) and (C.10) are ready to be combined. But before doing that, a divergent contribution needs to be salvaged from this vertex correction. It is handily isolated (and underlined in what follows) by writing 1 2 + s 2 n s 2 (q) = sgn(q) 1 2 + sgn(q)s 2 n s 2 |q| .
in (C.10), and similarly a term like sgn(r) 1 2 in (C.11). The other parts lead to momentum integrals that are rendered finite by the thermal weights. But the divergent vacuum result, using dimensional regularisation, is equal to where L 2 = (K − P ) 2 = 0 was needed. This is what the curly braces in (C.1) produce in vacuum. With a large cut-off Λ on the magnitude of q , it is easy to show how the same type of divergence arises. The restricted integral, where we assumed λ → 0 . In this limit e 3 = r + λ 2 2r + . . . , which enabled the arguments of the logarithms to be simplified. Hence using log 4Λ 2 = −1 + logμ 2 is the choice consistent with the two-point function in (C.12).
Returning to Eqs. (C.10) and (C.11), these two contributions can be drawn up against their real counterparts, by substituting · · · = 1 (4π) 2 dq n s 2 n s 3 n s 4 log λ 2 4 q + 1 2 + s 3 n s 3 log q 2 r 2 . (C.13) The first summand combines with the real corrections to remove any dependence on λ overall. In the second, we again single out the divergent part and calculate it with a cut-off regulator. Namely, from which was converted to dimensional regularisation. Now O( ) terms in Eq. (A.4) must be kept for the outer integration. The sequestered part of (C.13) therefore generates −1 (4π) 2 1 + 2 logμ 2 K 2 + log (C.14) Equation (C.14) contains an ultraviolet divergence, which will end up being multiplied by a function of the temperature in the final expression. The first term in (C.1) can be written, after using signed momenta to incorporate the sum over the {v i }, The formula for U must be recovered by our regularisation procedure. This can be checked by redoing the calculation with a cut-off, as before. But now one finds that log 4Λ 2 = −1 + logμ 2 + 2 is needed to consistently convert the restricted integral, to dimensional regularisation. The argument is otherwise unchanged, and altogether leads to the formula for U V (p ; λ) that was given in (5.8).
The special master integral, defined in (5.11), also follows along the lines above. Here we give some more details: The curly braces of Eq. (C.2) give, for λ → 0 , Anything in the integrand that is proportional to a distribution function n s 2 (|q|) or n s 3 (|r|) (evaluated at positive arguments) will be finite. The divergent terms are easily isolated using the same approach as for U V . They can be calculated by a cut-off regulator, and altogether must give Once again, the hard cut-off integral can be given explicitly: This part of (C.16) should give the complete vacuum result for the vertex correction and consequently, should be converted to dimensional regularisation using log 4Λ 2 = −1 +logμ 2 + 4k 0 /K 2 . Going back to (C.16) and rearranging the distribution functions, the quantity in curly braces becomes, for λ → 0 , The first line above will combine with the real corrections; cf. Eq. (5.13). A remnant, which is proportional to sgn(r) 1 2 in the second line, diverges and can now be calculated with a regulator, With that, inserted into (C.17), we are able to single out a fragment that is the same as V up to a factor of K 2 / . We thus arrive at Eq. (5.14) of the main text, after using the principal valued integral to drop some polynomial parts of the thermal proportion. Moving along to the function U (n) VI (p ; λ) needed in (6.2), we derive it from (C.3). Let us focus on the three terms in curly braces and consider them individually. The first, for To carry out the v i -sums, the momenta p and q were extended to negative values. [Here only the q integration is made plain, the impending outer p-integration takes the form of (A.4).] Similarly we find, for λ → 0 , where the integration variable v has been extended to negative values. For the term where 'particle-3' is on-shell, we use the integration variable r = v 3 E 3 . Then, for λ → 0 , The restriction E 3 ≥ λ assists to explicitly regulate the r ≈ 0 singularity. Unlike in (C. 18) and (C. 19), where the corresponding integral may be understood in the principle valued sense, the risk that s 3 = +1 in (C.20) would make this futile. Therefore we keep λ to control the exclusion of this integration interval.
Let us change integration variables to q = k 0 − v instead of v and q = k 0 − p − r instead of r in (C. 19) and (C.20) respectively. The distribution functions can be rewritten in a way that the singular terms coalesce with their real counterparts.
We do this by using an identity to rewrite the integrand of (C.20), The two lines above are related by the exchange p ↔ and q ↔ v (together with same swap in statistics). Combining the result with (C. 19) and (C. 20) gives, for the curly braces in The case n = 2 needs to be handled with some care, as the subsequent q-integral diverges.
[Equation (D. 19) of Appendix D also exposes this fact.] That can be seen from the ultraviolet behaviour of the first line in (C.21), since for q → ±∞ the difference (n 2 − n 5 ) becomes ±1 , and the logarithm Hence the integration in (C.21) will be log-divergent for n = 2 , due to the 1/r = 1/(k 0 −p−q) left over. It can be attributed to the vertex correction, i. e. the three-point function studied in Ref. [34]. The entry of the rank-2 tensor that we need, is given by The ' 1 2 '-terms next to the distribution functions in Eqs. (C.18)-(C.20) concomitantly ought to reproduce this vacuum result. Explicitely taking these equations together, with an imposed cut-off for large-q and an infrared regulator for q , we find, for λ → 0 , Terms that vanish as Λ → ∞ were omitted. This integral appears in U VI (p ; λ) , to give (C.22) once the dust has settled and all factors are collected. This means that the cut-off regulator should be replaced by log 4Λ 2 → −1 + 2 logμ 2 + 2 .

Appendix D: Large-K 2 expansions
For external energies k 0 , that are much larger than both the temperature T and momentum k , spectral functions can be studied by OPE techniques [18]. The resulting approximations are applicable in the deeply virtual regime K 2 T 2 and may also be obtained systematically from the master sum-integrals themselves [35].
Carrying out the Matsubara sums of (2.1) produces terms, besides the vacuum result, with different loop momenta put 'on-shell' and weighted by a thermal distribution. These For general a, b, d, c and d, These five terms are compatible with the original symmetry of the diagram, e.g. A 2 with P ↔ Q gives a result with Q (associated with ∆ b Q in the original labelling) being the on-shell momenta. In A 1 , φ a (p) denotes the residue of ∆ a P at its positive pole p 0 = p , viewing the scalar propagator ∆(P ) = (p 2 0 − p 2 ) −1 as a function the complex energy. It may be expressed using the gamma function as φ a (p) ≡ (−1) a−1 √ π Γ(a + 1 2 ) Γ(a)p 2a−1 . (D.2) For the vacuum-like Q integrations it is safe to expand the integrand in P = (p 0 , p) .
Doing just that, for what we need 11 (dropping any terms we can, thanks to P 2 = 0) ∆ a K−P 1 K 2a 1 + 2a K · P K 2 + 2a(a + 1) (K · P ) 2 K 4 + . . . , After the necessary expansion is inserted into the definitions of each A i , one finds a variety of ordinary vacuum integrals. These types of integrals are all derivable from a class of 1-loop tensors [34], viz. (not the same m, n as before) In particular, Lorentz invariance implies that they are each linear combinations of independent tensor (of appropriate rank) that can be constructed from just K µ and g µν . We need only those with up to four indices, denoted J µ m,n = A m,n K µ , J µν m,n = C m,n g µν + B m,n K µ K ν , J µνρ m,n = E m,n K m g νρ + sym. + D m,n K µ K ν K ρ , J µνρσ m,n = H m,n g µν g ρσ + sym. + G m,n K µ K ν g ρσ + sym. + F m,n K µ K ν K ρ K σ .
The coefficients A, B, C, D, E, F, G and H can all be related to the fundamental scalar integral J m,n (it has no powers of q 0 in the numerator).
Without loss of generality, we assume m ≥ n so that the pairs (m, n) of interest here are 11 We write four products X · Y = x 0 y 0 − x · y .
P µ J µ m,n = (K · P )A m,n J 0 m,n = k 0 A m,n (D.4) P µ P ν J µν m,n = (K · P ) 2 B m,n P µ J µ0 m,n = p 0 C m,n + k 0 (K · P )B m,n J 00 m,n = C m,n + k 2 0 B m,n (D.5) P µ P ν J µν0 m,n = 2p 0 (K · P )E m,n + k 0 (K · P ) 2 D m,n P µ J µ00 m,n = (K · P ) + 2k 0 p 0 E m,n + k 0 (K · P ) 2 D m,n (D.6) P µ P ν J µν00 m,n = 2p 2 0 H m,n + (K · P ) (K · P ) + 4k 0 p 0 G m,n + k 2 0 (K · P ) 2 F m,n (D.7) These expressions still depend on the relative angle between k and p . That is to be taken into account when performing the integrals in (D.1), which we carry out using e. g. into the form Im I ∼ ω 0 K 2 + ω 2 T 2 + ω 4 T 4 K 2 + . . . , (K 2 /T 2 → ∞) where ω i will depend on details of the master integral. The first coefficient, ω 0 , is exactly the vacuum result and therefore might contain a term ( −1 + 2 logμ 2 /K 2 ) . Coefficients of powers of T 2 , which are all finite, arise from moments of equilibrium distribution functions.
So, for instance, ω 4 stems from p pn s i (p) = O(T 4 ) . Only ω 0 is independent of the statistical nature (via s i = ±1) of the propagators.
With the abbreviation n i ≡ s i n s i (p) , we list some of those integrals now: The masters that factor into a product with a tadpole diagram are zero in vacuum: they only start at O(T 2 ) . It turns out that they also have no T 4 -term, Im I One may also obtain the masters I 11010 and I 11020 by replacing s 3 with s 2 in the above.
The sunset integral has the expansion Im I which is evidently symmetric in {s i } for the indices i = 1, 2 and 3 . It also has no T 4 -term, and the T 2 -term can equal zero if only one of the particles is bosonic, the rest fermionic.
For the spectacle diagram, which bears an ultraviolet divergence (from the 1-loop factor) that persists after taking the imaginary part, we find 3 1 + 2 logμ 2 K 2 + 4 (D.11) + p n 1 + n 2 + n 4 + n 5 16πp + p p(n 1 + n 2 + n 4 + n 5 ) 4π K 4 k 2 0 + The result is symmetric in {s i } for i = 1, 2, 3 and 4 . Moreover, given any combination of statistics, the T 2 -and T 4 -order corrections are not zero. The leading term is the vacuum result and thermal corrections would start at O(T 2 ), but this term is absent in accordance with [18].