The Sun Also Rises: the High-Temperature Expansion of the Thermal Sunset

We give a prescription for calculating the high-temperature expansion of the thermal sunset integral to arbitrary order. Results up to $\mathcal{O}(T^0)$ are given for both bosonic and fermionic thermal sunsets in dimensional regularisation, and for all odd powers of $T$ up to order $\epsilon^0$. The methods used generalize to non-zero external momentum. We verify the results with sundry analytical and numerical cross-checks. Intermediate steps involve integrals over three Bessel functions.


Introduction
Finite-temperature field theory has a wide range of applications, from phase transitions to the inner workings of neutron stars. Although the field was established last century, the techniques evolve constantly to deal with an ever-increasing demand for precision. And fresh applications. Finite temperature calculations are perturbative whenever possible, though lattice calculations are also viable.
Yet perturbative calculations are arduous and often involve many mass scales. Fortunately, many applications feature a hierarchy between particle masses and the temperature: T 2 X ; a high-temperature expansion is applicable. This idea works wonders at one loop, but general 2-loop high-temperature expansions are only known to leading order in T 2 -higher orders are intricate.
To circumvent this problem one can resort to numerical evaluations. Which is a viable strategy for some cases. Yet this leaves something to be desired when the temperature follows a strict power counting. Not to mention the computational complexity of numerical integrals; especially if numerous mass scales are present.
So in this paper we present a new method to calculate a high-temperature expansion of the bosonic and fermionic sunset diagrams to arbitrary order in T . As intermediate steps, sums and integrals of three Bessel functions are given. We apply these methods to sunsets without external momentum, and derive explicit formulas up to O(T 0 ). To ensure the validity of the results we perform sundry theoretical and numerical cross-checks.

Bosonic sunset
The thermal bosonic sunset for three arbitrary squared masses X , Y, Z is I(X , Y, Z) ≡ T 2 n p ,n q ,n l p, q, l δ( p + q + l)δ n p +n q +n l ,0 with the measure p = we use dimensional regularisation in the M S scheme, and take µ as the M S scale. More compactly,

2)
Expand the sunset as where the divergent contributions are known to all orders in T , and the finite piece is previously known to O(T 2 ) [1,2], (2.5)

7)
A(X ) is the finite piece of the bosonic 1-loop bubble integral, see equation (3.11).

Fermionic sunset
The fermionic sunset is n p ,n q ,n l p, q, l δ( p + q + l)δ n p +n q +n l ,0 × 1 p 2 + X + (π(2n p + 1)T ) 2 Here propagators involving P & X and Q & Y are due to fermions-evidenced by the odd Matsubara frequencies.
Again, expand in ε, Divergent pieces are known to all orders, and the finite piece is zero to O(T 2 ), where A(Z) and A F (X ) are the finite parts of the bosonic and fermionic 1-loop bubble integrals; see equations (3.11) and (3.13).
The fermionic result of this paper is (2.14) Note the different logarithms for a bosonic versus a fermionic mass. Sections 3-5 give detailed derivations of the above results.

Gaining Momentum
Temperature dependence arises as multiplicative factors from Feynman rules and from propagators. The latter is an involved sum over Matsubara frequencies. Our approach is to perform the hightemperature expansions before the Matsubara sums-not after. In this section we discuss this approach and introduce handy labels for the derivation to come. We give examples of how the high-temperature expansion works at one and two loops.

Hard/Soft split
Thermal integrals are often evaluated by first getting rid of all Matsubara sums. One advantage of this approach is that vacuum and thermal contributions are clearly separated. Also, because there are no sums left, remaining integrals can be performed numerically in the absence of analytical results.
To derive the high-temperature expansion of the sunset, we'll instead do the opposite: expand in T before doing the Matsubara sums-the high-temperature expansion gets done straight off the bat. The problem is reduced to sums over master integrals. Separate momenta into hard and soft [3] Hard p 2 ∼ T 2 , Where a hierarchy between T and the masses is assumed: T 2 X , Y, Z. Consider the bosonic propagator. There are two cases for hard momenta. First, for a finite Matsubara mode, Second, for a Matsubara zero-mode, Likewise for soft momenta with a finite Matsubara mode, and with a Matsubara zero-mode Similar relations hold for the fermionic propagator.

