Integrating by parts at ﬁnite density

: Both nonzero temperature and chemical potentials break the Lorentz symmetry present in vacuum quantum ﬁeld theory by singling out the rest frame of the heat bath. This leads to complications in the application of thermal perturbation theory, including the appearance of novel infrared divergences in loop integrals and an apparent absence of four-dimensional integration-by-parts (IBP) identities, vital for high-order computations. Here, we propose a new strategy that enables the use of IBP techniques in the evaluation of Feynman integrals, in particular vacuum or bubble diagrams, in the limit of vanishing temperature T but nonzero chemical potentials µ . The central elements of the new setup include a contour representation for the temporal momentum integral, the use of a small but nonzero T as an IR regulator, and the systematic application of both temporal and spatial diﬀerential operators in the generation of linear relations among the loop integrals of interest. The relations we derive contain novel inhomogeneous terms featuring diﬀerenti-ated Fermi-Dirac distribution functions, which severely complicate calculations at nonzero temperature, but are shown to reduce to solvable lower-dimensional objects as T tends to zero. Pedagogical example computations are kept at the one-and two-loop levels, but the application of the new method to higher-order calculations is discussed in some detail.


Introduction
Perturbation theory is undoubtedly the single most powerful technique for making quantitative predictions in quantum field theory [1][2][3].Importantly, its application is not limited to problems in vacuum, but perturbative methods are equally applicable to the study of extended systems in and even out of thermal equilibrium [4][5][6][7].In the latter context, the physical systems of relevance range from the early Universe to the extremely dense cores of neutron stars and the hot fireball created in heavy-ion collisions.In such applications, Quantum Chromodynamics (QCD) has proven a particularly challenging theory [8] due to a combination of complicated phenomenological problems and the slow convergence of weak-coupling expansions.For this reason, thermal perturbation theory calculations in QCD need to be pushed to particularly high loop orders, and, e.g., the equation of state (EoS) of deconfined QCD matter is currently known to partial four-loop level both at high temperatures T and small or vanishing chemical potentials µ [9][10][11][12][13][14] as well as at high µ and small or vanishing T [15][16][17][18].
In high-loop-order perturbative calculations, be that in vacuum or a thermal setting, various methods of automation become indispensable.In vacuum, one of the most crucial tools enabling, e.g., the determination of the five-loop running coupling constant of QCD [19,20], is the integration-by-parts (IBP) technique [21,22].It allows for the derivation of linear relations between different Feynman integrals, reducing the number of independent master integrals in need of explicit analytic or numerical evaluation.These identities are typically generated by multiplying the integrand of a Feynman integral, schematically P g(P ), by a momentum-dependent function f ν (P ) and using the fact that integrals over total derivatives vanish, i.e.P ∂ ν f ν (P )g(P ) = 0.
Unfortunately, the derivation of IBP relations relies heavily on Lorentz symmetry, which is broken by both a nonzero T and µ.In the imaginary-time formalism of thermal field theory, applicable in thermal equilibrium, the spacetime metric becomes Euclidean with integrals over the temporal momentum component being replaced by discrete sums over the so-called Matsubara frequencies [23].The latter take the form p 0 = ω B/F n with ω B n = 2nπT for bosons and ω F n = (2n + 1)πT + iµ for fermions, n ∈ Z, where an imaginary shift appears for nonzero chemical potentials.This implies that total derivatives in the temporal direction no longer lead to vanishing momentum integrals, which prevents the closing of traditional IBP relations.In existing literature, this issue has typically been resolved by only considering spatial derivative operators and integrations in the derivation of IBP relations [24,25], which however severely limits their power in practical computations. 1n a largely unrelated development, it was recently pointed out in [26] that even the simplest loop calculations involving chemical potentials require particular care in the low-temperature limit.A naive application of the residue theorem may lead to the missing of important contributions to physical quantities.This issue appears whenever fermionic propagators are raised to (integer) powers higher than 1 and can be most conveniently circumvented through a consistent use of a contour-integral formulation for the temporal integral, introduced in [26] and reviewed in sec. 2 below.In the strict T = 0 limit, this procedure gives rise to δ-function terms easily missed in standard residue calculations that utilize a single linear contour along the real axis, 2 which also highlights the fact that thermal field theory calculations performed at T = 0, µ = 0 should always be thought of as the T → 0 limit of a finite-temperature computation.As detailed in [26], even if infinitesimally small, the parameter T plays an important regulatory role in loop integrals.
Motivated by the phenomenological need to determine the EoS of cold and ultradense quark matter, stemming from attempts to determine the EoS of neutron-star matter in a model-agnostic manner (see e.g., [28][29][30] for recent studies), our aim in this paper is to generalize IBP methods to be at least partially available in the limit of nonzero chemical potentials but vanishing temperatures.Crucially, we are interested in IBPs in all D ≡ 4−2ǫ dimensions, and to this end insist on including also temporal derivatives in the operators generating linear relations.As we will detail below, the use of the contourintegral prescription of [26] will give rise to extra boundary terms or inhomogeneities in IBP relations that often feature derivatives of the Fermi-Dirac distribution function.While highly problematic at nonzero temperatures, in the strict T → 0 limit these derivative terms simplify considerably, leading to more easily computable entities that close the extended IBP relations although some differences to standard vacuum IBP relations are seen to arise. 3  In the practical derivation of the new IBP relations, our strategy is to use temporal derivatives to write the boundary terms in a form, where ideally all fermionic distribution functions have been differentiated, while all other terms are of a form familiar from vacuum-type IBP relations.While the latter can be further simplified using spatial IBP relations as discussed in sec.3 below, the integrals with differentiated distribution functions need to be evaluated explicitly.The techniques needed in such computations are outlined in sec.2, where we also introduce our notation and general formalism.In secs. 4 and 5, we then provide explicit examples of this procedure at the one-and two-loop levels, wherein we also describe the required regularization due to the complex nature of fermionic propagators at nonzero µ [26].As a byproduct of our calculations, we end up generalizing the factorization of the vacuum sunset diagram with collinear scales to finite µ, motivated by an observed factorization from our IBP relations; this topic is discussed in detail in appendix E. The final objective of our program is to create an automatable algorithm that, in combination with other tools such as the cutting rules of [31], will facilitate computations at high orders of perturbation theory.Although this important extension of our framework is largely left to future work, some related aspects are discussed in the concluding sec.6.

Formalism and setup
We work in the imaginary-time formalism of thermal field theory.Here, the starting point of any computation at nonzero temperature T involves discrete sums over Matsubara frequencies [23], amounting to ω B n = 2nπT for bosons and ω F n = (2n + 1)πT + iµ for fermions at a nonzero chemical potential µ, with n ∈ Z.As demonstrated in numerous textbooks [4][5][6], sums over these discrete frequencies can be easily converted to continuous integrals using the analyticity properties of the Bose-Einstein and Fermi-Dirac distribution functions, Concretely, we may write the (sum-)integration measures appearing in thermal Feynman integrals in the forms where P = (p 0 , p) stands for a momentum in D = d + 1 dimensions, p is a d-dimensional spatial vector, and p ≡ |p|.Sum-integrals are denoted by the shorthand Σ P = T pn p , where p n = ω B/F n and curly brackets are used to indicate the fermionic nature of the loop momentum.Above, we have also introduced a shorthand notation for the (d + 1)-dimensional integral and its temporal integration contours where the integration contour runs (anti-)clockwise for fermions (bosons) around the real axis of the complex p 0 -plane.For the dimensionally regulated (MS) spatial integrals in d = 3 − 2ǫ dimensions, we finally introduce the measure p ≡ e γ E Λ2 4π where Λ is the corresponding renormalization scale and γ E is the Euler-Mascheroni constant.Unless otherwise stated, we will provide results for general d dimensions whenever possible.Below, we will briefly discuss two technical issues that are both frequently used in our forthcoming presentation.The first issue has to do with subtleties in taking the zero-temperature limit in the above (sum-)integration measures; see sec.2.1.The second introduces two convenient deformations of the fermionic p 0 -integration contour defined above, which turn out to greatly simplify practical calculations; see sec.2.2.

Small-temperature limit and regularization
As discussed in detail in [26], taking the zero-temperature limit in the integrals defined above can be surprisingly non-trivial, and in particular, changing the order of the temporal and spatial integrations must be handled with care.To prepare for these subtleties, we identify the real-valued parameters x ≡ Re(p 0 ) and y ≡ Im(p 0 ) − µ for the Fermi-Dirac distribution as well as x ≡ Re(p 0 ) and y ≡ Im(p 0 ) for the Bose-Einstein distribution, whereby the low-temperature limits of the corresponding distribution functions become (2.6) In this result, the bosonic case is always shown in parenthesis, and the chosen notation for the two functions defined on the lower line corresponds to the respective behaviors of the real and imaginary parts of the distribution functions at small but non-vanishing T : lim The derivative of the real-valued fermionic distribution function finally gives rise to a well-known nascent delta distribution, referred to as the δ-sequence below: (2.7) Similarly, the complex generalizations of the differentiated Fermi-Dirac and Bose-Einstein distributions (often referred to as primed distribution functions) yield where we introduced two more functions with the respective limits of lim T →0 δ F(B) T (y) = δ(y) and lim T →0 0 F(B) T 2 = 0.The simplest example of an integrand, for which the use of T as a regulatory parameter (through the above relations) becomes important 5 , consists of a massless fermionic propagator raised to an arbitrary real-valued power α.Here, a naive way to proceed would be to first drop the lower part of the p 0contour in eq.(2.3) due to its exponential suppression at small T and then to set n (2.9) This strategy indeed yields a p 0 -integral that can be easily evaluated using the residue theorem, but the result turns out to be the physically correct one only for 1/2 < α ≤ 1; see appendix A.1 and [26].For α > 1, an explicit divergence occurs at p = µ, p 0 ∼ 0 (2.10) 4 Upper indices in eq.(2.6) indicate the particular θ-function and 0-sequences defined by this equation. 5Keeping T nonzero prevents the temporal components of fermionic momenta p f 0 = ω F n = (2n + 1)πT + iµ from vanishing, which in turn protects integrals with high powers of fermionic propagators from problematic divergences from the region p = µ, Re(p0) = 0.
which explains the fact that different integration orders are seen to yield differing results (2.11) As first pointed out in [26] (see also app.A.1), this discrepancy is correctly treated through the Fermi-Dirac distribution and its derivative, which induces an additional δ-function-type contribution when p 0 crosses the imaginary axis along the contour of eq.(2.3).Similar issues are absent for purely bosonic integrals, given the lack of scales like chemical potentials.However, graphs mixing bosonic and fermionic momenta may again feature nontrivial low-temperature limits.Like in the fermionic case, keeping T small but nonzero in eqs.(2.6) and (2.8) and utilizing (simplifications of) the bosonic δ-sequences ensure that the order of the temporal and spatial integrations are always interchangeable.

