Factorizing two-loop vacuum sum-integrals

We derive analytic results for scalar massless bosonic vacuum sum-integrals at two loops. Building upon a recent factorization proof of massive two-loop vacuum integrals, we are able to solve the corresponding Matsubara sums and map the result onto one-loop structures, thereby proving factorization also in the sum-integral setting. Analytic results are provided for generic integer-valued propagator- and numerator-powers of the class of sum-integrals under consideration, allowing to eliminate them from any perturbative expansion, dramatically simplifying the evaluation of some observables encountered e.g. in hot QCD.


Introduction
Due to their phenomenological importance in early-universe cosmology, compact-star astrophysics, and heavy-ion collisions, substantial effort has been spent on evaluations of quantities relevant in hot and/or dense systems in and near equilibrium.The main emphasis has been on understanding quantum chromodynamics (QCD) [1] or generalizations thereof in such extreme conditions, since the interplay of asymptotic freedom at high energies (temperatures T and/or densities as defined by the fermion chemical potentials µ f ) and strong coupling at low energies provides for a very rich scenario, in which an interplay of different theoretical methods becomes important.
Within finite-temperature quantum field theory (see [2][3][4][5] for textbook-level treatments), the evaluation of vacuum-type integrals plays a central role in the determination of equilibrium observables, such as the equation of state (EoS) of a thermal system [6,7].The EoS is typically given in terms of the pressure (or free energy) of the system under consideration, and constitutes a crucial ingredient in data analysis and modeling in the phenomenological applications mentioned above.First-principles determinations of the QCD EoS via Monte Carlo simulations of lattice gauge theory [8][9][10] are difficult in particular when fermions are present, leaving weak-coupling expansions as a viable alternative.
One very fruitful weak-coupling approach, for example, has been the application of effective field theory (EFT) methods to hot and dense systems, in particular in the nonabelian setting of QCD, where the gauge symmetry severely constrains (but also guides) methods that go beyond naive perturbative expansions 1 .This has led to the development and application of dimensionally reduced EFT's [17][18][19][20] that describe the strongly coupled low-energy QCD degrees of freedom in a systematic way, taking account of the correct high-energy behavior by mapping all EFT couplings onto QCD parameters (strong gauge coupling, temperature, group structure and field content, fermion masses and chemical potentials).Sometimes, and in particular for bosonic systems, comparison with lattice simulations are possible and provide valuable insight into questions such as convergence and range of validity of perturbative series.
For the QCD pressure, for example, almost all contributions up to the order of g 6 (in a formal power-counting in terms of the gauge coupling g, corresponding to the four-loop level in a diagrammatic expansion) that arise from QCD as well as two levels of EFT's that correspond to screened and unscreened 'soft' gluonic modes are known [21][22][23], except for a well-defined but technically challenging perturbative piece from the 'hard' QCD sector.While 3d lattice simulations already play a role in determining parts of the low-energy contributions [24], it has even been possible to obtain an order-of-magnitude estimate for the missing perturbative piece from comparisons with 4d lattice simulations at intermediate temperatures where both, weak-coupling and lattice approaches are justified [25,26].From a theoretical standpoint this is not a satisfactory situation, however, strongly motivating the development of techniques that allow for an evaluation of the missing perturbative contributions, see [27] for the state-of-the-art at µ f = 0, and [23] for µ f ̸ = 0.
The aim of the present paper is to advance perturbative techniques that directly apply to such contributions, necessitating the evaluation of many sum-integrals (for a definition, see section 2 below).As has been repeatedly observed in the past, many two-loop sumintegrals can be decomposed into one-loop factors that are known in terms of Zeta and Gamma functions, allowing for analytic solutions in the space-time dimension D = d + 1. See section 6.2 for a number of known examples from the literature.In the present paper, we prove that this decomposition is generic, and give an algorithm that constructs this decomposition for any scalar massless bosonic two-loop vacuum sum-integral.
In the fermionic case, the situation is less clear.While again some isolated two-loop sum-integral factorization formulas have been observed in the literature, when both scales T and µ f are present a recursive treatment might still apply to reduce the generic two-loop case to one-loop factors [28].
The structure of the paper is as follows.After introducing some basic notation in section 2, we introduce a generic massless two-loop vacuum sum-integral as the main object of our paper in section 3, and immediately represent it as infinite double (Matsubara-) sums over of massive two-loop (continuum-) integrals.For the latter, for generic positive values of propagator powers and in the specific mass-combinations of interest here, a factorization formula is known, whose details we recall in section 4. In section 5, we then show how to manipulate the double sums such as to generalize the factorization formula to the sumintegral setting, via four separate lines of proof.We finally present our factorized result as well as an extension to non-positive propagator powers in section 6, and close with a brief discussion of possible generalizations in section 7. Two short appendices collect details about single and double sums of interest, as well as sample Mathematica code that implements our reduction algorithm.