1-loop momentum split
Take the bosonic 1-loop bubble integral The traditional method proceeds by summing over Matsubara modes [4], This is swell since vacuum and temperature parts are separated; the Bose factor isolates the thermal part. The 4D part is readily evaluated where µ is the M S scale. There are various ways to evaluate the temperature integral; one is to z z 2 +(2πl) 2 to expand the Bose factor. After a fair bit of work one finds Adding the vacuum and thermal parts together gives the full result (3.10) Note that all ε-poles come from the vacuum part; thermal contributions can not diverge in the UV at one loop. There is a corresponding result for the fermionic 1-loop bubble, This method is fine at one loop, but unnecessarily complicated. First, Matsubara sums are done and then effectively reintroduced in the high-temperature expansion. Second, the expansion of the Bose factor is messy. This is where the hard/soft split comes in. The high-temperature expansion and the momentum integration are done before the sums. To wit where H and S stand for hard and soft momenta respectively. Start with the hard contribution, Only standard integrals and sums were necessary. The soft contribution is where all contributions from n = 0 modes vanish due to scaleless integrals. Only zero-mode terms survive when all momenta are soft, this stays true at higher loop orders. One can similarly apply the hard/soft-split to the fermionic bubble. Though in this case there is no zero-mode and hence no soft contribution.

Sunset Momentum Split
There are two momentum integrals for the sunset, and so a few more cases. First take all Matsubara modes to zero. We use the notation F ≡ P,Q (3.17) where F H HS is defined so that p, q are hard, and l is soft. There are no terms with two soft momenta due to momentum conservation. Furthermore, F H H H = F H HS = 0 to all orders. For example, take l soft, since all terms multiply a scaleless integral. Only the all-soft contribution is finite. Next consider one zero and two finite modes, say n p = n q = 0, n l = 0. Denote these as G l ≡ P,Q,L Again split the integral into different momentum regions, In this case only G l SSS vanishes to all orders.
The leading order (in T ) comes from G l H HS and is of order T . Explicitly, other hard/soft momenta assignments of G l H HS are zero to all orders of T ; G l H HS is given to all orders in the next section.
Finally there's the case when all Matsubara modes are finite. We use the notation D ≡ P,Q Note that D H HS = 0,1 and D SSS = 0 to all orders. It is well known that In the momentum split this statement means-to leading order in T -that F H H H +3G l H H H +D H H H = 0, which is confirmed in the next section.

The Bosonic Sunset
There are two families of terms corresponding to odd and even powers of T .

Odd powers of T
The starting point is where p ∼ q ∼ T , and l 2 ∼ Z T 2 . There are two sources of higher-order terms. The first type comes from using the delta function and expanding (p + l) 2 in powers of l 2 (odd l terms integrate 1This is a bit subtle, and care must be taken when evaluating the overall momentum delta function. For example, if we collapse the delta function as q = − p − l , then all integrals are scaleless since l is soft. If the delta function is collapsed the other way, l = − p − q, it seems like the result is finite. This is however only an illusion since for l to be soft we must have p = − q + k, where k is soft. And all integrals vanish with this change of variables. to zero). These terms give corrections of order Z T 2 . Remaining terms come from expanding the finite-mode propagators in orders of X and Y . Ignoring O (ε) corrections, the general result is where α = 0 is a special case, due to an ε-pole, and must be treated separately, see equation (3.21). Note that G l H HS is symmetric in X and Y as it must be. For example, the T −1 and T −3 contributions are respectively

Even powers of T
The derivations of the even T terms are way more involved, and arguably more interesting.
The momentum integrals are more lucid after Fourier transforming to coordinate space [5], The coordinate space representation is then Performing this replacement in the sunset eliminates all momentum integrals. Leaving a single R integral. The propagator's coordinate representation in any dimension is derived in appendix C.

Order T 2
There are a priori three different contributions at order T 2 . The all zero-mode contribution F SSS , and the two finite mode contributions G l H H H , D H H H . Finite mode contributions do not depend on the masses and are zero. And so F SSS gives the full order T 2 result. The integral can be performed using standard techniques [6]. Now for showing that the mass-independent T 2 term vanishes. Start with G l H H H , Ignore T 0 terms for now, G l H H H simply refers to T 2 terms in this section. From now on, unless otherwise specified, Matsubara modes n, m, and l always refer to the absolute value of said integer, and are all positive. Transforming to coordinate space we find The second integral is Rescaling the momenta and going to coordinate space, The integral is known in closed form [7,8], Yet there is a subtlety. The formula is technically not valid when l = n + m. Which is precisely the case of interest. All is not lost, though-the formula holds in dimensional regularisation. The trick is to enforce the Kronecker delta by setting l = n + m − δ 2 (l is here the absolute value), and then take the δ → 0 limit. All δ dependence cancels to O δ 2 , so the δ → 0 limit is valid.

So adding the two contributions D H H H + G l H H H gives
The numerical factors above come from different summation regions in the first case, and permutations of the zero-mode in the second. There are actually 6 summation regions in total: (4.16) Using the methods in appendix A the second sum is ∞ n,m=1 Add everything together to find So the mass independent T 2 term vanishes as promised.