Contour deformations
To no surprise, the technical complications that we encountered in the T → 0 limit of the one-loop integral (2.9) and that culminated in the isolation of a novel divergence in eqs.(2.10) and (2.11) continue to be present at higher loop orders.Although one may often bypass these issues by analytically continuing results derived for convergent values of parameters such as the dimension d and the exponent α above, in the application of IBP relations we will again encounter integrals containing differentiated Fermi-Dirac distribution functions that give rise to extra contributions from the momentum region that makes the argument of the n ′ F vanish, namely p = µ, p 0 ∼ 0. It turns out that for the two classes of p 0 -integrals featuring either undifferentiated or differentiated n F functions the optimal computational strategies differ somewhat and in particular feature slightly different ways of deforming the original p 0 -contour (2.4).Below, we briefly introduce these two deformations that are respectively summarized in fig. 1

(left) and (right).
The former class of integrals is initially characterized by such parameter values that allow one to take the step-function limit of the Fermi-Dirac distribution function without introducing new divergences, leading to special function solutions analytically continued to the relevant parameter values.Akin to eq. (2.9), it is then sufficient to only consider the first line integral in eq.(2.4), and η can furthermore be taken to 0. The practical evaluation of such integrals can be simplified by deforming the remaining integral in a fashion illustrated in fig. 1 (left) and detailed in [26].Using the Cauchy theorem, we can namely replace the infinite integral (L a , L b ) parallel to the real axis by three separate line integrals (L 1 , L 2 , L 3 ), of which L 2 and L 3 are of finite length and run along the imaginary axis while L 1 runs along the real axis, giving 12) The validity of this construction, detailed in appendix A.1, requires that the integrated function f (p 0 ) be regular inside the closed Cauchy contour.This can typically be achieved for certain ranges of Figure 1: Left: deformation of the line integral (2.12) for integrals involving n F -terms, i.e. (L a , L b ) to (L 1 , L 2 , L 3 ).The contours L a and L b indicate the original line integral parallel to the real axis.By mapping L a from (−∞ + iµ, iµ) to (−iµ, ∞ − iµ) and changing the sign of the integrand f (p 0 ) → −f (−p 0 ), the rectangular contour of interest can be closed, while the dashed contours can be neglected.In turn, the integral is evaluated in terms of the thick blue lines L 1 , L 2 as well as L 1 , L 3 .The lines L 2 and L 3 are oriented in opposite direction, while L 1 is used twice in this construction.Right: deformation of the fermionic contour (2.3) for integrals involving n ′ F -terms, i.e. (C a , C b ) to (C u , C l ).The initial contour, which is closed by two vertical segments at Re(p 0 ) = ±∞, is first split into two semi-infinite closed contours by introducing the vertical segment C u,1 and its inverse.The left of these closed contours is then mapped to the fourth quadrant of the complex plane upon setting p 0 → −p 0 , whereby the inverse of C u,1 becomes C l,1 and we end up with the final blue contours C u and C l .
propagator exponents, resulting in expressions that can be later analytically continued to the divergent parameter regions.
For the second class of integrals, a different strategy involving the full fermionic contour of eq.(2.3) is needed, owing to the sensitivity of the differentiated Fermi-Dirac distribution to those momentum regions where its argument vanishes.This can be understood by recognizing that eq.(2.8) describes a two-sided δ-distribution, which results in both line integrals of the fermionic contour (2.4) contributing to the results at T → 0. Assuming, as usual, that f (p 0 ) is regular within the fermionic contour, we can again split the original integrals, denoted now by C a and C b , at the imaginary axis as shown in fig. 1 (right).The contour on the left-hand side of the imaginary axis is then mapped to the fourth quadrant of the complex plane, so that (C a + C b ) → (C u + C l ) as detailed in the figure .Considering now the n ′ F -terms integrated over the contours C u and C l , we obtain where the only nonzero contributions originate from C u,1 and C l,1 of fig. 1 (right).The leading low-temperature behavior of the fermionic contour integral then becomes which demonstrates that the only nonzero contributions to the original integral originate from the neighborhood of p 0 = iµ of eq.(2.4).The deformed contours defined here greatly simplify analytic computations involving the δ-sequences but are nevertheless able to capture the novel thermal corrections from the differentiated distribution functions.This resembles, to some extent, the role of the finite line integrals in eq.(2.12).

Integration-by-parts relations
Moving on from the details of integration contours to the derivation of the desired IBP relations, let us first specify the precise form of the vacuum-type Feynman integrals we will be studying in this paper.To this end, we abbreviate the Fermi-Dirac distribution functions as and assume a scalarized numerator structure in the Feynman integrals to be considered.A generic vacuum-type Feynman integral then takes the form where we again denote momenta in d + 1 dimensions by P = (p 0 , p).In this expression, N f ℓ stands for the number of fermionic and N b ℓ for the number of bosonic loops, with P k = (p k,0 , p k ) denoting fermionic and Q l bosonic loop momenta.Finally, the momenta R (no indices) are linear combinations of the P k and Q l , picked from the sets R ± ≥2 (P k , Q l ). 6s two explicit examples of the (scalarized) integrals to be considered, the one-loop fermion bubble and the two-loop fermionic sunset are defined as where we henceforth use calligraphic letters for in-medium and Latin letters for d-dimensional T = 0 loop integrals and also identify I α (µ, T ) ≡ I 0 α (µ, T ) as well as S α 1 α 2 α 3 (µ, T ) ≡ S 00 α 1 α 2 α 3 (µ, T ) .In the graphical notation, wiggly lines always stand for bosonic and directed solid lines for fermionic propagators.We also introduced the corresponding d-dimensional vacuum integrals for the one-loop bubble and two-loop sunset while deferring special mass signatures of the latter to appendix D.1.
In practice, both spatial and temporal IBP relations are derived starting from integral expressions that vanish as total derivatives.Assuming that the integration orders are interchangeable, for one loop momentum, P µ , such a relation can be written as where the integrand f (P ) typically takes the form of a product of individual propagators and we have introduced a notation for the differential operator ∂ ∂P µ • P ν acting on an integral.The last term on the last line of eq.(3.7) is seen to contain a derivative of the Fermi-Dirac distribution, which is specific to a thermal setting.The explicit distribution functions are essential for ensuring that the total derivatives vanish and that integration orders are interchangeable.
Considering next two distinct loop momenta P and Q, the simplest differential operators that can be expected to give rise to non-trivial IBP relations are bilinear in both momenta and derivatives.They form two disjoint classes that are either diagonal or off-diagonal7 in Lorentz indices, When acting on the integrand apart from the distribution functions, the two classes are self-contained, such that integrals generated by an operator of a given class can be algebraically related to integrals generated by other operators in the same class.This property can be seen in the one-loop expression where the emergent quadratic numerator structure p 2 0 has been expressed in terms of a scalar multiplier and a total spatial derivative.
So far, only spatial differential operators have been considered in the derivation of IBP relations [24,25] for systems with broken Lorentz symmetry.The extension we propose here is to generalize the IBP relations to the full (d + 1)-dimensional spacetime, taking advantage of natural simplifications such as that present in eq.(3.9).A pedagogical example of this procedure at the one-loop level is obtained by acting on I s α (µ, T ) with both spatial and temporal derivatives.This leads to where on the last line we combined spatial and temporal total derivatives into the (d + 1)-dimensional bilinear operator In the relations (ibp.1a)-(ibp.1c),each term notably involves either one differentiated or one nondifferentiated distribution function.Proceeding to higher loop orders, the classification of integrals appearing in our generalized IBP identities will naturally become slightly more complicated, but all integrals nevertheless still fall into one of the following three disjoint classes: type A: integrals with no differentiated distribution functions, type B: integrals with only differentiated distribution functions, and type C: integrals with at least one differentiated and one non-differentiated distribution function.
Conventional spatial IBP identities such as (ibp.1a)only involve integrals of type A, which are also familiar from vacuum computations [21,22].There they typically reduce the numerator and denominator powers in type A integrals while introducing rational functions of the dimension d.In contrast, the novel terms of types B and C are generated by temporal derivatives in the thermal context, where they play the role of inhomogeneous terms closing the IBP relations between type A integrals.This can be seen already in the one-loop relation (ibp.1c),where the presence of the last term allows writing I α (µ, T ) in terms of an integral of reduced dimensionality in the zero-temperature limit due to the δ-sequence of eq.(2.7), thus constituting a new representation for the original integral.
At higher loop orders, the utility of the novel IBP relations similarly derives from a qualitative hierarchy between the computational workloads associated with evaluating type A, B, and C integrals: type A > type C ≫ type B. This is due to the appearance of differentiated distribution functions in the type C and B integrals, which according to the δ-sequence of eq.(2.7) amounts to a reduction in their dimensionality. 8As we will demonstrate in the following two sections of the article, results for type B integrals are typically immediately available through a straightforward analytic continuation of known massive vacuum (T = µ = 0) integrals, but even type C integrals are normally dramatically simpler to evaluate than the original type A entities.Accordingly, our general strategy will focus on deriving relations between type A and B terms, i.e. integrals involving either only differentiated or non-differentiated distributions functions, with type C terms remaining present only if absolutely necessary.