Notation and one-loop example
Let us start with some definitions, to set the stage for the rather technical nature of the present paper.Our setting is thermal field theory in the imaginary-time formalism, defined on a (d + 1) dimensional Euclidean space with a compact temporal coordinate τ (which emerges after Wick rotating it → τ , leading to e iS → e −S E in the path integral) that has period 1/T set by the inverse temperature.Dimensionally regularized integrals then employ d = 3 − 2ε in our conventions.Since we do not ε-expand but work fully analytically in this paper, the particular integer value of the space-time dimension will not play any role below.Bosonic (fermionic) fields are periodic (antiperiodic) functions of τ .This leads to discrete temporal Fourier transforms to momentum space, such that momentum integrals are traded for sum-integrals.We work with Euclidean four-momenta P = (P 0 , p) with bosonic2 Matsubara frequencies P 0 = 2n p πT where n p ∈ Z. Massless propagators are 1/[P 2 ] = 1/[(P 0 ) 2 + p 2 ], and we employ the sum-integral symbol defined as While some of the machinery used in organizing and reducing perturbative expansions, such as e.g.integration-by-parts (IBP) methods [29,30], carries over from zero to finite temperatures [31], evaluating master integrals is severely complicated by the Matsubara sums encountered at finite T .In particular, the number of analytically known sum-integrals (mostly simple one-and two-loop cases) is rather limited, while a number of phenomenologically relevant cases at three and four loops have been worked out up to a couple of terms in their ε expansion, typically relying on tedious subtractions of sub-divergences in one-loop substructures [6], and resorting to numerical evaluations for the constant terms.
As one concrete well-known example for an analytic solution, let us warm up by discussing the somewhat trivial 1-loop massless bosonic vacuum sum-integrals (η ∈ N 0 ) One way to evaluate them is to view the (d + 1)-dimensional sum-integral with massless propagators as a sum over massive d-dimensional 1-loop tadpole integrals (m ≥ 0) Taking advantage of the simple mass dependence of the continuum integral J, this gives where we have first separated the point n = 0 where propagator mass vanishes (leaving a massless tadpole that vanishes identically in dimensional regularization), mapped negative Matsubara indices onto positive ones, and used eq.( 2.3) to scale out the mass to finally identify the single Matsubara sum as a zeta function.As it should be, due to the sumintegral's symmetry under a momentum shift P → −P , all I η odd ν (d, T ) vanish identically, as is reflected by the prefactor [1 + (−1) η ].Therefore, we only encounter zeta values of the form ζ(n even − d) at one loop3 , a fact that will become important to recall in section 6.
The goal of the present paper is to expand the set of known generic-index analytic solutions to the two-loop level, going beyond the finite number of fixed-index cases that had been solved by IBP methods in the past.In analogy to the 1-loop case above, the main idea is to view a massless multi-loop sum-integral as a multiple sum over massive d-dimensional momentum integrals, where Matsubara frequencies such as P 0 play the role of propagator masses.