Order T 0
Order For clarity G l H H H refers to the T 0 term for the remainder of this section. The trick to evaluating these integrals in coordinate space is to rewrite powers of the propagators as derivatives acting on it, as in (4.20) Use this trick to suss out There are also contributions from permutations of the zero-mode. Adding them together and expanding in ε gives Note that there are no ε-divergences. This is not surprising. Before the sum, by dimensional reasons, the n power must be n −4ε−2 .2 This can never diverge; similar for higher-order terms. The situation is different however when three modes are finite. On dimensional grounds the sum must be of the form (4.24) 2Because Matsubara modes only appear in the combination nT .

Mimicking G l H H H , D H H H only denotes the T 0 term. In coordinate space the integral is
The integral can again be done in closed form. Rewrite the squared propagator using the trick above, and then use the relation where the second integral is given by equation (4.10). Again we imagine collapsing the Kronecker delta by setting l = n + m − δ 2 , and taking the δ → 0 limit. There are some new subtleties because of δ −2 terms, but these all cancel out, and the δ → 0 limit can be consistently taken. After some simplifications one finds n [mn(m + n)] 2ε+1 2n(ε + 1)m 2ε+1 + (2ε + 1)m 2ε+2 − m 2 (2ε + 1)(m + n) 2ε +n 2 (m + n) 2ε − n 2ε − 2mnε(m + n) 2ε , (4.28) B = − 2 −2(ε+4) e 2γε sec(πε)Γ (2ε + 1) This result can be double-checked order-by-order in ε by expanding the Bessel functions before integrating. Using this we have double-checked all integrals to O (ε). Again, same as for the T 2 term, there are 6 summation regions. Note that region (iii) gives the same contribution as (i) since they are related by relabelling dummy indices. Region (ii) is different-because the negative mode is in the P −4 propagator. We have to do this case separately. So yet another integral is needed: With the same B as in equation (4.28). After doing the sums we get the neat result

The result of adding D H H H and G l H H H + permutations is given in (5.26); or after expanding to
So all ε poles come from D H H H . While the full O T 0 term is new, all ε poles are previously known [1]. These poles agree with our results to O T 0 -giving a swell cross-check.

The Fermionic Sunset
The fermionic sunset's high-temperature expansion is similar to its cousin's. But now only the l propagator can have a zero Matsubara mode; and so only the l momentum can be soft. The book-keeping is also more complicated, because we have to keep track which Matsubara modes can be odd or even.
To make life simpler we use the same notation G, D. In this case

Odd T contributions
Odd T powers are calculated analogously to the bosonic sunset. All arise from the integral The O(T ) contribution is Ignoring O (ε) corrections, the general result for higher orders is

T 2 Contribution
As fortune would have it we don't have to do any work for the T 2 contribution; just write down the answer. First, two Matsubara modes are always finite, so the T 2 contribution can't depend on masses. Second, the remaining contribution is of the form And we can use the summation trick in [9] to reshuffle things as But the right-hand side is zero-so the fermionic sunset is zero at O T 2 . Nevertheless-since it introduces methods we will need later-let's explicitly show that the T 2 contribution vanishes to all orders in ε.
There are two contributions at order T 2 . First, In this case two modes are fermionic, and one is bosonic. Without loss of generality we'll take n as the bosonic mode. There are six summation regions: Again, the convention is that the Kronecker delta collapses upon the negative mode. So in region (i) l is odd and is l = m + n; m is odd and n is even. All flipped regions give an overall factor of two and region (i) and (iii) are identical. Up to an irrelevant prefactor, the result is Where eo 1 n a (n+m) b stands for the sum over (positive) even n and odd m. Adding G l H H H gives (5.14) This vanishes after using identities from appendix A. So the fermionic sunset is indeed zero at order T 2 .

T 0 Contribution
All contributions at order T 0 come from G l H H H + D H H H . First, note that if all squared masses are identical (= X ), where all contributions from region (i)-(iii) have been added. All sums can be done in closed form using various identities in appendix A. In terms of ζ-functions the result is neatly expressed as So after expanding to O ε 0 the bosonic mass contribution is Let's continue with the contribution of a fermionic mass, D X H H H . Following the same steps as before we suss out with B F as in equation (5.18). Including both fermionic masses, the proper prefactors, and expanding to O ε 0 , And for the fermionic mass terms,