One-loop integrals: new strategy in action
Having laid out our general strategy above, we will now take a closer look at explicit IBP calculations at the one-loop level, where we first return to the well-known integral (3.3) with s = 0 (for s ∈ N, see eq. (A.6)).It was shown in [26] that in the limit T → 0 and for any α ∈ R, the integral takes the form which is nothing but the T → 0 limit of the I 0 α function defined in [10], generalized to an arbitrary spatial dimensionality d.The integral converges for parameter values of α and d inside the triangleshaped region {2α − d > 1 | 1 2 < α < 1 ∧ d > 0}, but for strictly vanishing T features a d-independent novel divergence in the neighborhood of p = µ for α > 1 that requires analytic continuation from the convergent region.The result in eq.(4.1) is easiest derived starting from the spatial integral and utilizing a complex-valued generalization of Euler's beta function, 9 which ensures a consistent treatment of the problematic neighborhood in the low-temperature limit.A generalization of this result to non-vanishing masses is detailed in appendix B with additional context found in sec.III of [26].Further consistency checks are relegated to appendix C, including the verification of the fact that acting on the master integral eq.(4.1) with spatial total derivatives ∂ ∂p i • p i leads to a vanishing result as it should.
As alluded to in the previous section, acting on thermal integrals with temporal total derivative operators may give rise to non-trivial linear relations, where terms featuring differentiated distribution functions are expected to lead to simplifications in the low-temperature limit.To demonstrate this at the simplest possible level, we focus on the identity (ibp.1c), which clearly allows the evaluation of the master integral I α (µ) through a new representation To evaluate the right-hand side of this relation, we first apply a complex-valued generalization of the integral form of Euler's beta function, resulting in eq.(2.8) Subsequently, we utilize the contour prescription of fig. 1 (right), where nonzero contributions arise from the two finite line segments C u,1 and C l,1 .Since these two terms are related to each other via complex conjugation, we easily obtain lim Finally, to ensure that the monomial structure resulting from the final p 0 -integration is well-defined for arbitrary non-integer powers α, we regulate the imaginary unit inside the integral by defining with κ > 0 an infinitesimal positive real number.This allows to extract a result in terms of trigonometric functions, swiftly leading to the same result reported in eq.(4.1) upon taking κ → 0.
As demonstrated by the above example, the appearance of a differentiated distribution function reduces the number of integrations by one in the T → 0 limit, while the remaining ones correspond to complex-valued generalizations of standard d-dimensional integrals.It is often operationally useful to first perform the temporal integral but leave the complex-valued substitution (with the regulated imaginary unit Π + ) to be performed after the spatial integrations.In the above example, this would imply writing lim where the spatial integral can be identified as the standard one-loop vacuum integral from eq. (3.5).
To conclude this section, a few general remarks are in order.First, it should be apparent from above that the procedure applied in eq.(4.6) will greatly simplify the treatment of the inhomogeneous (type B) terms at higher loop orders, which should optimally contain the maximal number of differentiated distribution functions.Second, we have relegated several additional example calculations to appendix C; these include consistency checks of the one-loop IBP relation (4.2), an explicit computation involving linear p 0 -structures within the integrand numerator, and a brief discussion of the use of parametric µ-differentiation.

Proceeding to the two-loop level and beyond
The value of the proposed novel IBP machinery only becomes apparent in multi-loop computations.While this paper is mostly devoted to developing the new formalism, we will therefore next discuss its application at the lowest nontrivial order of perturbation theory, i.e. two loops, which amounts to the leading correction to the non-interacting limit for the QCD pressure.This order is somewhat special, as it allows the derivation of nontrivial IBP relations while keeping the practical computations at a pedagogical and easily tractable level.This is particularly true for integrals with unit numerators but sufficiently general denominators, so that the result does not trivially factorize into one-loop structures. 10or the aforementioned reasons, we shall concentrate on the fermionic sunset (cf.eq.(3.4)) where S 00 α 1 α 2 α 3 ≡ S α 1 α 2 α 3 and the non-calligraphic S α 1 α 2 α 3 denoted the corresponding d-dimensional vacuum sunset from eq. (3.6).The denominator exponents {α 1 , α 2 , α 3 } can be chosen such that the p 0 -and q 0 -integrals are regular both in the ultraviolet (UV) and at possible individual poles.As detailed in sec. 2 and appendix A, this allows for using the Cauchy theorem (see e.g.eq.(2.12)) to analytically continue the results to all integer-valued exponents.In addition, it allows choosing the integration order at will even in the low-temperature limit and conveniently isolates the novel boundary contributions to finite line integrals along the imaginary axis as summarized in fig. 1 (left).
Finally, we reiterate that the main goal of our new IBP program is to find relations between Feynman integrals that are of the vacuum-type A in the presence of additional non-vacuum type B terms that are needed to close the relations.The latter terms include a maximal number of differentiated distribution functions and are computationally less complex as demonstrated in the previous section at the oneloop level.At the two-loop level, scrutinized in the present section, we begin our discussion from these maximally primed terms in sec.5.1, derive and solve novel IBP relations in sec.5.2, and finally, dwell on the closely related factorization of d-dimensional two-loop vacuum integrals in the subsequent appendix E. While the implementation of the computational methods is explicitly discussed only at the two-loop level, all results are at least in principle straightforwardly generalizable to an arbitrary loop order.

Integrals with only differentiated distribution functions
Let us first examine the evaluation of the boundary or inhomogeneous terms containing only primed distribution functions n ′ F .Using the contour prescription of eq.(2.14), we can perform the Π + substitutions akin to eq. (4.6) independently for the T → 0 limit of each temporal integral, which generates 2 n independent terms from the contour integrals of an n-loop diagram.The sign of each term is determined by the even/odd symmetry of the integrand as indicated by eq.(2.14).Below, we detail this procedure at the two-loop level focusing on the sunset integral of eq. ( 5.1).
Acting on S α 1 α 2 α 3 with two temporal derivatives gives rise to a term with two primed distribution functions, i.e. an integral of the maximally primed type B. For the purpose of abbreviating such expressions, we introduce n F -differentiating operators D k that only act on the distribution functions of the corresponding temporal momentum via Consistently with the above discussion, we find altogether 2 2 = 4 separate terms in the lowtemperature limit, totaling Each of the terms can be mapped into a form where the propagator masses of underlying vacuum sunset S α 1 α 2 α 3 (m 1 , m 2 , m 3 ) obey the linear relation m 1 + m 2 = m 3 .This class of massive Feynman integrals is typically referred to as collinear [32].For their evaluation, we use standard Feynman parametrization and carefully regulate the vanishing bosonic mass scale11 related to (p 0 ± q 0 ) = O(κ) with the κ-regulator introduced in eq.(4.5).One curious detail of this procedure is that for positive integer values of the exponent α 3 , associated to the bosonic scale, each arising hypergeometric integral from eq. (D.7) can be shown to factorize into a product of one-loop vacuum integrals.This is apparent for all four integrals in eq. ( 5.3) wherein the complex-valued scales p 0 and q 0 allow for extracting contributions proportional to p 2 0 − q 2 0 = O(κ).This factorization property at the two-loop level will be addressed using a closed-form solution of the collinear vacuum sunset in appendix E and at the hypergeometric function level in appendix D.1.
Next, we demonstrate in detail, how the above maximally primed integrals of type B can be computed for S s 1 s 2 α 1 α 2 α 3 using the strategy established in sec.4. We do this by allowing the power of the temporal momentum components s 1 + s 2 in the numerator to take arbitrary even or odd values.To this end, we first write eq. ( 5.3) in the form where we have only used the fact that s 1 + s 2 is an integer but note in addition that the term corresponding to the imaginary (real) part here clearly vanishes for even (odd) values of this parameter.
Finally, the corresponding d-dimensional vacuum sunset integral is evaluated in eq.(D.10). 12or concreteness, let us next inspect the simplest nontrivial sunset S 111 .After directly evaluating the d-dimensional integral using (D.15), we observe that eq. ( 5.4) can be recast as a product of two one-loop integrals.The result reads We note that in this expression, S 111 (µΠ + , µΠ + , 0) resembles the standard result for the vacuum sunset with two real-valued mass scales [33], continued to complex scales, 13 while S 111 (µΠ + , µ/Π + , 0) needs to be evaluated with eq.(D.10).It is also worth pointing out that akin to eq. (5.5), all type B terms originating from eq. ( 5.1) are seen to naturally factorize into one-loop master integrals.For more details on these considerations and an independent calculation applying the factorization of vacuum bubble diagrams [32], we refer the reader to appendix D.1.
one another, so that e.g. the replacement of α3 ↔ α2 and q0 − p0 ↔ q0 yields the spatial integral Sα 1 α 3 α 2 (p0, q0 − p0, q0) but does not change the value of the original Sα 1 α 2 α 3 (p0, q0, q0 − p0).The spatial integrals in eq.(5.4) are written to make the collinearity of the scales explicit (e.g.p0 + (q0 − p0) = q0) which allows the direct application of the results of [32], further discussed in appendix E. 13 The mass scales appearing in the propagators can be kept general in such a computation, with the only necessary assumption being that the squared mass scales are not strictly negative (i.e. that they are either complex-valued or positive and real).As shown in [26], this allows the necessary extraction of beta functions on the complex plane [26].