Massless bosonic two-loop vacuum sum-integrals: Setup
The general scalar massless bosonic 2-loop vacuum sum-integral is defined as Imitating the strategy presented in the evaluation of eq.(2.4) above, one can view the massless two-loop sum-integral L as a double sum over massive continuum integrals B, with the Matsubara frequencies playing the role of the continuum integral's masses, as where In general, one would expect the massive integrals B to have a complicated dependence on (ratios of) its three masses.In fact, general results in terms of Appell's hypergeometric function F 4 can be found in [32].In our particular case, however, these masses are seen to be linearly related, due to 4-momentum conservation at the sum-integral's vertices.For linearly related masses there are significant simplifications, as recently pointed out in [33] by a detailed analysis of the corresponding IBP relations, leading to a general factorization formula of the type B → J • J.The idea is now to use results of [33] for the massive 3d vacuum tadpole B of eq.(3.2), and do the double-sum afterwards.Amazingly, this will allow us to present a proof of factorization for, schematically, L → I • I, where the I are the 1-loop sum-integrals of eq.(2.4).
Further following the successful 1-loop strategy as explained below eq.(2.4), let us start by sorting out the cases where masses of B can vanish.To this end, we decompose the double-sum into sectors where the masses are always positive (see figure 1 for region mapping), taking into account that B depends on squared masses only, and use the integral's symmetry to re-order some mass-index pairs in B. We obtain4 and Hνa,ν b ,νc ηa,η b ,ηc (d) ≡ Naturally, the sum-integral L vanishes for odd Ση i , as is immediately clear from the definition eq.(3.1), considering the integrand's symmetry under the simultaneous shifts in d dimensions [34], while can be extracted from [32] (cf.appendix C of [33]).Note however that these specific results are not necessarily relevant for our treatment of L, as can be appreciated from our final result eq.(6.2) below.