Comparison of Bosonic and Fermionic
Let's now compare the complete bosonic and fermionic sunsets at O(T 0 ). The bosonic sunset is And the fermionic sunset is (5.28) Comparing the bosonic and fermionic sunset we confirm that indeed I F (X , X , X ) = (4 1+2ε − 1)I(X , 0, 0) at O(T 0 ): giving a highly non-trivial cross-check. Finally, it is quite remarkable that all ζ(4ε + 2) terms cancel between D H H H and G H H H , in both the bosonic and the fermionic case.

Numerical Tests
With analytical calculations out of the way it is interesting to compare them to numerical results. Luckily the full result is known in terms of definite integrals [1]. The numerical integration is fast when the three masses are of similar order; not so much when there is a hierarchy between masses.
Start with the bosonic sunset. To ease comparisons we set all masses the same, and then plot I(X , X , X )/T 2 versus X /T in figure 1. Each subplot corresponds to a different scaling µ = c T , as noted in the figures.
In figure 2 we instead focus on the fermionic sunset. In this case we chose the boson mass to be half of the fermion mass, to mark that these are different particles.

Conclusion
In this paper we have introduced techniques to perform 2-loop high-temperature expansions in dimensional regularisation, applicable for zero and non-zero external momentum. We explicitly derived the O(T ) and O(T 0 ) contributions to the bosonic and fermionic thermal sunset integrals, and all contributions of odd powers of T to order ε 0 . As far as we know, this result is not previously available in the literature-though see [10] for some partial results. All terms, at a given order, are expressed in standard functions. These methods are likely useful even at higher loops. Testing them on 3-loop thermal integrals is an avenue of future research. Various analytical cross-checks and numerical comparisons validate the results. There is a stark agreement for a range of masses. The previously known approximation to O T 2 is, as noted in the past, quite poor for any sizeable mass. The next O (T ) correction is rather important and extends The bosonic sunset is described well at O T 0 , for masses up to ∼ 6T . The range of validity is shorter for the fermionic sunset, where the accuracy is reasonable for masses up to ∼ 3T . Though note that these conclusions depend slightly on the renormalization scale. And a clever choice for µ extends the range of validity somewhat.
There are a number of uses for the result in this paper. In high-temperature calculations where the size of T plays a role in the power counting, such as in EFT-and resummation-techniques, it's not convenient to use the numerical evaluation of the integrals. The method presented in this paper can extend the reach of perturbation theory.
These results can also be used to improve numerical calculations of the sunset integral itself. As noted in [1], care must be taken when there is a hierarchy between the masses. A numerical calculation can use the expansion we have derived in such a region of parameter space.
Finally, it is quaint that all Matsubara sums can be done in closed form in terms of ζ functions-at least up till O T 0 . This suggests that something deeper is going on, and could be interesting to explore in the future.

Acknowledgments
We would like to thank Renato Fonseca for many interesting discussions.

A Sums & ζ-functions
Many sums used in this paper are of the form ∞ n,m=1 Sums of this type are straightforward to evaluate [11]. For example, start with the more general sum Putting everything together, ∞ n,m=1 Choosing a = b gives the above sum. We also need another class of sums for fermionic sunsets. These are of the form And similar for other combinations of even and odd. These sums are straightforwardly evaluated by using [11,12]  This integral is intractable in general; yet it turns out that the leading terms, for all sums considered in this paper, come from the x → 0 region; for which the integral is readily evaluated. Yet higher order ε-corrections might be useful in the future. So let's outline how these corrections can be obtained. An example speaks more than a thousand words, so consider the sum ∞ n,m=1 n −1−ε m −1−ε (n + m) −2ε . Using the Feynman trick, Now, there're are two ways to evaluate this integral: the fast way, and the systematic way. The fast way only gives the leading terms in ε; the systematic way is more cumbersome but enables one to calculate the sum to an arbitrary order in ε.
The fast way uses that ε-poles arise from the x → 0 region; so let's introduce a cut-off R and isolate the poles, The systematic way first subtracts the small-x divergences: The integral is performed with a regulator; we choose the same as in [13]: where the exponential e −2x is chosen to reproduce the asymptotic behaviour of Li 2 ε+1 (e −x ). The divergences are then Using this method one can evaluate all possible sums arising in the high-temperature expansion.

C Coordinate space propagator
The propagator's Fourier transform is (C.1) We'll proceed in two steps. First, the angular integration, and then the "radial" integration. So jumping directly into the fray, recall the volume element of a d-dimensional sphere To make life easier, orient the coordinate system so that p · R = pR cos φ 1 ; for the angular integration then splits into one integration over φ 1 times the solid angle for a d − 1 sphere. Explicitly, where J is the Bessel function of the first kind. All that remains is with K being the modified Bessel function of the second kind, and the µ term added as usual in dimensional regularization.