Integration-by-parts at two loops
Armed with the knowledge of the explicit expressions of the maximally primed terms introduced above, we will next inspect the derivation of the (d+ 1)-dimensional IBP relations at the two-loop level.Here, we will again concentrate on the sunset integrals of eq. ( 5.1), with our aim being to relate them to simpler structures.To achieve this, we find it convenient to study S 111 in a slightly different manner from the strategy described in sec.3, which involves an IBP relation initially containing integrals of type C, containing both differentiated and non-differentiated distribution functions.While this leads to a simple factorizing result in this particular case, with more complicated propagator structures, we strongly recommend the use of IBP relations containing only terms of type A and B.
The first nontrivial IBP relation arises upon taking a diagonal total derivative along the temporal direction, which produces Applying next a diagonal total derivative in the spatial direction on the last term on the right-hand side, it can be further recast into where the spatial derivative term clearly vanishes.Following the same strategy as for the diagonal temporal p 0 -derivatives, all other purely temporal combinations from the bilinear set of eq.(3.8) generate a system of in total n 2 linear relations at the n-loop level.Together with the n 2 spatial IBP relations (ibp.2c)-(ibp.2d) in appendix D.4, this gives rise to the following (d + 1)-dimensional IBP relations for the two-loop fermionic sunset which are written in a more compact operator form, omitting S s 1 s 2 α 1 α 2 α 3 onto which each operator acts.The two remaining bilinear combinations corresponding to P ↔ Q can be found through substitutions 1 ↔ 2.
Here, we have employed standard IBP notation for the raising and lowering operators [34] for the denominator and {p 0 , q 0 } numerator powers via ) (5.9) The relations (ibp.2a)-(ibp.2b)represent the two-loop equivalent of the one-loop (d+1)-dimensional IBP relation (ibp.1c).Independent of the number of loops, the system of IBP equations is always linear, underdetermined, and inhomogeneous, as the finite-µ relations are characterized by their righthand sides containing the D k differential operators.While no general closed-form solution of such infinite systems of equations can be found with current methods, a typical approach to find a solution is to attempt solving for a finite number of master integrals [35] using a finite set of starting integrals in the parameter space of {α 1 , α 3 , α 3 ; s 1 , s 2 } for S s 1 s 2 α 1 α 2 α 3 and its n ′ F -version.The solution is then found via Gaussian elimination implemented in a Laporta-type algorithm [25].While this type of a systematic approach, implemented using suitable programs for symbolic manipulation (e.g.FORM [36]), typically quickly becomes indispensable at higher loop orders, below we approach the two-loop problem iteratively by hand.
For concreteness, let us now focus on the specific choice of indices α 1 = α 2 = α 3 = 1, s 1 = s 2 = 0, which constitutes the S 111 integral.In this case, the IBP relation (ibp.2a)gives rise to an aesthetically pleasing linear dependence between differentiated and non-differentiated integrals after applying the P ↔ Q symmetry between the loop momenta, (5.10) Next, we linearly combine IBP relations by operating on the above expression for S 111 with a diagonal total derivative in both the P µ -and (P µ − Q µ )-directions.This results in the relation where the term 3 + (2 − − 1 − )S 111 = S 102 − S 012 can be shown to cancel on the integrand level upon a relabeling of loop momenta.The integrals S 012 and S 102 also vanish individually in the T → 0 limit, since they factorize into a fermionic and a vanishing bosonic one-loop integral, as detailed in appendix D.2.By explicitly writing the integrals in the above result, we arrive at the relation where on the right-hand side we identify a sum of three terms.Among them, the first integral clearly factorizes into two one-loop expressions, the second contains a shifted bosonic propagator with momentum P − Q, and the third is a linear combination of type B terms in eq.(D.17) which vanishes in the T → 0 limit.Upon closer inspection, detailed in appendix D.2, even the second integral can be shown to vanish as T → 0, since it contains a µ-independent loop integral.
By removing the vanishing integrals from eq. (5.12) and using the one-loop (d + 1)-dimensional IBP relation (ibp.1c), the low-temperature limit of the sunset integral is seen to simplify to This result is in agreement with computations of the low-temperature pressure of QCD (see e.g.appendix B of [10]) and can moreover be computed directly using the cutting rules of [31], as we demonstrate in appendix D.3.As far as we are aware, the factorization formula is, however, new.It is notably of different structure than the factorizing integrals of odd-valued numerator exponents in eq.(D.29), which was derived using a set of spatial IBP operators listed in eqs.(ibp.2c)-(ibp.2d) of appendix D.4.As a consequence, we can relate the factorizing one-loop integrals with odd-valued numerator exponents to the even-valued ones.This is technically not a linear relation among one-loop integrals and curiously only gets generated at the two-loop level.At low temperatures, we find the relation which represents the T = 0 limit of a more complicated relation valid at all T and µ.Note that in the opposing limit of µ = 0 but T = 0, integrals with odd-valued numerator powers such as the left-hand side of (5.14) vanish, implying that the relation derived here will necessarily obtain corrections at nonzero T .While the above example led to a simple and aesthetically pleasing result, it was not perfectly in line with the strategy we outlined in sec.3, where we emphasized the importance of the maximally primed type B terms.A more systematic and more easily generalizable strategy in line with this idea is to consider exclusively differential operators that contain temporal derivatives with respect to all fermionic loop momenta.The vanishing of the total derivatives then leads to relations such as where only terms of type A and B are present.
To demonstrate the latter strategy in action, let us again consider the general fermionic sunset S α 1 α 2 α 3 with no numerators.While the full (d + 1)-dimensional diagonal IBP relations (ibp.2a)-(ibp.2b)already give rise to non-trivial relations, full IBP simplification is only possible via their combination with either the temporal or spatial IBP relations generated from the bilinear set (3.8).By combining two consecutive diagonal d + 1 and temporal derivatives, we indeed obtain which is the only relation we obtain given its symmetry under the exchange of P ↔ Q.After using the relations (ibp.2a) and the temporal version of (ibp.2b),we carefully expand all derivatives on the lefthand side of this equation and again combine the corresponding spatial and temporal total derivatives following eq.(ibp.1c).After multiple steps of straightforward algebra, the resulting expression (5.16) simplifies to a recurrence relation between various sunset integrals S s 1 s 2 α 1 α 2 α 3 that can be dressed in the finite-T form 14∂ ∂P µ where and we have used the operator definitions of eqs.(5.2), (5.8) and (5.9).The right-hand side of this expression can be studied using the factorization formulae derived in appendix D.1 or the collinear integrand formulae given in [32] and their extensions discussed in appendix E.
As an interesting special case of the above result, we note that by setting once again α 1 = α 2 = α 3 = 1 and s 1 = s 2 = 0, the identity reduces to After removing the vanishing bosonic integral S 202 and inserting the result of eq.(2.14), the T → 0 limit of eq.(5.19) reads given again in terms of eq.(4.1).By independently studying the lengthy spatial IBP relation (D.28), we find a factorization into the one-loop integrals of eq.(3.3) and recover a non-trivial result for general µ and T which also agrees with the T = 0, µ = 0 result of [37] 15 and wherein we use the superscript b to denote the bosonic finite-temperature integral I b α (T ).It is also straightforward to take the T → 0 limit from eq. (5.21), giving which is in agreement with eq.(5.20).We emphasize, though, that the derivation leading to the more general spatial IBP relation requires a lengthy ansatz with 10 bilinear operators [see eq.(D.28)] as opposed to the more compact one in eq.(5.16).We find this a very promising observation from the point of view of the future evaluation T = 0 Feynman integrals at higher loop levels.