Massive two-loop vacuum integrals: Factorization formula
Let us now examine the massive two-loop scalar vacuum integrals B ν 1 ,ν 2 ,ν 3 m 1 ,m 2 ,m 3 (d) that were defined in eq.(3.3), over whose masses we need to sum.As already mentioned above, general results in terms of Appell's hypergeometric function F 4 of two variables (mass ratios) can be found in eq. ( 4.3) of [32].This is not very practical for our purpose, since it involves four infinite sums.One can actually do (much) better by analyzing the integrals B from scratch, at the special kinematic point that is of relevance here.Stemming from the Matsubara sums, in the case of interest here the masses obey a 'sum-rule' at each vertex; in fact, according to the mapping in eq.(3.4) we only need m 3 = m 1 + m 2 .For this special "collinear" case, the general results of [32] can be simplified considerably.
Imposing the linear relation between the three masses, one can set up an integrationby-parts (IBP) recursion that reduces the integral B for all positive values of ν i ∈ N towards one or zero.Relevant explicit formulas can be found e.g. in eqs.(92,93,97) of [35].We have recently been able to solve these IBP recursion relations explicitly [33], and obtained (from now on we implicitly assume m 3 = m 1 + m 2 in all equations; also, Σν where the coefficients c νa,ν b ;j (d) are rational functions in d that obey a number of symmetries, such as c and are given explicitly by [33] c with Pochhammer symbols (a) ν ≡ Γ(a+ν) Γ(a) and integers n j = ⌈ Σν i +j 2 ⌉.Note that using eq.(4.1) in eq.(3.4) explicitly decomposes the two-loop sum-integral L into a sum over one-loop factors, a fact that we will exploit in section 5 below.
For later use, we spell out eq. (4.1) for two special cases that appear in eq.(3.4): with with and 5 Massless bosonic two-loop vacuum sum-integrals: Matsubara sums The final task is now to evaluate the remaining double sums H and H in eq.(3.4).Having eq. ( 4.1) at hand, the mass structure is fixed explicitly, and we can work on performing the sums without actually specifying the rational coefficient functions c νa,ν b ;j (d) yet.
We will now show how the specific linear combination of sums H and H that occurs in eq.(3.4)  all of which allows us to formulate our final result in a compact way.

Proof of statement (a)
Plugging eq. ( 4.1) into eq.(3.5), re-expressing the 'unwanted mass' (the one not already appearing in the j-sum) in the prefactor in terms of the two others (by introducing an additional binomial sum, as for example in , shifting summation indices appropriately (e.g.k → η 3 − k and j → −j) and using the symmetry eq.(4.2) to collect similar expressions, the last two terms of eq.(3.4) combine to where we have abbreviated ℓ p ≡ Σν i − η p − j − k, and the ηp k are binomial coefficients.The double sums in the last two lines of eq.(5.1) have the form To get to the third line, we have shifted n 1 → n 1 − n 2 and n 1 ↔ n 2 in the second and third terms of the second line, respectively, whereupon a sum that we cannot solve cancels and leaves us with (multiple) zeta values, see eq. (A.2).The double sum in the first line of eq. ( 5.1) has the form where the solution follows from the shuffle identity eq.(A.6).

Proof of statement (b)
In a similar way, plugging eq. ( 4.1) into eq.(3.6), re-expressing the 'unwanted mass' in the prefactor in terms of the two others, shifting summation indices appropriately and using the coefficients' symmetry, , (5.4) (5.5) The double sums needed here are given in eqs.(A.2), (A.3) and (A.4), respectively, resulting in single and double zeta values.Combining with the prefactors corresponding to eq. (3.4), we therefore have Hν Here, Note that the prefactor of eq.(3.4) guarantees that Ση i is even, so we can drop the (−1) Ση i in the last line above.Now, all multiple zeta values occur in symmetric pairs, and can hence be converted to single zeta values (and products thereof) by the shuffle relation eq.(A.6):

Proof of statement (c)
Looking at the structure of eq.(5.7), we are happy with the products of zetas (which we can in the end interpret as products of 1-loop sum-integrals eq.(2.4)), and note that all remaining single zeta values have arguments e p + d p = 2Σν i − 2d − Ση i that do not depend on the summation indices j, k.For these terms, the k-sum can be easily performed using producing Comparing with eqs.(4.4)-(4.7), the remaining sums over j can now be seen to be nothing but the coefficients of special cases of B, , the basis tensor contractions; leading to The 1-loop tadpoles I in eq.(6.5) can finally be mapped onto the basis Î via eqs.(2.4) and (6.1).On the one hand, this introduces Gamma factors in the denominators which correctly set all I η ν≤0 = 0. To exclude such vanishing terms from the outset, we can also modify the limits of the k-sum as min(−ν 3 −n,ν 2 −1) k=max(0,1−ν 1 −ν 3 −n) .On the other hand, the map introduces prefactors which implement that I η=odd ν = 0. Modifying the j-sum as where x = mod(η 2 + n, 2) then eliminates the need for such explicit prefactors [1 + (−1) η ] stemming from eq. (2.4).After some slight simplifications, this leaves us with a representation on par with eq. ( 6.2): where x ≡ mod(η 2 + n, 2) as well as y ≡ k − j − ⌊(η 2 + n)/2⌋ and z ≡ Σν i + n − k.As a check, we note that for ν 3 = 0 and positive values of ν 1 , ν 2 , eq. (6.6) as well as eq.( 6.2) with eq. ( 6.3) both reduce to single sums 5 , and agree in their factorization of L η 1 ,η 2 ,η 3 ν 1 ,ν 2 ,0 (d, T ).We can therefore use the much more general eq.(6.6) for reducing sum-integrals involving non-positive indices, without ever applying eq.(6.3).In order to check eq.(6.6), we have verified it against a pre-existing in-house IBP code, for a large number of integer values for ν i and η i , finding complete agreement.
For sample Mathematica code implementing this strategy, see appendix B.

Checks and examples
For checks, let us reproduce some known IBP results, using our main result eq.(6.2).
As another check, we can expand the numerator of eq.(3.1) for positive values of the index η 3 , and verify that each L can be expressed in terms of instances with η 3 = 0 as Indeed, testing eq.(6.12) for a wide range of integers {ν 1 , ν 2 , ν 3 , η 1 , η 2 , η 3 } by applying eq.(6.2) on both sides, we find perfect agreement.Going beyond published results, one could now use eq.(6.2) to generate a large table of sum-integral reductions, to facilitate future calculations and/or serve as a benchmark.We refrain from doing so, however, but list a couple of interesting cases that exhibit reductions to linear combinations of different master integrals, a feature not yet seen in the examples above.A few of these multi-term cases (at weight 4, 6 and 8) are On the other hand, from eqs. (6.4) and (6.5) some examples for non-positive indices read

Conclusions and outlook
We have presented an analytic solution for two-loop scalar massless bosonic vacuum sumintegrals L, with general integer values of propagator and numerator powers.Our key result presented in eq. ( 6.2) provides a mapping to analytically known one-loop sum-integrals, based on a similar decomposition that had been derived earlier for a class of closely related massive two-loop continuum integrals [33].While some fixed-index cases of such a two-loop onto one-loop sum-integral mapping had already been observed in the literature, our general-index proof exposes the mechanism behind this most helpful mathematical fact, and allows to eliminate all such sum-integrals from perturbative expansions, on an algebraic level.
For convenience, we have extended our result to be applicable also for non-positive propagator indices, essentially eliminating the need for tensor reduction when inverse propagators are present in the numerator.In practice, the reader might find the algorithmic implementation presented in appendix B useful, where all cases are combined.
Going beyond the scalar sum-integrals discussed here, one natural future generalization would be to consider tensors the numerator of eq.(3.1).If such tensors are fully symmetric in the µ i , one can trivially apply the decomposition presented in eqs.( 8), ( 16) of [27].For the general case, however, some additional work is required.This particular generalization might become interesting when evaluating four-loop vacuum sum-integrals, since then product-like structures ∼ L • L appear, with a numerator potentially containing scalar products between momenta from different L, necessitating a tensor decomposition.In practice, for fixed-index cases this can also be achieved by IBP relations, such that we postpone the analytic generalization to future work.
Other interesting avenues for future work would be generalizations to massless fermions (with and without chemical potential µ f ), and higher loops.Including (bosonic as well as fermionic) propagator masses deserves to be investigated as well, but seems to be a different game altogether.Note that not even the massive one-loop sum-integral is known analytically; what is known are only high-temperature expansions (as presented in textbooks, such as in the appendix of [2]) or practical integral representations and evaluations of some terms in their ε-expansion (see e.g.appendix A of [39] as well as [40]).
As for sum-integral factorization properties at higher loops, there seems to be little hope, if our observations of some particular fixed-index cases of higher-loop IBP reductions are of any guidance.One could hope, however, that for some special integral sectors (i.e.particular propagator subsets) such general simplifications might exist.

Figure 1 .
Figure 1.Different regions in the Matsubara double-sum, as used for deriving eq.(3.4).The regions in the first diagram on the right are one-dimensional sums only and can be immediately written as Riemann zeta functions, while we choose to map all eight regions of the last diagram to the lower-right region in the first quadrant, i.e. to double-sums of type n1>n2>0 .
combine to (a) cancel all Z (for a definition, see eq. (A.7)) in the last two terms of eq.(3.4),(b) cancel all ζ(i, j) (see eq. (A.2)) in the sum of all four terms of the last line of eq.(3.4), (c) cancel all remaining single ζ(i) in eq.(3.4) to (d) leave us with the expected tadpole reductions that contain products ζ(i) ζ(j) containing only ζ(n even −d) allowing us to map onto the 1-loop sum-integrals of eq.(2.4),