Conclusions and outlook
As recently pointed out in [38], to improve from the present uncertainties in the model-agnostic determination of the neutron-star-matter equation of state [28][29][30], the thermodynamic properties of cold and dense quark matter need to be known at unprecedented precision.In the absence of nonperturbative tools, perturbation theory plays a pronounced role in efforts to reach this goal, with the most prominent individual challenge being the extension of the weak-coupling expansion of the pressure of cold and dense unpaired quark matter.While the current state of the art in this problem lies at an incomplete O(α 3 s ) order [17,18,39], all qualitative issues hindering the completion of this full order have been resolved by now, so that only one numerical coefficient is lacking.This is the contribution of the hard momentum scale µ B to the O(α 3 s ) pressure, encoded in the sum of all four-loop vacuum, or bubble, diagrams of full QCD.
While conceptually a straightforward task, the practical evaluation of four-loop vacuum diagrams in a thermal setting presents an extremely complicated technical challenge.Only a handful of simple individual diagrams have been successfully completed by now [40][41][42], and they all correspond either to a scalar field theory or to very specific gauge-theory topologies that greatly simplify the evaluation of the diagrams.The particular case of dense fermionic matter in the T = 0 limit can, however, be seen to exhibit some simplifications, including the vanishing of diagrams with no fermion loops and the presence of a convenient computational tool, the so-called cutting rules of [31].It is, however, evident that a successful evaluation of all four-loop vacuum diagrams in the theory necessitates automated tools of computation, both in the identification and eventual evaluation of the master integrals.Given the complexity of the task at hand, even minor advances, perhaps leading to the successful evaluation of a handful of four-loop master integrals, would be extremely valuable.
In the paper at hand, we have taken first preliminary steps towards generalizing one of the most powerful methods of vacuum (T = µ = 0) perturbation theory to the context of nonzero chemical potentials: integration-by-parts, or IBP, techniques.They provide linear relations between master integrals, which are typically derived using the vanishing of (loop) integrals over total derivatives and lead to a reduction in the overall number of masters.In the thermal context, a problem can be seen to arise from additional boundary contributions that are generated due to the breaking of Lorentz invariance and are related to the special nature of the temporal direction.Our novel idea, presented in this article, is to generate the full (d + 1)-dimensional IBP identities starting from a contour integral formulation of thermal sums.This leads to the extra boundary terms being tractable in the small-T limit due to the simple limiting behaviors of the bosonic and fermionic distribution functions and their derivatives.Another key element of our approach is the use of the temperature T as a natural regulator for specific types of infrared divergences that are encountered in the strict T = 0 limit.
As part of our new framework, we introduce two deformations of the usual complex contour of fermionic p 0 -integrals: eq.(2.12) or fig. 1 (left) for convergent line integrals involving non-differentiated distribution functions, and eq.(2.14) or fig. 1 (right) for thermal corrections involving δ-function limits.Their use requires the order of various integrations to be interchangeable and leads to many simplifications as demonstrated in secs.3-5.For the simplicity of presentation, the example calculations we present here are restricted to one-and two-loop calculations involving massless propagators.The proposed method is, however, more general, and we anticipate its true utility to manifest at the three-and four-loop orders.
We also note that in practical computations, our approach can (and should) be complemented by purely spatial IBP relations, which have been studied in the past both in vacuum and in thermal systems lacking chemical potential.At the two-loop order in particular, the integrands of (d + 1)dimensional bubble diagrams with chemical potentials experience factorization in a manner similar to how d-dimensional real-valued vacuum bubbles can be expressed as linear combinations of products of their collinear mass scales, m 1 + m 2 = m 3 (see [32] and appendix E), with complex-valued temporal momenta now playing the role of masses.An explicit example of the application of spatial IBP relations is given in appendix D.4, with emphasis on how relations derived at µ = 0 may differ from the µ = 0 ones due to symmetries broken by a nonzero chemical potential.
Finally, we note that while the spatial operators of eq.(D.28) can be used to derive closed expressions for S 11s for s ∈ {1, 2}, their (d + 1)-dimensional counterparts in eqs.(5.11) and (5.16) are noticeably simpler.Furthermore, the explicit one-loop factorizations arising from temporal differentiations [eqs.(5.12)-(5.13)and (5.19)] differ from those obtained with spatial IBPs, cf.eqs.(D.29) and (D.30).These features imply that our new framework is capable of producing novel and useful results, a property that we expect to be particularly pronounced at higher loop orders.While formally not a closed group action with respect to differentiation, unlike the purely spatial differential operators, the temporal derivatives are seen to complement the existing IBP framework.

Future directions
To conclude our presentation, let us make a few remarks on the potential extensions and applications of our main results, including some details of calculations at higher perturbative orders.So far, our discussion has relied on the presence of only one physical scale (the chemical potential µ) in the integrals considered.While this special case often suffices for phenomenological applications in neutron-star physics, it should be noted that our generic method can be straightforwardly generalized to the presence of flavor-dependent chemical potentials, nonzero fermion (or boson) masses, and even to n-point functions, i.e. the presence of external momenta.The complexity of the novel δ-function integrals, however, rapidly increases with the number of scales present, which we demon-strate in appendix B.2 where the one-loop integrals of sec. 4 are generalized to the case of nonzero masses.While this exercise and its two-loop extension in eq.(E.7) can still be completed in a fairly straightforward manner, introducing external momenta with non-vanishing temporal components or considering topologies more complicated than the sunset one would invalidate all the computational tools relying on collinearity.
At higher loop orders, practical calculations will inevitably rely on a large-scale automation of the IBP reduction algorithm.While we have not yet performed extensive studies in this direction, we expect an automated approach to be fully compatible with our strategy to generate relations where all new thermal terms are of type B, i.e. contain a maximal number of differentiated distribution functions that reduce to delta functions as T → 0. While the calculations will surely be technically more demanding than in our treatment of the fermionic sunset integral in sec.5, we expect important simplifications to occur from the fact that the extra boundary terms have the form of complexified d-dimensional massive integrals that can be evaluated with the standard methods of vacuum field theory.
In an algorithmic implementation of the IBP approach, one faces an infinite, underdetermined system of inhomogeneous linear equations for the unknown integrals. 16In a practical approach, we are typically only interested in a finite number of master integrals, so instead of solving an infinite system at once, it is sufficient to only solve for a finite number of integrals that also constitute a finite linear system of unknowns.This system can in turn be solved via Gaussian elimination -or a variant of the so-called Laporta algorithm [25].For the fermionic sunset considered in sec.5, this implies applying both the diagonal (d + 1)-dimensional identities (ibp.2a)-(ibp.2b),the spatial identities (ibp.2c)-(ibp.2d),and their off-diagonal variations to a sufficiently large but finite set of integrals such that all initial Feynman diagrams are expressed in terms of the masters.The complexity of this process rapidly increases with the number of loops given that at the n-loop level, the IBP-system of vacuum integrals is generated by 2n 2 diagonal and 2n 2 off-diagonal identities.
As we want the master integrals to be as simple as possible, we propose an ordering prescription that penalizes more complicated integrals during the reduction.In contrast to the standard Laporta algorithm [25], and basing on our treatment of the fermionic sunset S s 1 s 2 α 1 α 2 α 3 (cf.eq.(5.1)), we propose the following ordering in a decreasing degree of complexity: (i) lowest number of primed distribution functions, (ii) highest total power of temporal momenta in the numerator (s 1 + s 2 ), (iii) highest total number of propagator powers in the denominator (α With this type of an ordering, it is conceivable to implement a generalized Laporta algorithm that respects the presence of the novel type B boundary terms at nonzero T and µ.A practical implementation of the algorithm is beyond the scope of this preparatory study, but we plan to tackle it in the near future as part of the evaluation of the four-loop vacuum diagrams of QCD in the limit of vanishing temperature and nonzero chemical potentials.

A. Applying the Cauchy theorem to loop integrals
This appendix collects additional details on the constraints related to the Cauchy theorem as it appears in secs.2-5.We also further detail the application of the Cauchy theorem at the one-loop level.

A.1. Single propagator
The simplest integral encountered at the one-loop level is given by eq.(4.1).As an example, we first revisit this integral, working in the parameter region where it converges: A subset of the above parametric conditions for convergence arises from both the UV behavior and novel µ-specific divergences occurring at p ∼ µ and p 0 → 0. To track down the regularity of the integral, we apply the low-temperature step function limit as in eq.(2.6) to only consider the upper line integral.After a change of variables to move to the real axis, we remove the regulator η, perform the d-dimensionally regularized integral, and then generalize the beta function on the complex plane following [26].As a result, we obtain the intermediate regulated expression On the last line, we have clearly separated the UV convergence condition 2α − d > 0 arising from Euler's Gamma function and d+ 1− 2α < 0 from the remaining integration.Conversely, when starting with the p 0 -integration, the UV convergence requires that α > 1 2 .The remaining parametric condition arises from eq.(2.10), which occurs at p = µ and p 0 ∼ 0 on the first line of eq.(A.2).Let us reiterate the equation and note that it only converges for α < 1.17 Combining these conditions, both the full integral in the region of eq.(A.1) and also the p 0 -integral converge.The latter allows us to use the Cauchy theorem in the computations.
It is worth noting that it would be possible to apply the Cauchy theorem already for the original contour integral of eq.(A.1).While this integral does not play a role in our computations involving δ-function entities, it is quite demonstrative, and to this end, we next reverse the order of integration starting with the temporal integral.To start, we, however, need to again remove the lower line integral and replace the distribution function by unity for the upper one.Since everything apart from the above-mentioned occupation-function limit is exponentially suppressed by e −βη (for η > 0, see eq. (2.6)), we can neglect subleading terms, giving the thermally leading integral the convergent structure Given that the above integral is convergent, the η-regulator in the line integral could be dropped, followed by changing p 0 → −p 0 in the first integral; this is the setup described in sec.2. As the integrand is even in p 0 , it does not visibly experience the sign change as described in eq.(2.12).To find the most suitable alternative representation for both of line integrals, we close the contours along the imaginary and real axes as in fig. 1 (left).Neglecting the finite segment at real positive infinity Re(p 0 ) = ∞ is possible due to the convergence for α > 1  2 .This way, we find the minimal integral expression The real line integral vanishes in dimensional regularization, and we were able to simplify the second line since all values p > µ explicitly yield imaginary values for the integral inside the brackets.While the result is formally limited to a small subset of dimensions and exponents, we can analytically continue it to almost the full parameter space, as suggested in eq.(4.1).At the one-loop level also another class of master integrals can appear depending if their temporal momentum numerator powers are even (σ = 0) or odd (σ = 1).By combining the Cauchy theorem and starting with the spatial integral, we obtain where I 0 α (µ) ≡ I α (µ).

A.2. Multiple propagators
As discussed in the main text, the Cauchy theorem approach generalizes to integrals with multiple propagators and loops.In this section, we indeed extract convergence constraints for multi-loop and multi-propagator structures similar to the one-loop case in eq.(A.1).To this end, we consider an integral with a (d + 1)-dimensional external momentum K = (k 0 , k) with a complex-valued temporal component, relevant to thermal field theory calculations: The novel problematic divergence should now only occur at |p − k| 2 ∼ [µ − Im(k 0 )] 2 for α 2 ≥ 1.Therefore, we can even deal with multiple propagators and analytically continue the result by applying the exponent regulation of appendix A.1.
Each loop order introduces one additional contour integral, akin to the two-loop level discussion in sec. 5.As at the one-loop level, we study the convergence properties of eq.(5.1), and write the strict (naive) low-temperature limits in terms of the upper line integrals of the contours: The first two (fermionic) propagators are similar in their behavior when p = µ for α 1 ≥ 1 and q = µ for α 2 ≥ 1.The third (bosonic) propagator on the other hand diverges in the hyperplane |p − q| = 0 for α 3 ≥ 1, similar to the other two cases but shifted from the position near the origin by the length of p 0 or q 0 depending on the chosen integration variable.Again, by neglecting these divergences, not only the integral becomes dependent on the integration order but also the dimensionally regularized loop momenta.Thus, the conditions α 1 , α 2 , α 3 < 1 are necessary (but not sufficient) to be able to apply the Cauchy theorem.Further conditions are needed to treat the UV behavior of both zero-components of the (d + 1)-loop momenta.For this purpose, we use the Feynman parametrization which is valid for α 1 , α 2 > 0. Parametrizing the second and third propagator in this fashion yields , where we apply the compact notation α {i} = j∈{i} α j .
The UV behavior of the integral is contained in the q 0 -integral if 1 2 − α 23 < 0. To consider this behavior explicitly, we change the integration variables according to u(q 0 ) = q 0 − xp 0 for fixed values of x and p 0 , and focus on the innermost integral The UV behavior of this expression can be studied easier by moving the line integral to run along the real axis which gives rise to a beta function as in our Cauchy theorem prescription.Hence, the computation reduces to evaluating where we used a complex-valued generalization of Euler's beta function to find the given representation (see e.g.sec.III of [26]).Extracting all of the elements contributing to the remaining p 0 -integration above, we need to consider another Feynman parametrization: The UV behavior of this integral can again be extracted after a reflection of the p 0 -integral to the real axis and the evaluation of the corresponding beta function.The resulting hypergeometric integral structure is well-defined for α 123 − 1 > 0, α i<j > 1 2 and α i > 0 for {i, j} ∈ {1, 2, 3}.These constraints are complemented by the above-mentioned conditions α i < 1, associated with the IR behavior of the propagator structure.
Applying total derivatives to any full expression such as eq.( 5.1) can organically generate different propagator powers.However, differentiation can not be successfully applied directly to formulae, in which the zero-temperature limit has been previously taken such that the distribution function is fully removed; see eq. (A.5).Indeed, differentiated distribution functions generate essential corrections arising from the finite-temperature regulation as implied in eq.(2.8).In the presence of such δ-function limits, the Cauchy theorem does not need to be applied to find a closed-form result.A more convenient prescription involves a clever splitting of the original contour in two (see sec. 2) and the completion of the two line integral halves.This approach is faithful to the original conditions of the fermionic sum-integral and allows a simple substitution procedure to replace highly non-trivial complex contour computations.

B. One-loop master integral with a mass scale
In [26], two methods were introduced for approaching the computation of the fermion one-loop master integral involving a single imaginary scale iµ.The integration order, where one starts from the spatial integral, excels in directly treating the (potential) p ≃ µ divergence for general α > 1 through a generalization of Euler's beta function.Adding another (mass) scale to the propagator would merely shift the location of the (potential) novel divergence while other properties stay similar, so that in the T = 0 limit, the massive variant where I α (µ, m, T ) ≡ I 0 α (µ, m, T ).The massive integral has the same convergence properties as the massless integral, with convergence achieved in the triangle-shaped subregion of eq.(A.1).In the opposing case of m > µ, we on the other hand automatically avoid not only the p ≃ µ divergence, but also any µ dependence on the result as will be explicitly demonstrated below.
In this appendix, we explicitly compute the low-temperature limit of eq.(B.1) using the regulatory methods of [26] as well as hypergeometric algebra.Subsequently, we construct a generalization of the IBP relation of eq. ( 4.2) in the low-temperature limit, including the additional mass scale.

B.1. Direct computation
In analogy with the zero-mass case of our eq.( 4.3) (see also [26]), computing first the d-dimensional dimensionally regularized spatial integral in eq.(B.1) yields wherein we simultaneously scaled out the chemical potential and took the low-temperature limit by discarding the lower part of the contour as it is exponentially suppressed.The remaining integral is most straightforwardly evaluated using the Cauchy theorem, assuming the usual convergence region of eq.(A.1).To this end, we apply eq. ( 2.12) and write First inspecting the case m > µ, we note that the last term above is purely imaginary.Inserting the real part of the remaining expression into eq.(B.2) yields then which agrees with the massive vacuum integral I α (m) from eq. (3.5), in (d + 1) dimensions, with no chemical potential dependence present.For m < µ, we proceed somewhat differently by first simplifying the last term in eq.(B.3) according to i 2 where on the last line we applied the so-called Pfaff transformation of [43]. 18The motivation for this transformation stems from the series expansion of the hypergeometric function 2 F 1 [a, b; c; z] only converging when its last argument satisfies z < 1.To obtain an expression satisfying this constraint, we further apply the transformation [43][44][45] 2 which allows us to write i 2 In the second iteration, we regulated the complex-valued power functions using Π + ≡ exp iπ 2 (1 + κ) as defined in eq.(4.5).
Taking now the real part of the above expression and the κ-regulator to vanish, we can write for m < µ

Re
i 2 where we extracted the real part via Re(Π + ) s ∼ cos(π(1 + κ)s/2) for s ∈ R. As a final step, we use Euler's reflection formula on the trigonometric function and combine this expression with eqs.(B.2)-(B.4),which upon further simplifying the result with eq. ( 4.1) leads to This expression allows studying the m ≪ µ regime in detail, and in particular, leads us to verify that the m → 0 limit indeed yields eq. ( 4.1) and that the result for I α (µ, m, T ) is continuous in the limit m → µ − .The latter follows from the result provided convergence is enforced by the parameter constraints in eq.(A.1).In combination with eq.(B.10), this shows the continuous limit The above result has a logical interpretation even outside of the parametric values, for which the limit m → µ remains formally convergent.Denoting χ = µ 2 − m 2 , we recognize that eq.(B.11) is the O(χ 0 ) term of a series expansion around χ = 0, with all other terms being of order O χ k or O χ d−2α 2 +k with k ∈ Z + .Whenever I α (m > µ) is UV-convergent, both subsets only contain positive powers, and similarly to our analytical continuation of the Euler gamma functions in I α (m > µ), we may regulate both the dimension and denominator exponent parameters to remove all terms except for the one shown in eq.(B.11).

B.2. Integration-by-parts with two scales
Integration-by-parts relations can also be derived for massive integrals as in eq.(B.1).Given the full solution in eq.(B.10), we expect to find inhomogeneous equations for master integral structures for m < µ.For eq. (B.1), the IBP approach aligns with the discussion of sec.4. By acting on the massive master integral with diagonal spatial and temporal total derivatives, we find where in the last line we combined spatial and temporal total derivatives into a d + 1 dimensional total derivative bilinear akin to identity (ibp.1c).After setting s = 0, we again discern two hierarchy regimes: Hierarchy m > µ.In the low-temperature limit, the right-hand side vanishes such that where we have used the results of sec.4.This naturally implies the familiar recursion from zerotemperature scalar field theory, i.e.
Hierarchy m < µ.In the opposite regime, we regulate the power-law term and the imaginary unit as above to find where we again used the Π + convention of eq.(4.5).This implies the validity of the identity

C. One-loop consistency checks and additional examples
This appendix complements the results of sec. 4 by providing consistency checks related to known standard properties of dimensionally regularized loop integrals.In particular, we focus on the oneloop result given in eq. ( 4.1) as well as its finite-temperature generalization Along the way, we also discuss the role of chemical potential in parametric differentiation and its relation to symmetries at the one-loop order.We further emphasize that our results should be formally considered analytical continuations of mathematically convergent integrals, which in the lowtemperature regime again implies parametric constraints akin to eq. (A.1).

C.1. Linear algebra and loop-momentum differentiation
For the sake of completeness, let us first demonstrate that the fermionic one-loop master integral allows the application of various linear decompositions while keeping the result intact.For this purpose, we multiply eq.(C.1) by unity 1 = (p 2 0 + p 2 )/[p 2 0 + p 2 ] and first evaluate the spatial integral, producing After moving to numerator structures containing exclusively spatial structures, all required symmetries arise from the dimensionally regularized beta function.
In sec.4, we also discussed the effect of the diagonal total derivative operator ∂ ∂p 0 • p 0 , where one simplification was found in the vanishing of total spatial derivatives, or conversely the generated boundary terms.This is a central element in dimensionally regularized integral algebra, and here we apply the one-loop master integral to confirm this result explicitly.By first performing the spatial integration, we retain all differential and linear algebra within the innermost integral.The diagonal operator can be re-written such that Upon performing the innermost integral, this operator structure introduces two terms: one proportional to the master integral and another with increased radial power, of which the latter can be geometrically interpreted as a master integral in a higher dimension.To this end, the total spatial derivative acting on the integrand of eq.(C.1) results in 19 We emphasize that a naive application of the residue theorem at T = 0 (cf.[26]) would not agree with this relation. 20 The results listed above and in sec.4 give great insight into expressions with trivial or quadratic numerator structures.To address the subset of cases with numerator structures linear in the zerocomponents of momenta, let us next discuss the case of the temporal derivative ∂ ∂p 0 , which acts on the distribution function.Unlike the bilinear operators introduced in sec.3, such an operator allows for relating primed (n ′ F ) integrals to ones with linear numerator structures that cannot be further simplified using linear algebra from other total derivatives.The simplest such structure is the temporal integro-differential relation in eq.(ibp.1b) for a fixed numerator parameter s = −1 which can again be confirmed to vanish by direct computation of both expressions on the right-hand side.The first integral on the right-hand side is given by eq.(A.6), which was evaluated using the methods introduced in sec.III of [26].By combining the Cauchy theorem and starting with the spatial integral, we obtain 19 To simplify the computation, we vary the spatial dimensionality of the master integral, Iα(d, µ, T ), and write it explicitly in the superscript. 20A naive application of the residue theorem to eq. (4.1) results in d+1−2α , which leads to a non-vanishing right-hand side of eq.(C.4).
The second expression in eq.(C.5) can be computed using methods introduced in secs. 2 and 4, yielding which indeed agrees with eq.(C.6), causing eq.(C.5) to vanish as expected.This is another consistency check between the two contour methods introduced in sec. 2 for the evaluation of the loop integrals in eqs.(2.12) and (2.14).

C.2. Parametric differentiation
Parametric differentiation or Feynman's trick can also be used to associate most single-scale (µ) oneloop integrals to the master integral I α (µ).The practical value of such an operation relies heavily on the assumption that the differentiated structure can be dealt with more easily than the nondifferentiated one, and that there are no scale-invariant contributions removed by the differentiation.
As for the latter concern, we note that the integrals discussed in this work contain two independent scales, µ and T , the latter of which is taken to an unessential limit (and hence should not contribute in a meaningful way to the results we are seeking).To this end, the µ-differentiation can only miss fully scale-independent, O(µ 0 ), contributions to the integral in question. 21 In the low-temperature limit, the distribution functions bring no scale to the energy dimension of the full integral while their derivatives lower it by one, so that all loop integrals of interest take the form of a power of µ.The η-regulators are on the other hand suppressed in our explicit results, whose low-temperature limits take forms such as lim Accordingly, the presence of dimensional regularization should be sufficient to avoid any loss of information in parametric differentiation.At the one-loop level, Feynman's trick can be used both as another consistency check as well as in the derivation of µ-dependent symmetries extending beyond the low-temperature limit.While the full solution is explicitly temperature-dependent, the only parameter we can differentiate is the chemical potential, which leads to non-unique identities involving the null space of ∂ µ , such as (C.10) 21 The situation is different for integrals with, e.g., both mass and chemical potential scales.In these cases, each parametric differentiation carries the risk of canceling out a part of the full solution.Such a case can be found in appendix B for the scale hierarchy m > µ.
For multi-loop computations and when applied directly, the trick generates differential equations to replace the inhomogeneous recursive equations, corresponding to the power-law behavior of the result in µ. 22 While computations involving parametric differentiation are not central to our work, we add a few examples below to demonstrate their connection to the formulae used here and in [26].
Similarly to many previous computations presented in this work, we need to first keep the temperature nonzero to ensure that we do not miss important contributions to our integrals that might vanish upon differentiation at exactly T = 0. Without any loss of information, we can then write the partial differentiation of the distribution function as (C.11) and subsequently, move the differentiation outside the p 0 -integral due to the closed nature of the integration contour.By differentiating the low-temperature limit corresponding to the rest of the integrand on the right-hand side, we find using methods discussed above (C.12) As stated above, this implies that the parametric differentiation encompasses all terms containing explicit µ-dependence.Furthermore, we note that the order of the zero-temperature limit and parametric differentiation does not affect the result.Using this flexibility with the limits, we can also confirm the result derived for master integrals equipped with linear p 0 -structures: which indeed agrees with eq.(C.7).Finally, we note that we can rather easily derive additional formulae for similar master integrals, including e.g.lim and that we can discern from eqs. (4.1) and (4.2) -or eq.(C.13) -that the p 0 -numerator in eq.(C.12) can be replaced by µ such that Here, we applied eq.(ibp.1c) on the last line to re-write the relation in terms of the operator P µ • ∂ ∂Pµ which is reversed to operators of the set in eq.(3.8).This relation further demonstrates the non-trivial nature of the p 0 = ±iµ substitutions arising from eq. (2.14) in the low-temperature limit.The relevant details of the delta function limits and contours have been discussed already in secs. 2 and 4.

D. Additional details of two-loop IBP computations
In this appendix, we present several extensions and consistency checks for the computations presented in sec 5.In the first part, we introduce conventions for the spatial integrations taking place in eq. ( 5.4) and similar relations, while in the second subsection, we study in more detail the vanishing integrals that contain propagators with shifted bosonic momenta.The third subsection then evaluates the group of two-loop integrals {S 11s | s ∈ N} at T = 0, using the cutting rule toolkit from [31] and providing a comparison point for our independent evaluation of such diagrams.The last subsection finally utilizes spatial IBP operators from [37] to derive IBP relations for the integrals S 111 and S 112 that explicitly differ from what one would expect at µ = 0.These results act as further cross-checks of our full (d + 1)-dimensional IBP relations (5.13) and (5.22).

D.1. Vacuum integrals with complex mass scales
Let us begin by considering the one-loop integral of eq.(4.6) and the two-loop integrals of eq.(5.4) without explicit replacements.Akin to sec. 4 of [26], we may regulate the quadratic scales using iµ → Π + µ from eq. (4.5), which allows using the generalized Euler beta functions in the dimensionally regularized radial integrations.At the one-and two-loop levels, the initial spatial momentum integration corresponds to vacuum integrals of the type where in the first integral one assumes the hierarchy m 2 1 ≥ m 2 2 and in the latter two cases the collinear relation m 3 = m 1 + m 2 that gives rise to the factorization of the vacuum sunset reported in [32,46].
In the evaluation of the maximally primed type B integrals of sec.5.1, we may apply the above collinear results because one of the three complex mass scales in the corresponding d-dimensional integrals always vanishes due to the substitutions of p 0 and q 0 in relations such as eq.(5.4).The specific mapping between the masses m i appearing above and the momenta p 0 , q 0 , and p 0 ± q 0 is determined by the integral in question, typically leading to only one of the three terms in eqs.(D.2) and (D.3) obtaining a nonzero value.To remove the vanishing mass scale from the denominator, we typically use the one-loop massive vacuum IBP relation I α+1 (m) = −(d − 2α)/(2αm 2 )I α (m) (see appendix E), which leads to two cases p 0 I 2 (p 0 ) q 0 I 2 (q 0 ) , for p 0 + q 0 = 0 .(D.5) To finally arrive at results such as eq.( 5.5), we must still connect the temporal momentum components to the substitutions µ(Π + ) ±1 .In this way, eq.(D.4) is associated with the p 0 = q 0 = µΠ + substitution in the first term on the right-hand side of eq.(5.4), while eq.(D.5) is associated with the second term there.
One more subtlety, however, remains: for S 111 (µΠ + , µ/Π + , O(κ)), a straightforward generalization of the hypergeometric expression (D.1) would provide a complex-valued result for an integral that is by construction its own complex conjugate.To remedy this, we must return to the defining integral (obtained after shifting p → p + q) where we can use standard Feynman parametrization if α 1 , α 2 > 0. By regulating again the complexvalued mass scales such that i → Π + , we may straightforwardly evaluate the two spatial loop integrations, obtaining a hypergeometric integral of the form where again α {i} = j={i} α j .
To proceed from here, let us next assume that α 3 ∈ N. Since the denominator and numerator powers differ by (α 23 − d 2 ) + (α 13 − d 2 ) − (α 123 − d) = α 3 for the remaining hypergeometric integral, we can rewrite its integrand by multiplying with unity where the subtraction p 2 0 − q 2 0 = O(κ) is proportional to the κ-regulator scale.Writing the numerator of the multiplier in terms of a binomial expansion, we may then remove the contribution proportional to the regulator, obtaining where we have reversed the Feynman parametrization of eq.(A.9) to find the factorization explicitly manifested.
The result we have arrived at is equivalent to performing the binomial expansion symmetrically in terms of both p 0 and q 0 , but cannot be straightforwardly extended to α 3 ∈ R + .This can be seen by carefully studying Newton's extension to the binomial theorem in the context of eq.(D.9).Indeed, for α 3 ∈ (0, 1), p 0 = µΠ + and q 0 = µ(Π + ) −1 , we obtain where we have paid close attention to the signs of the O(κ) terms, according to which the integral has been split in two.In the end, we can collect the leading contributions from both integration intervals such that (D.13) = µ 2α 3 cos(πα 3 ) which retains the familiar factorizing structure of eq.(D.10) when α 3 ∈ N.
A natural generalization to eq. (D.10) is found by adding a mass m ∈ R + to the propagators.Such a scale can be inserted to both fermionic propagators of eq.(D.6) (associated to p 0 and q 0 ) while still conserving the binomial expansion condition (p 2 0 + m 2 ) − (q 2 0 + m 2 ) = O(κ).The regulator κ can be fitted such that for µ > m we may replace Π + µ → Π + µ 2 − m 2 in the thermal integral akin to the one-loop example of appendix B. For m > µ, the result on the other hand trivializes to the classical real-valued sunset with only one scale p 2 0 = q 2 0 = m 2 − µ 2 > 0. We emphasize, though, that the m > µ hierarchy leads to the vanishing of the type B term due to subtractions akin to eq. (5.3), making the corresponding IBP relation independent of µ.
Next, we apply eq.(D.10) to an explicit evaluation of the sunset S 111 with both Fermi-Dirac distributions differentiated.In practice, we write the integrand in terms of the spatial integral S 111 (p 0 , q 0 , 0) and use the contour deformations from eq. (2.14) to recover the solution to eq. (5.5), where standard trigonometric formulae and Euler's reflection formula have been applied.Similarly, it is straightforward to extend the result to the full class of maximally differentiated two-loop vacuum bubbles, for which α 1 , α 2 ∈ R + and α 3 ∈ N. Applying again eq.(D.10), we obtain Finally, we note that numerators with integer powers of temporal momentum components {p 0 } can be implemented through eq. ( 5.4) and hence do not introduce any further complications.Following the conventions of sec.5, the two type B terms required in the derivation of S 111 read (D.17)

D.2. Momentum-shifted bosonic propagators
Next, we take a closer look at the zero-temperature limits of the integrals S 102 and S 012 , appearing in eq.(5.7).Following the convention of eq. ( 2.3), we may write where we have used the sum-integral definitions of eq. ( 2.3) and noted the bosonic nature of the momentum R ≡ Q − P .To proceed to the low-temperature limit, we first perform the spatial integral and then compute the remaining zero-component integral through the step function limit of n B in eq.(2.6) (cf.sec.III of [26]): The above expression formally corresponds to the fermionic one-loop integral of eq.(4.1).Its dependence on the powers of the regulator scale η (from the contour) indicates the factorizing of the R-integral in eq.(D.18), leading to the entire diagram vanishing in dimensional regularization 23 .This result is retained also with a numerator linear in q 0 after splitting the computation into two upon writing q 0 = r 0 + p 0 .An integral with a p 0 -numerator also follows the same steps as above while an integral with an r 0 -numerator involves an imaginary part taken during the remaining r 0 -integral, similar to eq. (C.6).
Next, we inspect the closely related integral of eq.(5.12), which contains one primed and one unprimed distribution function, marking our first example of a simple type C integral.Using the formulae of sec.5, we obtain for it where the right-hand side follows after performing the d-dimensional spatial integrals.Focusing next on the remaining two integrals, we change the integration variable of the q 0 -integral via q 0 → r 0 = q 0 − p 0 and also split the p 0 -integral using the line integral description.This leads to where η p is the regulator on the p-integral contour and the ± symbols correspond to the ±iη p shifts of the bosonic contour along the imaginary axis.In the T → 0 limit, the leading behaviour of distribution functions is solely determined by the imaginary part of loop momenta [cf.eq.(2.8)], ñ′ which we can be immediately applied given that all the momentum dependence outside of distribution functions is of the form of generalized monomials 24 and where the p 0 line integrals fix the value of Im(p 0 ).Thus, the r 0 integration akin to eq. (2.14) becomes (D.23) 23 Taking the limit η → 0 before letting the dimensional regulator vanish results in a vanishing scaleless integral as expected for T = 0 bosonic integrals. 24Monomials do not introduce problematic poles and allow a straightforward taking of the naive T → 0 limit in both differentiated and non-differentiated distributions.
which shows the proprtionality to the regulator in the form ∝ η d−2α 2 p .Similarly, each term in eq.(D.21), through the evaluation of the r 0 integrals, contains a factorizing dependence on the regulator powers.As in the previous example of eq.(D. 19), this indicates that the full integral vanishes in dimensional regularization -a property that is trivially retained for unit-valued numerators.Should we take the additional regulator η p to zero beforehand, the integral containing δ-sequence expressions in the T → 0 limit would on the other hand be identical to one with a bosonic δ-sequence (as generated by differentiation).

D.3. Evaluation of the sunset integral with cutting rules
For the sake of completeness, we present next an alternative evaluation of the standard sunset integral at T = 0 using the so-called cutting rules of [31].By setting α 1 = α 2 = 1 and α 3 = s ∈ N in eq. ( 5.1), we are indeed able to take the zero-temperature limit and compute the remaining integral using the residue theorem [26,31], evaluating multiple integrals at once.We emphasize, however, that raising either of the other two exponent parameters α 1 , α 2 to values higher than unity would require the introduction of corrections akin to the differentiated distribution functions discussed in sec.2.
In the application of the cutting rules, we write the residues as cuts along the fermionic propagators, obtaining thereby the expression where we again use capital letters to signify Euclidean (d + 1)-momenta integrated over R d+1 using dimensional regularization.The first of the three terms on the right-hand side of this expression is said to correspond to a 0-cut, the second to a 1-cut, and the third to a 2-cut contribution to the Feynman integral.Given that all scaleless integrals vanish in dimensional regularization, the 0-cut contribution on the second line above clearly vanishes, and the same can be seen to be true for the 1-cut term given that the substitution enforces P 2 = 0.For a general exponent, s, the remaining integral and also the full result therefore becomes where the angular parameter in the dot product reads z = p • q/(pq).To evaluate the z-integral, we next change the integration variable as z → y = (1 + z)/2 to obtain a scaled Euler's beta function where we applied the one-loop result of eq.(4.1).It is easy to verify that the full result indeed reproduces our earlier expressions for S 111 in eq.(5.13) and for S 112 in eq.(5.22).

D.4. Sunset integral and spatial IBP relations
In sec.5.2, we derived factorizing low-temperature results for the fermionic sunset integrals S 111 and S 112 in eqs.(5.13) and (5.22), respectively.At nonzero temperature but vanishing chemical potentials, the sunset S 111 is on the other hand known to vanish identically, as demonstrated in [24,37] using spatial IBP relations and initially argued up to O(ǫ) in [47].Given the non-vanishing nature of our corresponding T = 0, µ = 0 results, it is interesting to compare these calculations and attempt to pinpoint the reason for the observed difference.
To gain further insights into the sunset integral, we apply the diagonal spatial IBP identities from the set (3.8) using an in-house Laporta-type reduction [25], implemented in FORM [36].The relations are homogeneous and can be written compactly in operator form (again omitting the two loop master integral S s 1 s 2 α 1 α 2 α 3 on which each operator acts) as where the two remaining bilinear combinations can be found via the substitution 1 ↔ 2 corresponding to p i ↔ q i .These relations can be combined to yield operator ansätze leading to an explicit one-loop factorization of two-loop expressions with collinear temporal scales.The ansatz we choose reads (cf.[37]) where we dubbed bosonic integrals with the superscript b and utilized the master integrals of eq.(A.6) and where the fully bosonic integrals I b α (T ) vanish here in the T → 0 limit.By inspecting eqs.(D.29) and (D.30), we immediately recover the high-temperature fermionic result of [37] using the fact that sum-integrals of the form

E. Two-loop factorization at finite density
In sec.5, we encountered the partially unexpected factorization of two-loop thermal vacuum integrals into a sum of products of one-loop integrals, the immediate consequences of which we will now investigate further.Such a property has been known to hold for d-dimensional T = µ = 0 sunset integrals with a collinear mass signature [32,46], for collinear masses m 3 = m 1 + m 2 . (E.1) The coefficients, c ij , and exponents {f ij , g ij } have recently been determined in a closed form for integer-valued exponents α n ∈ N in [32] and agree with values previously suggested by recursion formulae such as two-loop vacuum IBP reduction.In our present context, it is a nontrivial and very 25 The quadratic dependence on the dimension arises from contractions of the form interesting question, to which extent such a factorization extends to the description of the thermal sunset S α 1 α 2 α 3 (µ, T ) in terms of products of the one-loop integrals I α (µ, T ).In the following, we present such a generalization using the formalism of secs. 2 and 3.
Finally, we note that due to the factorization (via collinear scales) of the d-dimensional vacuum sunset S α 1 α 2 α 3 for α 1 , ...α 3 ∈ N [32], it is in fact all thermal two-loop bubbles S α 1 α 2 α 3 and not just the above special cases that factorize, and this appears to happen for arbitrary values of T and µ.Specific examples of this generic result include the low-temperature regime of T → 0 but µ = 0, where the factorization emerges explicitly as observed in sec.5.1, as well as the high-temperature regime T = 0 but µ = 0, where the factorization was demonstrated already in [24,48] for a few special cases.The factorization becomes particularly transparent in the IBP approach, where it is automatically built in the relations (ibp.2a)-(ibp.2d);see sec.5.2 and the appendix D.1.We emphasize, however, that it is not a universal property of multi-loop thermal integrals nor do we expect it to generalize to higher loop orders or massive propagators.Even in these cases, the maximally primed type B integrals and the corresponding δ-functions nevertheless remain tractable at low temperatures, thus simplifying the corresponding computations.
One example of possible extensions of the above results is the two-loop fermionic sunset with equal positive mass scales m 2 and chemical potential µ in the fermionic propagators, Evaluating this integral directly using factorization is not possible, since the d-dimensional integral is no longer collinear.Its maximally primed expression, however, still has this property and can thus be evaluated directly using Feynman parametrization and the corresponding hypergeometric integral (D.10) or using the factorization observed in [32].The result can be given readily for µ > m in the form D p D q S 111 (µ, m, T ) while for m > µ it is observed to vanish; see appendices B and D.1 for more context and details.Such differentiated integrals are again seen to appear as a part of (d + 1)-dimensional IBP reduction akin to eq. (5.16) or (5.18), thus contributing to the determination of the full massive thermal integral.