aainstitutetext: Centro de Ciencias Exactas, Departamento de Ciencias Básicas,
Universidad del Bío-Bío, Avenida Andrés Bello 720, Chillán,
Chilebbinstitutetext: Grenfell Campus, Memorial University of Newfoundland, 20
University Dr., Corner Brook, NL, A2H 5G4, Canadaccinstitutetext: Department
of Physics and Helsinki Institute of Physics, P.O. Box 64, University of
Helsinki, Finland

# Factorizing two-loop vacuum sum-integrals

Andrei I. Davydychev a,c Pablo Navarrete a and York Schröder
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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.

###### Keywords:

Higher-Order Perturbative Calculations, Effective Field Theories of QCD,
Renormalization and Regularization, Thermal Field Theory

##  1 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) Gross:2022hyw  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 ${\mu_{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 Kapusta:2023eix ;
Bellac:2011kqa ; Laine:2016hma ; Ghiglieri:2020dpq  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 Arnold:1994ps ; Arnold:1994eb . 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 Boyd:1996bx ; Borsanyi:2012ve ; Giusti:2016iqr  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 non-abelian setting of QCD, where the gauge symmetry
severely constrains (but also guides) methods that go beyond naive
perturbative expansions111Another approach to account for low-energy screening
effects in a gauge theory is hard-thermal-loop perturbation theory
Braaten:1989kk ; Braaten:1989mz ; Andersen:1999fw , which has been developed
up to and including three-loop contributions Haque:2014rua ; Andersen:2015eoa
, superseding simpler resummation schemes that had been developed for scalar
theories, such as screened perturbation theory Karsch:1997gj .. This has led
to the development and application of dimensionally reduced EFT’s
Ginsparg:1980ef ; Appelquist:1981vg ; Kajantie:1995dw ; Braaten:1995jr  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 Kajantie:2002wa ; Vuorinen:2003fs ;
Gorda:2023mkk , 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
DiRenzo:2006nh , 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 Laine:2006cp ; Schroder:2006vz . 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 Navarrete:2022adz  for the state-of-the-art at
${\mu_{f}}=0$, and Gorda:2023mkk  for ${\mu_{f}}\neq 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 sum-integrals 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 ${\mu_{f}}$ are present a recursive
treatment might still apply to reduce the generic two-loop case to one-loop
factors Osterman:2023tnt .

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 sum-integral 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.

##  2 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 $\tau$ (which emerges after Wick rotating
$it\rightarrow\tau$, leading to $e^{iS}\rightarrow e^{-S_{\rm E}}$ in the path
integral) that has period $1/T$ set by the inverse temperature. Dimensionally
regularized integrals then employ $d=3-2\varepsilon$ in our conventions. Since
we do not $\varepsilon$-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 $\tau$.
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 bosonic222For fermions one would have
$P_{0}=(2n_{p}+1)\pi T+i{\mu_{f}}$ where ${\mu_{f}}$ is the chemical potential
of the respective fermion flavor. In this work we focus on the bosonic case,
however, and only briefly comment on fermionic generalizations in the
conclusions. Matsubara frequencies $P_{0}=2n_{p}\pi T$ where
$n_{p}\in\mathbb{Z}$. Massless propagators are
$1/[P^{2}]=1/[(P_{0})^{2}+p^{2}]$, and we employ the sum-integral symbol
defined as

$\displaystyle{\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int\,}_{\\!\\!\\!\\!\\!\raise-0.90417pt\hbox{$\scriptstyle{P}$}}\equiv
T\sum_{n_{p}\in\mathbb{Z}}\int\frac{{\rm d}^{d}p}{(2\pi)^{d}}\;.$ (1)

While some of the machinery used in organizing and reducing perturbative
expansions, such as e.g. integration-by-parts (IBP) methods Chetyrkin:1981qh ;
Laporta:2000dsw , carries over from zero to finite temperatures
Nishimura:2012ee , 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
$\varepsilon$ expansion, typically relying on tedious subtractions of sub-
divergences in one-loop substructures Arnold:1994ps , 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
($\eta\in\mathbb{N}_{0}$)

$\displaystyle I_{\nu}^{\eta}(d,T)$ $\displaystyle\equiv$
$\displaystyle{\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int\,}_{\\!\\!\\!\\!\\!\raise-0.90417pt\hbox{$\scriptstyle{P}$}}\frac{(P_{0})^{\eta}}{[P^{2}]^{\nu}}\;.$
(2)

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\geq 0$)

$\displaystyle J_{\nu}(d,m)$ $\displaystyle\equiv$
$\displaystyle\int\frac{{\rm
d}^{d}p}{(2\pi)^{d}}\frac{1}{[m^{2}+p^{2}]^{\nu}}\;=\;\big{[}m^{2}\big{]}^{\frac{d}{2}-\nu}\frac{\Gamma(\nu-\tfrac{d}{2})}{(4\pi)^{d/2}\,\Gamma(\nu)}\;.$
(3)

Taking advantage of the simple mass dependence of the continuum integral $J$,
this gives

$\displaystyle I_{\nu}^{\eta}(d,T)$ $\displaystyle=$ $\displaystyle
T\sum_{n\in\mathbb{Z}}(2n\pi T)^{\eta}\,J_{\nu}(d,|2n\pi T|)$ (4)
$\displaystyle=$ $\displaystyle
T\delta_{\eta,0}J_{\nu}(0)+\big{[}1+(-1)^{\eta}\big{]}T\sum_{n=1}^{\infty}(2n\pi
T)^{\eta}J_{\nu}(d,2n\pi T)$ $\displaystyle=$ $\displaystyle
0+\big{[}1+(-1)^{\eta}\big{]}T\sum_{n=1}^{\infty}(2n\pi
T)^{\eta+d-2\nu}\,J_{\nu}(d,1)$ $\displaystyle=$
$\displaystyle\frac{[1+(-1)^{\eta}]\,T\,\zeta(2\nu-\eta-d)}{(2\pi
T)^{2\nu-\eta-d}}\,\frac{\Gamma(\nu-\tfrac{d}{2})}{(4\pi)^{d/2}\,\Gamma(\nu)}\;,$

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. (3) to scale out the mass to finally identify the single Matsubara
sum as a zeta function. As it should be, due to the sum-integral’s symmetry
under a momentum shift $P\rightarrow-P$, all $I_{\nu}^{\eta_{\rm odd}}(d,T)$
vanish identically, as is reflected by the prefactor $[1+(-1)^{\eta}]$.
Therefore, we only encounter zeta values of the form $\zeta(n_{\rm even}-d)$
at one loop333At finite ${\mu_{f}}$, the symmetry argument breaks down, such
that the corresponding one-loop tadpoles do not vanish for odd values of
$\eta$. Furthermore, at finite $T$ and ${\mu_{f}}$ one encounters Hurwitz
(instead of Riemann) zeta functions., 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.

##  3 Massless bosonic two-loop vacuum sum-integrals: Setup

The general scalar massless bosonic 2-loop vacuum sum-integral is defined as

$\displaystyle L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$
$\displaystyle\equiv$
$\displaystyle{\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int\,}_{\\!\\!\\!\\!\\!\raise-0.90417pt\hbox{$\scriptstyle{PQ}$}}\frac{(P_{0})^{\eta_{1}}\,(Q_{0})^{\eta_{2}}\,(P_{0}-Q_{0})^{\eta_{3}}}{[P^{2}]^{\nu_{1}}\,[Q^{2}]^{\nu_{2}}\,[(P-Q)^{2}]^{\nu_{3}}}\;.$
(5)

Imitating the strategy presented in the evaluation of eq. (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

$\displaystyle L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$
$\displaystyle=$ $\displaystyle\frac{T^{2}}{(2\pi
T)^{2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-2d}}\sum_{n_{1},n_{2}\in\mathbb{Z}}n_{1}^{\eta_{1}}\,n_{2}^{\eta_{2}}\,(n_{1}-n_{2})^{\eta_{3}}\,B^{\nu_{1},\nu_{2},\nu_{3}}_{n_{1},n_{2},n_{1}-n_{2}}(d)\;,$
(6) $\displaystyle B^{\nu_{1},\nu_{2},\nu_{3}}_{m_{1},m_{2},m_{3}}(d)$
$\displaystyle\equiv$ $\displaystyle\int\\!\\!\frac{{\rm
d}^{d}p}{(2\pi)^{d}}\int\\!\\!\frac{{\rm
d}^{d}q}{(2\pi)^{d}}\frac{1}{[m_{1}^{2}+p^{2}]^{\nu_{1}}\,[m_{2}^{2}+q^{2}]^{\nu_{2}}\,[m_{3}^{2}+(p-q)^{2}]^{\nu_{3}}}\;,$
(7)

where ${\Sigma\nu_{i}}\equiv\nu_{1}+\nu_{2}+\nu_{3}$ and
${\Sigma\eta_{i}}\equiv\eta_{1}+\eta_{2}+\eta_{3}$. 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 Davydychev:1992mt . 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
Davydychev:2022dcw  by a detailed analysis of the corresponding IBP relations,
leading to a general factorization formula of the type $B\rightarrow\sum
J\cdot J$. The idea is now to use results of Davydychev:2022dcw  for the
massive 3d vacuum tadpole $B$ of eq. (6), and do the double-sum afterwards.
Amazingly, this will allow us to present a proof of factorization for,
schematically, $L\rightarrow\sum I\cdot I$, where the $I$ are the 1-loop sum-
integrals of eq. (4).

Further following the successful 1-loop strategy as explained below eq. (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 obtain444Using the shortcut
$\delta_{j}\equiv\delta_{j,0}$ for the Kronecker delta function.

$\displaystyle L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$
$\displaystyle=$ $\displaystyle\frac{T^{2}\,[1+(-1)^{\Sigma\eta_{i}}]}{(2\pi
T)^{2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-2d}}\,\Big{\\{}\tfrac{1}{2}\,\delta_{\eta_{1}}\delta_{\eta_{2}}\delta_{\eta_{3}}\,{\color[rgb]{0,0,0}B_{0,0,0}^{\nu_{1},\nu_{2},\nu_{3}}}+\zeta(2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-2d)\\!\times$
(8)
$\displaystyle{}\times\\!\Big{[}(-1)^{\eta_{3}}\delta_{\eta_{1}}{\color[rgb]{0,0,0}B_{0,1,1}^{\nu_{1},\nu_{2},\nu_{3}}}+\delta_{\eta_{2}}{\color[rgb]{0,0,0}B_{0,1,1}^{\nu_{2},\nu_{1},\nu_{3}}}+\delta_{\eta_{3}}{\color[rgb]{0,0,0}B_{0,1,1}^{\nu_{3},\nu_{1},\nu_{2}}}+2^{\eta_{3}}(-1)^{\eta_{2}}{\color[rgb]{0,0,0}B_{1,1,2}^{\nu_{1},\nu_{2},\nu_{3}}}\Big{]}$
$\displaystyle{}+{\color[rgb]{0,0,0}{\bar{H}}_{\eta_{3},\eta_{2},\eta_{1}}^{\nu_{3},\nu_{2},\nu_{1}}}+(-1)^{\eta_{3}}{\color[rgb]{0,0,0}{\bar{H}}_{\eta_{3},\eta_{1},\eta_{2}}^{\nu_{3},\nu_{1},\nu_{2}}}+(-1)^{\eta_{2}}{\color[rgb]{0,0,0}{H}_{\eta_{2},\eta_{1},\eta_{3}}^{\nu_{2},\nu_{1},\nu_{3}}}+(-1)^{\eta_{2}}{\color[rgb]{0,0,0}{H}_{\eta_{1},\eta_{2},\eta_{3}}^{\nu_{1},\nu_{2},\nu_{3}}}\Big{\\}}\;,\qquad$
with
$\displaystyle{H}_{\eta_{a},\eta_{b},\eta_{c}}^{\nu_{a},\nu_{b},\nu_{c}}(d)\equiv\sum_{n_{1}>n_{2}>0}n_{1}^{\eta_{a}}n_{2}^{\eta_{b}}(n_{1}+n_{2})^{\eta_{c}}\,B_{n_{1},n_{2},n_{1}+n_{2}}^{\nu_{a},\nu_{b},\nu_{c}}(d)\;,$
(9) and
$\displaystyle{\bar{H}}_{\eta_{a},\eta_{b},\eta_{c}}^{\nu_{a},\nu_{b},\nu_{c}}(d)\equiv\sum_{n_{1}>n_{2}>0}(n_{1}-n_{2})^{\eta_{a}}n_{2}^{\eta_{b}}n_{1}^{\eta_{c}}\,B_{n_{1}-n_{2},n_{2},n_{1}}^{\nu_{a},\nu_{b},\nu_{c}}(d)\;.$
(10)

Naturally, the sum-integral $L$ vanishes for odd ${\Sigma\eta_{i}}$, as is
immediately clear from the definition eq. (5), considering the integrand’s
symmetry under the simultaneous shifts $P\rightarrow-P$ and $Q\rightarrow-Q$.
The special-mass cases $B^{\nu_{1},\nu_{2},\nu_{3}}_{0,0,0}(d)$ and
$B^{\nu_{1},\nu_{2},\nu_{3}}_{0,1,1}(d)$ are known analytically in $d$
dimensions Vladimirov:1979zm ,

$\displaystyle B_{0,0,0}^{\nu_{1},\nu_{2},\nu_{3}}(d)$ $\displaystyle=$
$\displaystyle 0\mbox{~{}(massless tadpoles vanish in dimensional
regularization)}\;,$ (11) $\displaystyle
B_{0,1,1}^{\nu_{1},\nu_{2},\nu_{3}}(d)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(\nu_{1}+\nu_{2}-\frac{d}{2})\Gamma(\nu_{1}+\nu_{3}-\frac{d}{2})\Gamma(\frac{d}{2}-\nu_{1})\Gamma({\Sigma\nu_{i}}-d)}{\Gamma(\nu_{2})\Gamma(\nu_{3})\Gamma^{2}(1-\frac{d}{2})\Gamma(\frac{d}{2})\Gamma(\nu_{1}+{\Sigma\nu_{i}}-d)}\,B_{1,1,0}^{1,1,0}(d)$
(12) with $\displaystyle
B_{1,1,0}^{1,1,0}(d)\equiv\frac{\Gamma^{2}(1-\frac{d}{2})}{(4\pi)^{d}}=\big{[}J_{1}(d,1)\big{]}^{2}\quad\mbox{(see
eq.~{}(\ref{eq:1loopJ}))}\;,$ (13)

while $B_{1,1,2}^{\nu_{1},\nu_{2},\nu_{3}}(d)$ can be extracted from
Davydychev:1992mt  (cf. appendix C of Davydychev:2022dcw ). Note however that
these specific results are not necessarily relevant for our treatment of $L$,
as can be appreciated from our final result eq. (34) below.

$n_{1}$$n_{2}$$=$$n_{1}$$n_{2}$$+$$n_{1}$$n_{2}$ Figure 1:  Different regions
in the Matsubara double-sum, as used for deriving eq. (8). 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 $\sum_{n_{1}>n_{2}>0}$.

##  4 Massive two-loop vacuum integrals: Factorization formula

Let us now examine the massive two-loop scalar vacuum integrals
$B^{\nu_{1},\nu_{2},\nu_{3}}_{m_{1},m_{2},m_{3}}(d)$ that were defined in eq.
(7), 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 Davydychev:1992mt . 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. (8) we only need
$m_{3}=m_{1}+m_{2}$. For this special “collinear” case, the general results of
Davydychev:1992mt  can be simplified considerably.

Imposing the linear relation between the three masses, one can set up an
integration-by-parts (IBP) recursion that reduces the integral $B$ for all
positive values of $\nu_{i}\in\mathbb{N}$ towards one or zero. Relevant
explicit formulas can be found e.g. in eqs. (92,93,97) of Tarasov:1997kx . We
have recently been able to solve these IBP recursion relations explicitly
Davydychev:2022dcw , and obtained (from now on we implicitly assume
$m_{3}=m_{1}+m_{2}$ in all equations; also,
${\Sigma\nu_{i}}=\nu_{1}+\nu_{2}+\nu_{3}$)

$\displaystyle B^{\nu_{1},\nu_{2},\nu_{3}}_{m_{1},m_{2},m_{3}}(d)$
$\displaystyle=$ $\displaystyle
B_{1,1,0}^{1,1,0}(d)\,\Bigg{\\{}\sum_{j=1-\nu_{1}}^{\nu_{2}-1}(-1)^{\Sigma\nu_{i}}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\,m_{1}^{d-{\Sigma\nu_{i}}+j}\,m_{2}^{d-{\Sigma\nu_{i}}-j}$
(14)
$\displaystyle~{}~{}~{}~{}~{}+\sum_{j=1-\nu_{1}}^{\nu_{3}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\,m_{1}^{d-{\Sigma\nu_{i}}+j}\,m_{3}^{d-{\Sigma\nu_{i}}-j}$
$\displaystyle~{}~{}~{}~{}~{}+\sum_{j=1-\nu_{2}}^{\nu_{3}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\,m_{2}^{d-{\Sigma\nu_{i}}+j}\,m_{3}^{d-{\Sigma\nu_{i}}-j}\;\Bigg{\\}}\Bigg{|}_{\scriptsize
m_{3}=m_{1}+m_{2}}\;,$

where the coefficients $c^{({\Sigma\nu_{i}})}_{\nu_{a},\nu_{b};j}(d)$ are
rational functions in $d$ that obey a number of symmetries, such as

$\displaystyle c^{({\Sigma\nu_{i}})}_{\nu_{a},\nu_{b};j}(d)$ $\displaystyle=$
$\displaystyle c^{({\Sigma\nu_{i}})}_{\nu_{b},\nu_{a};-j}(d)\;,\quad
c^{({\Sigma\nu_{i}})}_{\nu_{a},\nu_{b};0}(d)=c^{({\Sigma\nu_{i}})}_{\nu_{b},\nu_{a};0}(d)\;,$
(15)

and are given explicitly by Davydychev:2022dcw

$\displaystyle c^{({\Sigma\nu_{i}})}_{\nu_{a},\nu_{b};j}(d)$ $\displaystyle=$
$\displaystyle\frac{(-1)^{{\Sigma\nu_{i}}-n_{j}+1}\,\big{(}1-\frac{d}{2}\big{)}_{n_{j}-j-1}}{2\,\big{(}\frac{1}{2}\big{)}_{n_{j}-\nu_{b}-\nu_{c}}\,\big{(}\frac{1}{2}\big{)}_{n_{j}-j-\nu_{a}-\nu_{c}}\,(\nu_{c}-1)!}\times$
$\displaystyle\times$ $\displaystyle\\!\\!\\!\\!\\!\\!\sum_{k={\rm
max}(1\\!+\\!j,1)}^{{\rm
min}(\nu_{a}\\!+\\!j,\nu_{b})}\frac{\big{(}\frac{d}{2}-n_{j}+1\big{)}_{k-1}\,(n_{j}-k-1)!}{\big{(}\frac{d+3}{2}\\!-\\!{\Sigma\nu_{i}}\big{)}_{n_{j}-k}\big{(}\frac{1}{2}\big{)}_{\nu_{c}-n_{j}+k}\,(k\\!-\\!1)!\,(k\\!-\\!j\\!-\\!1)!\,(\nu_{b}\\!-\\!k)!\,(\nu_{a}\\!-\\!k\\!+\\!j)!}$

with Pochhammer symbols $(a)_{\nu}\equiv\frac{\Gamma(a+\nu)}{\Gamma(a)}$ and
integers $n_{j}=\lceil\frac{{\Sigma\nu_{i}}+j}{2}\rceil$. Note that using eq.
(14) in eq. (8) 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. (14) for two special cases that appear in eq.
(8):

$\displaystyle B_{0,m,m}^{\nu_{1},\nu_{2},\nu_{3}}(d)$ $\displaystyle=$
$\displaystyle
m^{2d-2{\Sigma\nu_{i}}}\,B_{1,1,0}^{1,1,0}(d)\,\Big{\\{}\beta^{\nu_{1},\nu_{2},\nu_{3}}(d)\Big{\\}}\;,$
(17) with
$\displaystyle\beta^{\nu_{1},\nu_{2},\nu_{3}}(d)=\sum_{j=1-\nu_{2}}^{\nu_{3}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\;,$
(18) $\displaystyle B_{m,m,2m}^{\nu_{1},\nu_{2},\nu_{3}}(d)$ $\displaystyle=$
$\displaystyle
m^{2d-2{\Sigma\nu_{i}}}\,B_{1,1,0}^{1,1,0}(d)\,\Big{\\{}r_{1}^{\nu_{1},\nu_{2},\nu_{3}}(d)+2^{d-2}\,r_{2}^{\nu_{1},\nu_{2},\nu_{3}}(d)\Big{\\}}\;,$
(19) with $\displaystyle
r_{1}^{\nu_{1},\nu_{2},\nu_{3}}(d)=(-1)^{\Sigma\nu_{i}}\sum_{j=1-\nu_{1}}^{\nu_{2}-1}c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\;,$
(20) and $\displaystyle
r_{2}^{\nu_{1},\nu_{2},\nu_{3}}(d)=2^{2-{\Sigma\nu_{i}}}\Big{[}\sum_{j=1-\nu_{1}}^{\nu_{3}-1}(-\tfrac{1}{2})^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)+\sum_{j=1-\nu_{2}}^{\nu_{3}-1}(-\tfrac{1}{2})^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\Big{]}\,.\qquad$
(21)

##  5 Massless bosonic two-loop vacuum sum-integrals: Matsubara sums

The final task is now to evaluate the remaining double sums ${H}$ and
${\bar{H}}$ in eq. (8). Having eq. (14) at hand, the mass structure is fixed
explicitly, and we can work on performing the sums without actually specifying
the rational coefficient functions
$c^{({\Sigma\nu_{i}})}_{\nu_{a},\nu_{b};j}(d)$ yet.

We will now show how the specific linear combination of sums ${H}$ and
${\bar{H}}$ that occurs in eq. (8) combine to

  * (a)

cancel all $Z$ (for a definition, see eq. (51)) in the last two terms of eq.
(8),

  * (b)

cancel all $\zeta(i,j)$ (see eq. (46)) in the sum of all four terms of the
last line of eq. (8),

  * (c)

cancel all remaining single $\zeta(i)$ in eq. (8) to

  * (d)

leave us with the expected tadpole reductions that contain products
$\zeta(i)\,\zeta(j)$ containing only $\zeta(n_{\rm even}-d)$ allowing us to
map onto the 1-loop sum-integrals of eq. (4),

all of which allows us to formulate our final result in a compact way.

###  5.1 Proof of statement (a)

Plugging eq. (14) into eq. (9), 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
$n_{1}^{\eta_{1}}=\sum_{k=0}^{\eta_{1}}{\eta_{1}\choose
k}(-n_{2})^{k}(n_{1}+n_{2})^{\eta_{1}-k}$), shifting summation indices
appropriately (e.g. $k\rightarrow\eta_{3}-k$ and $j\rightarrow-j$) and using
the symmetry eq. (15) to collect similar expressions, the last two terms of
eq. (8) combine to

$\displaystyle\frac{{H}_{\eta_{2}\eta_{1}\eta_{3}}^{\nu_{2}\nu_{1}\nu_{3}}+{H}_{\eta_{1}\eta_{2}\eta_{3}}^{\nu_{1}\nu_{2}\nu_{3}}}{B_{1,1,0}^{1,1,0}(d)}\;=$
(22)
$\displaystyle\hphantom{+}\sum_{j=1-\nu_{1}}^{\nu_{2}-1}(-1)^{\Sigma\nu_{i}}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\sum_{k=0}^{\eta_{3}}{\eta_{3}\choose
k}\sum_{n_{1}>n_{2}>0}\Big{(}n_{2}^{d-\ell_{1}}n_{1}^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{1}}+n_{1}^{d-\ell_{1}}n_{2}^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{1}}\Big{)}$
$\displaystyle+\sum_{j=1-\nu_{1}}^{\nu_{3}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\sum_{k=0}^{\eta_{2}}{\eta_{2}\choose
k}(-1)^{k}\sum_{n_{1}>n_{2}>0}\Big{(}n_{1}^{d-\ell_{1}}+n_{2}^{d-\ell_{1}}\Big{)}\Big{(}n_{1}+n_{2}\Big{)}^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{1}}$
$\displaystyle+\sum_{j=1-\nu_{2}}^{\nu_{3}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\sum_{k=0}^{\eta_{1}}{\eta_{1}\choose
k}(-1)^{k}\sum_{n_{1}>n_{2}>0}\Big{(}n_{1}^{d-\ell_{2}}+n_{2}^{d-\ell_{2}}\Big{)}\Big{(}n_{1}+n_{2}\Big{)}^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{2}}$

where we have abbreviated $\ell_{p}\equiv{\Sigma\nu_{i}}-\eta_{p}-j-k$, and
the ${\eta_{p}\choose k}$ are binomial coefficients. The double sums in the
last two lines of eq. (22) have the form

$\displaystyle\sum_{n_{1}>n_{2}>0}\Big{(}\frac{1}{n_{1}^{\alpha}}+\frac{1}{n_{2}^{\alpha}}\Big{)}\frac{1}{(n_{1}+n_{2})^{\beta}}$
$\displaystyle\quad=\;\sum_{n_{1}>n_{2}>0}\frac{1}{n_{1}^{\alpha}}\frac{1}{(n_{1}+n_{2})^{\beta}}+\Big{\\{}\sum_{n_{1},n_{2}\in\mathbb{N}}-\sum_{n_{2}>n_{1}>0}-\sum_{n_{1}>0}\delta_{n_{2}-n_{1}}\Big{\\}}\frac{1}{n_{2}^{\alpha}}\frac{1}{(n_{1}+n_{2})^{\beta}}$
$\displaystyle\quad=\;\sum_{n_{1}>n_{2}>0}\Big{[}\frac{1}{n_{1}^{\alpha}}\frac{1}{(n_{1}+n_{2})^{\beta}}+\frac{1}{n_{2}^{\alpha}}\frac{1}{n_{1}^{\beta}}-\frac{1}{n_{1}^{\alpha}}\frac{1}{(n_{1}+n_{2})^{\beta}}\Big{]}-\sum_{n>0}\frac{1}{n^{\alpha}}\frac{1}{(2n)^{\beta}}$
$\displaystyle\quad=\;\zeta(\beta,\alpha)-\frac{1}{2^{\beta}}\,\zeta(\alpha+\beta)\;.$
(23)

To get to the third line, we have shifted $n_{1}\rightarrow n_{1}-n_{2}$ and
$n_{1}\leftrightarrow 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. (46). The double sum in the first line of eq.
(22) has the form

$\displaystyle\sum_{n_{1}>n_{2}>0}\Big{(}\frac{1}{n_{1}^{\alpha}}\frac{1}{n_{2}^{\beta}}+\frac{1}{n_{2}^{\alpha}}\frac{1}{n_{1}^{\beta}}\Big{)}$
$\displaystyle=$
$\displaystyle\zeta(\alpha)\,\zeta(\beta)-\zeta(\alpha+\beta)\;,$ (24)

where the solution follows from the shuffle identity eq. (50).

###  5.2 Proof of statement (b)

In a similar way, plugging eq. (14) into eq. (10), re-expressing the ‘unwanted
mass’ in the prefactor in terms of the two others, shifting summation indices
appropriately and using the coefficients’ symmetry,

$\displaystyle\frac{{\bar{H}}_{\eta_{3}\eta_{2}\eta_{1}}^{\nu_{3}\nu_{2}\nu_{1}}}{B_{1,1,0}^{1,1,0}(d)}$
$\displaystyle=$
$\displaystyle\sum_{j=1-\nu_{1}}^{\nu_{2}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\sum_{k=0}^{\eta_{3}}{\eta_{3}\choose
k}(-1)^{k-\eta_{3}}\sum_{n_{1}>n_{2}>0}n_{2}^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{1}}n_{1}^{d-\ell_{1}}$
(25) $\displaystyle+$
$\displaystyle\sum_{j=1-\nu_{1}}^{\nu_{3}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\sum_{k=0}^{\eta_{2}}{\eta_{2}\choose
k}(-1)^{k-\eta_{2}}\sum_{n_{1}>n_{2}>0}(n_{1}-n_{2})^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{1}}n_{1}^{d-\ell_{1}}$
$\displaystyle+$
$\displaystyle\sum_{j=1-\nu_{2}}^{\nu_{3}-1}(-1)^{\Sigma\nu_{i}}\,c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\sum_{k=0}^{\eta_{1}}{\eta_{1}\choose
k}\sum_{n_{1}>n_{2}>0}(n_{1}-n_{2})^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{2}}n_{2}^{d-\ell_{2}}\;,$
$\displaystyle\frac{{\bar{H}}_{\eta_{3}\eta_{1}\eta_{2}}^{\nu_{3}\nu_{1}\nu_{2}}}{B_{1,1,0}^{1,1,0}(d)}$
$\displaystyle=$
$\displaystyle\sum_{j=1-\nu_{1}}^{\nu_{2}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\sum_{k=0}^{\eta_{3}}{\eta_{3}\choose
k}(-1)^{k}\sum_{n_{1}>n_{2}>0}n_{2}^{d-\ell_{1}}n_{1}^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{1}}$
$\displaystyle+$
$\displaystyle\sum_{j=1-\nu_{1}}^{\nu_{3}-1}(-1)^{\Sigma\nu_{i}}\,c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\sum_{k=0}^{\eta_{2}}{\eta_{2}\choose
k}\sum_{n_{1}>n_{2}>0}(n_{1}-n_{2})^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{1}}n_{2}^{d-\ell_{1}}$
$\displaystyle+$
$\displaystyle\sum_{j=1-\nu_{2}}^{\nu_{3}-1}(-1)^{j}\,c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\sum_{k=0}^{\eta_{1}}{\eta_{1}\choose
k}(-1)^{k-\eta_{1}}\sum_{n_{1}>n_{2}>0}(n_{1}-n_{2})^{d-2{\Sigma\nu_{i}}+{\Sigma\eta_{i}}+\ell_{2}}n_{1}^{d-\ell_{2}}\;.$

The double sums needed here are given in eqs. (46), (47) and (48),
respectively, resulting in single and double zeta values. Combining with the
prefactors corresponding to eq. (8), we therefore have

$\displaystyle\frac{{\bar{H}}_{\eta_{3}\eta_{2}\eta_{1}}^{\nu_{3}\nu_{2}\nu_{1}}+(-1)^{\eta_{3}}{\bar{H}}_{\eta_{3}\eta_{1}\eta_{2}}^{\nu_{3}\nu_{1}\nu_{2}}+(-1)^{\eta_{2}}[{H}_{\eta_{2}\eta_{1}\eta_{3}}^{\nu_{2}\nu_{1}\nu_{3}}+{H}_{\eta_{1}\eta_{2}\eta_{3}}^{\nu_{1}\nu_{2}\nu_{3}}]}{(-1)^{\Sigma\nu_{i}}\,B_{1,1,0}^{1,1,0}(d)}$
(27)
$\displaystyle=\\!\\!\sum_{j=1-\nu_{1}}^{\nu_{2}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\sum_{k=0}^{\eta_{3}}{\eta_{3}\choose
k}\Big{\\{}(-1)^{\ell_{3}}\Big{[}\zeta(d_{1},e_{1})\\!+\\!\zeta(e_{1},d_{1})\Big{]}+(-1)^{\eta_{2}}\Big{[}\zeta(e_{1})\zeta(d_{1})\\!-\\!\zeta(e_{1}\\!+\\!d_{1})\Big{]}\Big{\\}}$
$\displaystyle+\\!\\!\sum_{j=1-\nu_{1}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\sum_{k=0}^{\eta_{2}}{\eta_{2}\choose
k}\Big{\\{}(-1)^{\ell_{2}}\Big{[}\zeta(d_{1},e_{1})\\!+\\!\zeta(e_{1},d_{1})\\!-\\!\frac{\zeta(e_{1}\\!+\\!d_{1})}{2^{e_{1}}}\Big{]}+(-1)^{\eta_{3}}\zeta(e_{1})\zeta(d_{1})\Big{\\}}$
$\displaystyle+\\!\\!\sum_{j=1-\nu_{2}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\sum_{k=0}^{\eta_{1}}{\eta_{1}\choose
k}\Big{\\{}\zeta(d_{2})\zeta(e_{2})+(-1)^{\ell_{2}}\Big{[}(-1)^{\Sigma\eta_{i}}\zeta(d_{2},e_{2})\\!+\\!\zeta(e_{2},d_{2})\\!-\\!\frac{\zeta(e_{2}\\!+\\!d_{2})}{2^{e_{2}}}\Big{]}\Big{\\}}.$

Here, $\ell_{p}={\Sigma\nu_{i}}-\eta_{p}-j-k$, $d_{p}=\ell_{p}-d$ and
$e_{p}=2{\Sigma\nu_{i}}-d-{\Sigma\eta_{i}}-\ell_{p}$. Note that the prefactor
of eq. (8) guarantees that ${\Sigma\eta_{i}}$ is even, so we can drop the
$(-1)^{\Sigma\eta_{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. (50):

$\displaystyle\frac{{\bar{H}}_{\eta_{3}\eta_{2}\eta_{1}}^{\nu_{3}\nu_{2}\nu_{1}}+(-1)^{\eta_{3}}{\bar{H}}_{\eta_{3}\eta_{1}\eta_{2}}^{\nu_{3}\nu_{1}\nu_{2}}+(-1)^{\eta_{2}}[{H}_{\eta_{2}\eta_{1}\eta_{3}}^{\nu_{2}\nu_{1}\nu_{3}}+{H}_{\eta_{1}\eta_{2}\eta_{3}}^{\nu_{1}\nu_{2}\nu_{3}}]}{(-1)^{\Sigma\nu_{i}}\,B_{1,1,0}^{1,1,0}(d)}$
(28)
$\displaystyle=\\!\\!\sum_{j=1-\nu_{1}}^{\nu_{2}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\sum_{k=0}^{\eta_{3}}{\eta_{3}\choose
k}\Big{\\{}\Big{[}(-1)^{\ell_{3}}+(-1)^{\eta_{2}}\Big{]}\Big{[}\zeta(e_{1})\zeta(d_{1})-\zeta(e_{1}+d_{1})\Big{]}\Big{\\}}$
$\displaystyle+\\!\\!\sum_{j=1-\nu_{1}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\sum_{k=0}^{\eta_{2}}\\!{\eta_{2}\choose
k}\\!\Big{\\{}\\!(-1)^{\ell_{2}}\Big{[}\zeta(d_{1})\zeta(e_{1})\\!-\\!\zeta(e_{1}\\!+\\!d_{1})\\!-\\!\frac{\zeta(e_{1}\\!+\\!d_{1})}{2^{e_{1}}}\Big{]}\\!+\\!(-1)^{\eta_{3}}\zeta(e_{1})\zeta(d_{1})\Big{\\}}$
$\displaystyle+\\!\\!\sum_{j=1-\nu_{2}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\sum_{k=0}^{\eta_{1}}{\eta_{1}\choose
k}\Big{\\{}\zeta(d_{2})\zeta(e_{2})+(-1)^{\ell_{2}}\Big{[}\zeta(d_{2})\zeta(e_{2})-\zeta(e_{2}+d_{2})-\frac{\zeta(e_{2}+d_{2})}{2^{e_{2}}}\Big{]}\Big{\\}}.$

###  5.3 Proof of statement (c)

Looking at the structure of eq. (28), we are happy with the products of zetas
(which we can in the end interpret as products of 1-loop sum-integrals eq.
(4)), and note that all remaining single zeta values have arguments
$e_{p}+d_{p}=2{\Sigma\nu_{i}}-2d-{\Sigma\eta_{i}}$ that do not depend on the
summation indices $j$, $k$. For these terms, the $k$-sum can be easily
performed using

$\displaystyle\sum_{k=0}^{N}{N\choose k}x^{k}$ $\displaystyle=$
$\displaystyle(1+x)^{N}\;,\quad\sum_{k=0}^{N}{N\choose
k}(-1)^{k}=\delta_{N}\;,$ (29)

producing

$\displaystyle-\frac{{\bar{H}}_{\eta_{3}\eta_{2}\eta_{1}}^{\nu_{3}\nu_{2}\nu_{1}}+(-1)^{\eta_{3}}{\bar{H}}_{\eta_{3}\eta_{1}\eta_{2}}^{\nu_{3}\nu_{1}\nu_{2}}+(-1)^{\eta_{2}}[{H}_{\eta_{2}\eta_{1}\eta_{3}}^{\nu_{2}\nu_{1}\nu_{3}}+{H}_{\eta_{1}\eta_{2}\eta_{3}}^{\nu_{1}\nu_{2}\nu_{3}}]}{(-1)^{\Sigma\nu_{i}}\,B_{1,1,0}^{1,1,0}(d)\,\zeta(2{\Sigma\nu_{i}}-2d-{\Sigma\eta_{i}})}\Big{|}_{\rm{single~{}zetas}}$
(30)
$\displaystyle\qquad=\sum_{j=1-\nu_{1}}^{\nu_{2}-1}c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\Big{[}(-1)^{{\Sigma\nu_{i}}-\eta_{3}-j}\delta_{\eta_{3}}+(-1)^{\eta_{2}}2^{\eta_{3}}\Big{]}$
$\displaystyle\qquad+\sum_{j=1-\nu_{1}}^{\nu_{3}-1}c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\Big{[}(-1)^{{\Sigma\nu_{i}}-\eta_{2}-j}\delta_{\eta_{2}}+(-1)^{{\Sigma\nu_{i}}-\eta_{2}-j}(\tfrac{1}{2})^{{\Sigma\nu_{i}}-d-\eta_{3}+j}\Big{]}$
$\displaystyle\qquad+\sum_{j=1-\nu_{2}}^{\nu_{3}-1}c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\Big{[}(-1)^{{\Sigma\nu_{i}}-\eta_{2}-j}\delta_{\eta_{1}}+(-1)^{{\Sigma\nu_{i}}-\eta_{2}-j}(\tfrac{1}{2})^{{\Sigma\nu_{i}}-d-\eta_{3}+j}\Big{]}\;.$

Comparing with eqs. (17)–(20), the remaining sums over $j$ can now be seen to
be nothing but the coefficients of special cases of $B$,

$\displaystyle{\bar{H}}_{\eta_{3}\eta_{2}\eta_{1}}^{\nu_{3}\nu_{2}\nu_{1}}+(-1)^{\eta_{3}}{\bar{H}}_{\eta_{3}\eta_{1}\eta_{2}}^{\nu_{3}\nu_{1}\nu_{2}}+(-1)^{\eta_{2}}[{H}_{\eta_{2}\eta_{1}\eta_{3}}^{\nu_{2}\nu_{1}\nu_{3}}+{H}_{\eta_{1}\eta_{2}\eta_{3}}^{\nu_{1}\nu_{2}\nu_{3}}]\Big{|}_{\rm{single~{}zetas}}$
(31)
$\displaystyle=-(-1)^{\Sigma\nu_{i}}\,B_{1,1,0}^{1,1,0}(d)\,\zeta(2{\Sigma\nu_{i}}-2d-{\Sigma\eta_{i}})\Big{\\{}(-1)^{{\Sigma\nu_{i}}-\eta_{3}}\delta_{\eta_{3}}\beta^{\nu_{3}\nu_{1}\nu_{2}}+(-1)^{\eta_{2}}2^{\eta_{3}}(-1)^{\Sigma\nu_{i}}r_{1}^{\nu_{1}\nu_{2}\nu_{3}}$
$\displaystyle\quad+(-1)^{{\Sigma\nu_{i}}-\eta_{2}}\delta_{\eta_{2}}\beta^{\nu_{2}\nu_{1}\nu_{3}}+(-1)^{{\Sigma\nu_{i}}-\eta_{2}}\delta_{\eta_{1}}\beta^{\nu_{1}\nu_{2}\nu_{3}}+(-1)^{{\Sigma\nu_{i}}-\eta_{2}}(\tfrac{1}{2})^{{\Sigma\nu_{i}}-d-\eta_{3}}\frac{2^{\Sigma\nu_{i}}}{4}r_{2}^{\nu_{1}\nu_{2}\nu_{3}}\Big{\\}}$
$\displaystyle=-\zeta(2{\Sigma\nu_{i}}-2d-{\Sigma\eta_{i}})\Big{[}\delta_{\eta_{3}}B_{0,1,1}^{\nu_{3}\nu_{1}\nu_{2}}+\delta_{\eta_{2}}B_{0,1,1}^{\nu_{2}\nu_{1}\nu_{3}}+(-1)^{\eta_{2}}\delta_{\eta_{1}}B_{0,1,1}^{\nu_{1}\nu_{2}\nu_{3}}+(-1)^{\eta_{2}}2^{\eta_{3}}B_{1,1,2}^{\nu_{1}\nu_{2}\nu_{3}}\Big{]}\;,$

which (after rewriting
$(-1)^{\eta_{2}}\delta_{\eta_{1}}=(-1)^{\Sigma\eta_{i}}(-1)^{\eta_{3}}\delta_{\eta_{1}}$
and noticing that ${\Sigma\eta_{i}}$ is even due to the prefactor of eq. (8))
cancels exactly against the second line of eq. (8). This completes the proof
that only products of single zeta values remain in the full result for our
sum-integral $L$.

###  5.4 Proof of statement (d)

All that we are left with for the curly brackets of eq. (8) are the products
of zeta values in eq. (28). Using that ${\Sigma\eta_{i}}$ is even (or
$(-1)^{\Sigma\eta_{i}}=1$), we can re-write the specific combinations of
prefactors in a more convenient form, to finally obtain
($\ell_{p}\equiv{\Sigma\nu_{i}}-\eta_{p}-j-k$)

$\displaystyle L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$
$\displaystyle=$ $\displaystyle\frac{T^{2}[1+(-1)^{\Sigma\eta_{i}}]}{(2\pi
T)^{2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-2d}}\,B_{1,1,0}^{1,1,0}(d)\,(-1)^{\Sigma\nu_{i}}\,\times$
$\displaystyle\Big{\\{}$
$\displaystyle\sum_{j=1-\nu_{1}}^{\nu_{2}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\sum_{k=0}^{\eta_{3}}\\!{\eta_{3}\choose
k}(-1)^{\eta_{2}}\,[1+(-1)^{\ell_{1}}]\,\zeta(\ell_{1}\\!-\\!d)\,\zeta(2{\Sigma\nu_{i}}\\!-\\!{\Sigma\eta_{i}}\\!-\\!\ell_{1}\\!-\\!d)$
$\displaystyle+$
$\displaystyle\sum_{j=1-\nu_{1}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\sum_{k=0}^{\eta_{2}}\\!{\eta_{2}\choose
k}(-1)^{\eta_{3}}\,[1+(-1)^{\ell_{1}}]\,\zeta(\ell_{1}\\!-\\!d)\,\zeta(2{\Sigma\nu_{i}}\\!-\\!{\Sigma\eta_{i}}\\!-\\!\ell_{1}\\!-\\!d)$
$\displaystyle+$
$\displaystyle\sum_{j=1-\nu_{2}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\sum_{k=0}^{\eta_{1}}\\!{\eta_{1}\choose
k}[1+(-1)^{\ell_{2}}]\,\zeta(\ell_{2}-d)\,\zeta(2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-\ell_{2}-d)\Big{\\}}\;.$

We recall from eq. (13) that the normalization factor
$B_{1,1,0}^{1,1,0}(d)=\frac{\Gamma^{2}(1-d/2)}{(4\pi)^{d}}=[J_{1}(d,1)]^{2}$
is a 1-loop massive tadpole squared (taken at $m=1$), cf. eq. (3). Note that
the factors $[1+(-1)^{\Sigma\eta_{i}}]$ and $[1+(-1)^{\ell_{p}}]$ ensure that
only $\zeta(n_{\rm even}-d)$ at even integers $n$ contribute to the right-hand
side of eq. (5.4).

While the symmetries of the integral $L$ are not made explicit in eq. (5.4),
it is easy to verify that they are indeed satisfied. For example,
interchanging the first and second set of indices,
$(\nu_{1},\eta_{1})\leftrightarrow(\nu_{2},\eta_{2})$, induces an overall sign
depending on the parity of $\eta_{3}$, see eq. (36). To check this, one first
observes that under this index permutation the roles of two last two lines
interchange, while for the remaining line in the curly brackets one applies
the first of eq. (15) followed by a change of summation indices
$j\rightarrow-j$ and $k\rightarrow\eta_{3}-k$.

##  6 Result and extension

Let us finally re-write our result as products of 1-loop sum-integrals
$I_{\nu}^{\eta}(d,T)$ from eq. (4). Each of the zeta functions in eq. (5.4)
will map onto one such 1-loop sum-integral $I_{\nu}^{\eta}(d,T)$, and we need
to distinguish the cases $\zeta(2n-d)$ with $n>0$ (giving master sum-integrals
$I_{n}^{0}$ without numerator) and $n\leq 0$ (giving master sum-integrals
$I_{1}^{2-2n}$ with numerators). Both cases can be covered by the single
expression

$\displaystyle\zeta(2n\\!-\\!d)=\frac{(2\pi
T)^{2n-d}}{T\,J_{1}(d,1)}\,\hat{I}_{n}(d,T)\,,\;\;\hat{I}_{n}(d,T)\equiv\frac{\Gamma(\sigma_{n})\,I_{\sigma_{n}}^{2\sigma_{n}-2n}(d,T)}{2\,\big{(}1-\tfrac{d}{2}\big{)}_{\sigma_{n}-1}}\,,\;\;\sigma_{n}\equiv{\rm
max}(n,1)\,,\quad$ (33)

with $J_{1}(d,1)$ the 1-loop tadpole of eq. (3). Doing these identifications,
we get (recall the abbreviations ${\Sigma\nu_{i}}=\nu_{1}+\nu_{2}+\nu_{3}$,
${\Sigma\eta_{i}}=\eta_{1}+\eta_{2}+\eta_{3}$ and
$\ell_{p}={\Sigma\nu_{i}}-\eta_{p}-j-k$)

$\displaystyle L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$
$\displaystyle=$
$\displaystyle[1+(-1)^{\Sigma\eta_{i}}]\,(-1)^{\Sigma\nu_{i}}\,\times$ (34)
$\displaystyle\Big{\\{}$
$\displaystyle\\!\\!\\!\sum_{j=1-\nu_{1}}^{\nu_{2}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{2};j}(d)\sum_{k=0}^{\eta_{3}}\\!{\eta_{3}\choose
k}(-1)^{\eta_{2}}\,[1\\!+\\!(-1)^{\ell_{1}}]\,\hat{I}_{{\ell_{1}}/2}(d,T)\,\hat{I}_{({2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-\ell_{1}})/2}(d,T)$
$\displaystyle+$
$\displaystyle\\!\\!\\!\sum_{j=1-\nu_{1}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{1},\nu_{3};j}(d)\sum_{k=0}^{\eta_{2}}\\!{\eta_{2}\choose
k}(-1)^{\eta_{3}}\,[1\\!+\\!(-1)^{\ell_{1}}]\,\hat{I}_{{\ell_{1}}/2}(d,T)\,\hat{I}_{({2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-\ell_{1}})/2}(d,T)$
$\displaystyle+$
$\displaystyle\\!\\!\\!\sum_{j=1-\nu_{2}}^{\nu_{3}-1}\\!\\!\\!c^{({\Sigma\nu_{i}})}_{\nu_{2},\nu_{3};j}(d)\sum_{k=0}^{\eta_{1}}\\!{\eta_{1}\choose
k}[1\\!+\\!(-1)^{\ell_{2}}]\,\hat{I}_{{\ell_{2}}/2}(d,T)\,\hat{I}_{({2{\Sigma\nu_{i}}-{\Sigma\eta_{i}}-\ell_{2}})/2}(d,T)\Big{\\}},\quad$

which, together with eq. (4), constitutes our final result. Sample Mathematica
code that implements our final result is given in appendix B.

###  6.1 Non-positive indices and full reduction algorithm

Having the decomposition eq. (34) for scalar bosonic two-loop vacuum sum-
integrals $L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$ with
positive propagator powers $\nu_{i}\in\mathbb{N}$ at hand, we can easily
extend it to be valid for all non-negative integers $\nu_{i}\in\mathbb{N}_{0}$
by analyzing the boundary cases when (at least) one $\nu_{i}=0$. Analyzing eq.
(3.24) of Davydychev:2022dcw  with boundary conditions eqs. (2.7)–(2.9)
therein, it is sufficient to amend the definition of the coefficient functions
given in eq. (4) above with the specific cases

$\displaystyle c^{(\nu_{a}+\nu_{b})}_{\nu_{a},\nu_{b};\nu_{b}-\nu_{a}}(d)$
$\displaystyle\equiv$
$\displaystyle(-1)^{\nu_{a}+\nu_{b}}\frac{\big{(}1-\frac{d}{2}\big{)}_{\nu_{a}-1}}{\Gamma(\nu_{a})}\frac{\big{(}1-\frac{d}{2}\big{)}_{\nu_{b}-1}}{\Gamma(\nu_{b})}\;,$
(35)

while the remaining boundary coefficients
($c^{({\Sigma\nu_{i}})}_{\nu_{a},0;j}=0=c^{({\Sigma\nu_{i}})}_{0,\nu_{b};j}$
and $c^{(\nu_{a}+\nu_{b})}_{\nu_{a},\nu_{b};j\neq\nu_{b}-\nu_{a}}=0$) are
already correctly reproduced by eq. (4).

With a little bit of additional work, a further useful extension to all
integer powers $\nu_{i}\in\mathbb{Z}$ can be done as follows. Given an
$L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$ as defined in
eq. (5), with $\eta_{i}\in\mathbb{N}_{0}$ and $\nu_{i}\in\mathbb{Z}$ with at
least one non-positive index $\nu_{i}\leq 0$, apply

$\displaystyle
L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)\;=\;L_{\nu_{1},\nu_{3},\nu_{2}}^{\eta_{1},\eta_{3},\eta_{2}}(d,T)\quad,\quad
L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)\;=\;(-1)^{\eta_{3}}\,L_{\nu_{2},\nu_{1},\nu_{3}}^{\eta_{2},\eta_{1},\eta_{3}}(d,T)$
(36)

(following from momentum shifts $Q\rightarrow P-Q$ and $P\leftrightarrow Q$,
respectively) to order the lower indices such that
$\nu_{1}\geq\nu_{2}\geq\nu_{3}$. Then, simplify and factorize those cases by
(a) realizing that integrals with odd upper-index sum ${\Sigma\eta_{i}}=$odd
vanish due to the denominator’s invariance under the shift
$\\{P,Q\\}\rightarrow\\{-P,-Q\\}$;
(b) if $\eta_{3}>0$, expanding the corresponding numerator binomially as in
eq. (44);
(c) if $\nu_{3}<0$, expand the corresponding numerator double-binomially to
map onto 1-loop tensor tadpoles
$T_{s,t}^{\mu_{1}\dots\mu_{n}}={{\Sigma}\\!\\!\\!\\!\raise
1.29167pt\hbox{$\int$}_{\mbox{\tiny\rm{${P}$}}}}\frac{(U\cdot
P)^{t}}{[P^{2}]^{s}}\,P^{\mu_{1}}\cdots P^{\mu_{n}}$, use the 1-loop tensor
tadpole decomposition into orthogonal bases made out of $U^{\mu}=(1,\bf{0})$
and $g^{\mu\nu}$ as given in Navarrete:2022adz  and compute the basis tensor
contractions; leading to

$\displaystyle L_{{\nu_{1},\nu_{2},\nu_{3}}\leq
0}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$ $\displaystyle=$
$\displaystyle\frac{1}{2}[1+(-1)^{\Sigma\eta_{i}}]\sum_{j=0}^{\eta_{3}}\sum_{n=0}^{-\nu_{3}}\sum_{k=0}^{-\nu_{3}-n}\sum_{\ell=0}^{\lfloor
n/2\rfloor}\frac{(-1)^{j}\,{\eta_{3}\choose
j}\,(-2)^{n}\,(-\nu_{3})!}{4^{\ell}\,\ell!\,(n-2\ell)!\,k!\,(-\nu_{3}-n-k)!}\;\times$
(37) $\displaystyle\times$
$\displaystyle\frac{\big{(}\frac{d}{2}\big{)}_{\ell}\;I_{\nu_{1}+\nu_{3}+n+k}^{\eta_{1}+\eta_{3}-j+n}(d,T)\cdot
I_{\nu_{2}-k}^{\eta_{2}+j+n}(d,T)}{\big{(}\frac{d}{2}+1-\nu_{2}+k\big{)}_{\ell}\big{(}\frac{d}{2}+1-\nu_{1}-\nu_{3}-k-n\big{)}_{\ell}}\;.$

The 1-loop tadpoles $I$ in eq. (37) can finally be mapped onto the basis
$\hat{I}$ via eqs. (4) and (33). On the one hand, this introduces Gamma
factors in the denominators which correctly set all $I_{\nu\leq 0}^{\eta}=0$.
To exclude such vanishing terms from the outset, we can also modify the limits
of the $k$-sum as $\sum_{k={\rm max}(0,1-\nu_{1}-\nu_{3}-n)}^{{\rm
min}(-\nu_{3}-n,\nu_{2}-1)}$. On the other hand, the map introduces prefactors
which implement that $I_{\nu}^{\eta={\rm odd}}=0$. Modifying the $j$-sum as
$\sum_{j=0}^{\eta_{3}}f(j)\rightarrow\sum_{j=x}^{\lfloor(\eta_{3}+x)/2\rfloor}f(2j-x)$
where $x={\rm mod}(\eta_{2}+n,2)$ then eliminates the need for such explicit
prefactors $[1+(-1)^{\eta}]$ stemming from eq. (4). After some slight
simplifications, this leaves us with a representation on par with eq. (34):

$\displaystyle L_{{\nu_{1},\nu_{2},\nu_{3}}\leq
0}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$ $\displaystyle=$
$\displaystyle\frac{1}{2}[1+(-1)^{\Sigma\eta_{i}}]\sum_{n=0}^{-\nu_{3}}\sum_{j=x}^{\lfloor(\eta_{3}+x)/2\rfloor}\sum_{k={\rm
max}(1,n+{\Sigma\nu_{i}}-\nu_{1})}^{{\rm
min}(\nu_{2},n+{\Sigma\nu_{i}}-1)}\sum_{\ell=0}^{\lfloor n/2\rfloor}\;\times$
(38) $\displaystyle\times$
$\displaystyle\frac{(-1)^{\eta_{2}}\,(-\nu_{3})!\,2^{2+n-2\ell}\,{\eta_{3}\choose
2j-x}}{\ell!\,(n-2\ell)!\,(\nu_{2}-k)!\,(\nu_{1}-z)!\,\Gamma(z)\,\Gamma(k)}\;\times$
$\displaystyle\times$
$\displaystyle\frac{\big{(}\frac{d}{2}\big{)}_{\ell}\,(1-\frac{d}{2})_{z-1}\,(1-\frac{d}{2})_{k-1}}{\big{(}\frac{d}{2}+1-z\big{)}_{\ell}\,\big{(}\frac{d}{2}+1-k\big{)}_{\ell}}\;\hat{I}_{{\Sigma\nu_{i}}-{\Sigma\eta_{i}}/2-y}(d,T)\cdot\hat{I}_{y}(d,T)\;,$

where $x\equiv{\rm mod}(\eta_{2}+n,2)$ as well as $y\equiv
k-j-\lfloor(\eta_{2}+n)/2\rfloor$ and $z\equiv{\Sigma\nu_{i}}+n-k$. As a
check, we note that for $\nu_{3}=0$ and positive values of $\nu_{1},\nu_{2}$,
eq. (38) as well as eq. (34) with eq. (35) both reduce to single sums555In
particular, the quadruple sum of eq. (38) collapses to a single sum since only
the terms at $n=0$, $k=\nu_{2}$ and $\ell=0$ contribute; while for eq. (34),
only the first of the three double-sums contributes, with its $j$-sum
collapsing to a single term $j=\nu_{2}-\nu_{1}$., and agree in their
factorization of $L_{\nu_{1},\nu_{2},0}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$. We
can therefore use the much more general eq. (38) for reducing sum-integrals
involving non-positive indices, without ever applying eq. (35). In order to
check eq. (38), we have verified it against a pre-existing in-house IBP code,
for a large number of integer values for $\nu_{i}$ and $\eta_{i}$, finding
complete agreement.

We are now in a position to formulate an algorithm that maps every massless
scalar bosonic two-loop vacuum sum-integral to its 1-loop factorized form:
Given an $L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$ as
defined in eq. (5), with $\eta_{i}\in\mathbb{N}_{0}$ and
$\nu_{i}\in\mathbb{Z}$,

  1. 1.

If all $\nu_{i}$ are positive integers, apply eq. (34) to map
$L\rightarrow\sum\hat{I}\cdot\hat{I}$, where the $\hat{I}$ are elements of the
1-loop basis defined in eq. (33).

  2. 2.

If any $\nu_{i}\leq 0$, put the lower indices into descending order by eq.
(36) and then apply eq. (38) to map $L\rightarrow\sum\hat{I}\cdot\hat{I}$.

For sample Mathematica code implementing this strategy, see appendix B.

###  6.2 Checks and examples

For checks, let us reproduce some known IBP results, using our main result eq.
(34). There are short IBP proofs in eq. (5) of Schroder:2008ex  (see also eq.
(A.3) of Nishimura:2012ee ) and eq. (E.2) of Ghisoiu:2012kn , respectively, of
the (propagator-) weight 3 examples

$\displaystyle L_{1,1,1}^{0,0,0}$ $\displaystyle=$ $\displaystyle 0\;,\quad
L_{2,1,1}^{2,0,0}\;=\;0\;.$ (39)

In appendix B of Ghisoiu:2012yk , we gave the IBP reductions (at weight 3, 4,
5 and 7, respectively)

$\displaystyle L_{2,1,1}^{0,2,0}$ $\displaystyle=$
$\displaystyle\frac{(d-3)}{(d-5)}\,I_{2}^{0}\,I_{1}^{0}\;,\quad
L_{3,1,1}^{2,2,0}\;=\;-\frac{(d-4)(d^{2}-8d+19)}{4(d-7)(d-5)}\,I_{2}^{0}\,I_{1}^{0}\;,$
(40) $\displaystyle L_{2,1,1}^{0,0,0}$ $\displaystyle=$
$\displaystyle-\frac{I_{2}^{0}\,I_{2}^{0}}{(d-5)(d-2)}\;,$ (41) $\displaystyle
L_{3,1,1}^{0,0,0}$ $\displaystyle=$
$\displaystyle-\frac{4\,I_{3}^{0}\,I_{2}^{0}}{(d-7)(d-2)}\;,\;\;L_{2,2,1}^{0,0,0}\;=\;0\;,\;\;L_{3,1,2}^{1,1,0}\;=\;\frac{(d-6)(d-5)(d-3)}{2(d-9)(d-7)(d-2)}\,I_{3}^{0}\,I_{2}^{0}\;,\qquad$
(42) $\displaystyle L_{3,3,1}^{0,0,0}$ $\displaystyle=$
$\displaystyle-\frac{12(d-8)(d-5)\,I_{4}^{0}\,I_{3}^{0}}{(d-11)(d-9)(d-4)(d-2)}\;.$
(43)

Using the Mathematica code of appendix B, we confirm all IBP reductions listed
above.

As another check, we can expand the numerator of eq. (5) for positive values
of the index $\eta_{3}$, and verify that each $L$ can be expressed in terms of
instances with $\eta_{3}=0$ as

$\displaystyle L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1},\eta_{2},\eta_{3}}(d,T)$
$\displaystyle=$ $\displaystyle\sum_{j=0}^{\eta_{3}}{\eta_{3}\choose
j}\,(-1)^{j}\,L_{\nu_{1},\nu_{2},\nu_{3}}^{\eta_{1}+\eta_{3}-j,\eta_{2}+j,0}(d,T)\;.$
(44)

Indeed, testing eq. (44) for a wide range of integers
$\\{{\nu_{1},\nu_{2},\nu_{3}},{\eta_{1},\eta_{2},\eta_{3}}\\}$ by applying eq.
(34) on both sides, we find perfect agreement.

Going beyond published results, one could now use eq. (34) 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

$\displaystyle L_{1,1,4}^{2,2,0}$ $\displaystyle=$
$\displaystyle-\frac{(d-4)\,I_{2}^{0}\,I_{2}^{0}}{4(d-9)(d-7)(d-2)}-\frac{(d-5)(d-4)\,I_{1}^{0}\,I_{3}^{0}}{2(d-9)(d-7)}+\frac{(2d-13)\,I_{1}^{2}\,I_{4}^{0}}{(d-9)}\;,$
$\displaystyle L_{1,1,4}^{0,0,0}$ $\displaystyle=$
$\displaystyle-\frac{4\,I_{3}^{0}\,I_{3}^{0}}{(d-9)(d-7)(d-4)(d-2)}-\frac{6\,I_{2}^{0}\,I_{4}^{0}}{(d-9)(d-2)}\;,$
$\displaystyle L_{2,2,4}^{0,0,0}$ $\displaystyle=$
$\displaystyle\frac{36\left(7d^{2}-105d+382\right)\,I_{4}^{0}\,I_{4}^{0}}{(d-13)(d-11)(d-9)(d-6)(d-4)(d-2)}+\frac{16(d-10)(d-5)\,I_{3}^{0}\,I_{5}^{0}}{(d-13)(d-11)(d-4)(d-2)}\;.$

On the other hand, from eqs. (36) and (37) some examples for non-positive
indices read

$\displaystyle L_{1,2,-2}^{2,1,1}$ $\displaystyle=$ $\displaystyle
2(d-3)\,I_{1}^{2}\,I_{1}^{4}\;,$ $\displaystyle L_{2,-3,3}^{1,2,3}$
$\displaystyle=$
$\displaystyle(d^{3}-14d^{2}+62d-94)\,I_{1}^{2}\,I_{1}^{4}+3(d-2)\,I_{1}^{0}\,I_{1}^{6}\;,$
$\displaystyle L_{-4,2,3}^{2,5,3}$ $\displaystyle=$ $\displaystyle
8(d^{2}-4d+10)\,I_{1}^{6}\,I_{1}^{6}-2(d-2)(d^{2}-6d+33)\,I_{1}^{4}\,I_{1}^{8}+4(d-2)(d-1)\,I_{1}^{2}\,I_{1}^{10}\;.$

##  7 Conclusions and outlook

We have presented an analytic solution for two-loop scalar massless bosonic
vacuum sum-integrals $L$, with general integer values of propagator and
numerator powers. Our key result presented in eq. (34) 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 Davydychev:2022dcw .

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 $P^{\mu_{1}}\cdots
P^{\mu_{k}}\,Q^{\mu_{k+1}}\cdots Q^{\mu_{n}}$ in the numerator of eq. (5). If
such tensors are fully symmetric in the $\mu_{i}$, one can trivially apply the
decomposition presented in eqs. (8), (16) of Navarrete:2022adz . 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 $\sim L\cdot 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 ${\mu_{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 Kapusta:2023eix ) or practical integral representations and
evaluations of some terms in their $\varepsilon$-expansion (see e.g. appendix
A of Laine:2019uua  as well as Laine:2017hdk ).

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.

###### Acknowledgements.

We wish to thank I. Kondrashuk for discussions at early stages of this work.
The work of A.D. was partially supported by CONICYT PCI/MEC 80180071. P.N. has
been supported by an ANID grant Magíster Nacional No. 22211544, and by the
Academy of Finland grants No. 347499, 353772 and 1322507. Y.S. acknowledges
support from ANID under FONDECYT projects No. 1191073 and 1231056.

##  Appendix A Single and double sums

Other than single sums

$\displaystyle\sum_{n=1}^{\infty}\frac{1}{n^{s}}$ $\displaystyle=$
$\displaystyle\zeta(s)\;,$ (45)

we encounter double sums of the following types:

$\displaystyle\sum_{n_{1}>n_{2}>0}\frac{1}{n_{1}^{s_{1}}\,n_{2}^{s_{2}}}$
$\displaystyle\equiv$ $\displaystyle\zeta(s_{1},s_{2})\;,$ (46)
$\displaystyle\sum_{n_{1}>n_{2}>0}\frac{1}{n_{1}^{s_{1}}\,(n_{1}-n_{2})^{s_{2}}}$
$\displaystyle=$ $\displaystyle\zeta(s_{1},s_{2})\;,$ (47)
$\displaystyle\sum_{n_{1}>n_{2}>0}\frac{1}{(n_{1}-n_{2})^{s_{1}}\,n_{2}^{s_{2}}}$
$\displaystyle=$ $\displaystyle\zeta(s_{1})\,\zeta(s_{2})\;.$ (48)

For the multiple zeta function $\zeta(s_{1},s_{2})$ occurring above, the
shuffle identity (see also figure 2)

$\displaystyle\sum_{n_{1}=1}^{\infty}\,\sum_{n_{2}=1}^{\infty}\frac{1}{n_{1}^{s_{1}}\,n_{2}^{s_{2}}}$
$\displaystyle=$
$\displaystyle\Big{(}\sum_{n_{1}>n_{2}>0}+\sum_{n_{2}>n_{1}>0}+\sum_{n_{1}=1}^{\infty}\delta_{n_{2}-n_{1}}\Big{)}\frac{1}{n_{1}^{s_{1}}\,n_{2}^{s_{2}}}$
(49) $\displaystyle\Leftrightarrow\qquad\quad\zeta(s_{1})\,\zeta(s_{2})$
$\displaystyle=$
$\displaystyle\zeta(s_{1},s_{2})+\zeta(s_{2},s_{1})+\zeta(s_{1}+s_{2})$ (50)

can be used to sort the two arguments.

We also encounter two types of unknown double sums

$\displaystyle\sum_{n_{1}>n_{2}>0}\frac{1}{(n_{1}+n_{2})^{s_{1}}\,n_{2}^{s_{2}}}$
$\displaystyle\equiv$ $\displaystyle Z(s_{1},s_{2})\;,$ (51)
$\displaystyle\sum_{n_{1}>n_{2}>0}\frac{1}{(n_{1}+n_{2})^{s_{1}}\,n_{1}^{s_{2}}}$
$\displaystyle=$
$\displaystyle\zeta(s_{1},s_{2})-Z(s_{1},s_{2})-\frac{1}{2^{s_{1}}}\,\zeta(s_{1}+s_{2})\;,$
(52)

which in eq. (8) fortunately always occur in combinations where the unknown
part $Z$ drops out, as proven in section 5.1.

$n_{1}$$n_{2}$$=$$n_{1}$$n_{2}$$+$$n_{1}$$n_{2}$$+$$n_{1}$$n_{2}$ Figure 2:
For the proof of the shuffle relation eq. (50), the sum over all pairs of
positive integers in the first quadrant is replaced by the sum over the three
regions on the right-hand side.

##  Appendix B Code snippet: Reduction of 2-loop sum-integral $L$

Here, we show sample Mathematica code that implements our eqs. (33), (34),
(36) and (38). It provides the six-index function ${\tt L}[{\tt n}_{1},{\tt
n}_{2},{\tt n}_{3},{\tt e}_{1},{\tt e}_{2},{\tt
e}_{3}]=L_{n_{1},n_{2},n_{3}}^{e_{1},e_{2},e_{3}}(d,T)$ which evaluates to a
linear combination of products of one-loop master sum-integrals ${\tt
i[s,t]}=I_{s}^{t}(d,T)$, with coefficients that are rational functions in the
number of spatial dimensions $d$.

    ClearAll[d,i,c,h,L];
    c[nu_,na_,nb_,j_] := Module[ {nj=Ceiling[(nu+j)/2], nc=nu-na-nb},
         (-1)^(nu-nj+1) Pochhammer[1-d/2,nj-j-1] / (2 Pochhammer[1/2,nj-nb-nc]
         Pochhammer[1/2,nj-j-na-nc] (nc-1)!) Sum[ Pochhammer[d/2-nj+1,k-1]
         (nj-k-1)! / (Pochhammer[d/2+3/2-nu,nj-k] Pochhammer[1/2,nc-nj+k]
         (k-1)! (k-j-1)! (nb-k)! (na-k+j)!), {k,Max[1+j,1],Min[na+j,nb]}] ];
    h[n_] := Module[ {s=Max[n,1]},
         Gamma[s] i[s,2s-2n]/2/Pochhammer[1-d/2,s-1] ];
    L[n1_,n2_,n3_,e1_,e2_,e3_] := Which[ n1<n2, (-1)^e3 L[n2,n1,n3,e2,e1,e3],
         n2<n3, L[n1,n3,n2,e1,e3,e2], True, Collect[ Module[ {nu=n1+n2+n3,
         eta=e1+e2+e3}, (1+(-1)^eta) If[ n3>0,
         (*all indices n1..n3 positive: apply 6.2*)
         Sum[ c[nu,n1,n2,j] Sum[ With[ {el=nu-e1-j-k}, Binomial[e3,k] (-1)^e2
         (1+(-1)^el) h[el/2] h[nu-eta/2-el/2] ], {k,0,e3}], {j,1-n1,n2-1}] +
         Sum[ c[nu,n1,n3,j] Sum[ With[ {el=nu-e1-j-k}, Binomial[e2,k] (-1)^e3
         (1+(-1)^el) h[el/2] h[nu-eta/2-el/2] ], {k,0,e2}], {j,1-n1,n3-1}] +
         Sum[ c[nu,n2,n3,j] Sum[ With[ {el=nu-e2-j-k}, Binomial[e1,k]
         (1+(-1)^el) h[el/2] h[nu-eta/2-el/2] ], {k,0,e1}], {j,1-n2,n3-1}],
         (*not all indices n1..n3 positive: apply 6.12*)
         2 (-1)^e2 (-n3)! Sum[ With[ {x=Mod[e2+n,2]},
         Sum[ With[ {y=k-j-Floor[(e2+n)/2], z=nu+n-k},
         Binomial[e3,2j-x] 2^n/(l!(n-2l)!(n2-k)!(n1-z)!Gamma[z]Gamma[k]4^l)
         Pochhammer[d/2,l] Pochhammer[1-d/2,z-1] Pochhammer[1-d/2,k-1]/
         (Pochhammer[d/2+1-z,l] Pochhammer[d/2+1-k,l]) h[nu-eta/2-y] h[y]],
         {j,x,Floor[(e3+x)/2]}, {k,Max[1,n+nu-n1],Min[n2,n+nu-1]},
         {l,0,Floor[n/2]} ] ], {n,0,-n3}]] ], i[__], Factor ] ];

After loading these definitions, one can simply call the function L to obtain
the factorized version of the corresponding sum-integral. As an example, the
result for $L_{3,1,1}^{2,2,0}(d,T)$ shown in eq. (40) is generated by
executing

    (*example use*)
    L[3,1,1,2,2,0]

## References

  * (1) F. Gross, E. Klempt, S. J. Brodsky, A. J. Buras, V. D. Burkert, G. Heinrich, K. Jakobs, C. A. Meyer, K. Orginos and M. Strickland, et al. 50 Years of Quantum Chromodynamics, [arXiv:2212.11107].
  * (2) J. I. Kapusta and C. Gale, Finite-Temperature Field Theory, Cambridge University Press, 2023.
  * (3) M. L. Bellac, Thermal Field Theory, Cambridge University Press, 2011.
  * (4) M. Laine and A. Vuorinen, Basics of Thermal Field Theory, Lect. Notes Phys. 925, Springer, 2016, [arXiv:1701.01554].
  * (5) J. Ghiglieri, A. Kurkela, M. Strickland and A. Vuorinen, Perturbative Thermal QCD: Formalism and Applications, Phys. Rept. 880 (2020), 1 [arXiv:2002.10188].
  * (6) P. B. Arnold and C. X. Zhai, The three loop free energy for pure gauge QCD, Phys. Rev. D 50 (1994), 7603 [arXiv:hep-ph/9408276].
  * (7) P. B. Arnold and C. X. Zhai, The three loop free energy for high temperature QED and QCD with fermions, Phys. Rev. D 51 (1995), 1906 [arXiv:hep-ph/9410360].
  * (8) G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lütgemeier and B. Petersson, Thermodynamics of SU(3) lattice gauge theory, Nucl. Phys. B 469 (1996), 419-444 [arXiv:hep-lat/9602007].
  * (9) S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Precision SU(3) lattice thermodynamics for a large temperature range, JHEP 07 (2012), 056 [arXiv:1204.6184].
  * (10) L. Giusti and M. Pepe, Equation of state of the SU(3) Yang–Mills theory: A precise determination from a moving frame, Phys. Lett. B 769 (2017), 385-390 [arXiv:1612.00265].
  * (11) E. Braaten and R. D. Pisarski, Resummation and Gauge Invariance of the Gluon Damping Rate in Hot QCD, Phys. Rev. Lett. 64 (1990), 1338.
  * (12) E. Braaten and R. D. Pisarski, Soft Amplitudes in Hot Gauge Theories: A General Analysis, Nucl. Phys. B 337 (1990), 569.
  * (13) J. O. Andersen, E. Braaten and M. Strickland, Hard thermal loop resummation of the free energy of a hot gluon plasma, Phys. Rev. Lett. 83 (1999), 2139 [arXiv:hep-ph/9902327].
  * (14) N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland and N. Su, Three-loop HTLpt thermodynamics at finite temperature and chemical potential, JHEP 05 (2014), 027 [arXiv:1402.6907].
  * (15) J. O. Andersen, N. Haque, M. G. Mustafa and M. Strickland, Three-loop hard-thermal-loop perturbation theory thermodynamics at finite temperature and finite baryonic and isospin chemical potential, Phys. Rev. D 93 (2016) no.5, 054045 [arXiv:1511.04660].
  * (16) F. Karsch, A. Patkos and P. Petreczky, Screened perturbation theory, Phys. Lett. B 401 (1997), 69 [arXiv:hep-ph/9702376].
  * (17) P. H. Ginsparg, First Order and Second Order Phase Transitions in Gauge Theories at Finite Temperature, Nucl. Phys. B 170 (1980), 388.
  * (18) T. Appelquist and R. D. Pisarski, High-Temperature Yang-Mills Theories and Three-Dimensional Quantum Chromodynamics, Phys. Rev. D 23 (1981), 2305.
  * (19) K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, Generic rules for high temperature dimensional reduction and their application to the standard model, Nucl. Phys. B 458 (1996), 90 [arXiv:hep-ph/9508379].
  * (20) E. Braaten and A. Nieto, Free energy of QCD at high temperature, Phys. Rev. D 53 (1996), 3421 [arXiv:hep-ph/9510408].
  * (21) K. Kajantie, M. Laine, K. Rummukainen and Y. Schröder, The Pressure of hot QCD up to $g^{6}\ln(1/g)$, Phys. Rev. D 67 (2003), 105008 [arXiv:hep-ph/0211321].
  * (22) A. Vuorinen, The Pressure of QCD at finite temperatures and chemical potentials, Phys. Rev. D 68 (2003), 054017 [arXiv:hep-ph/0305183].
  * (23) T. Gorda, R. Paatelainen, S. Säppi and K. Seppänen, Equation of State of Cold Quark Matter to $O(\alpha_{s}^{3}\ln\alpha_{s})$, Phys. Rev. Lett. 131 (2023) no.18, 181902 [arXiv:2307.08734].
  * (24) F. Di Renzo, M. Laine, V. Miccio, Y. Schröder and C. Torrero, The Leading non-perturbative coefficient in the weak-coupling expansion of hot QCD pressure, JHEP 07 (2006), 026 [arXiv:hep-ph/0605042].
  * (25) M. Laine and Y. Schröder, Quark mass thresholds in QCD thermodynamics, Phys. Rev. D 73 (2006), 085009 [arXiv:hep-ph/0603048].
  * (26) Y. Schröder, Weak-coupling expansion of the hot QCD pressure, PoS JHW2005 (2006), 029 [arXiv:hep-ph/0605057].
  * (27) P. Navarrete and Y. Schröder, Tackling the infamous $g^{6}$ term of the QCD pressure, PoS LL2022 (2022), 014 [arXiv:2207.10151].
  * (28) J. Österman, P. Schicho and A. Vuorinen, Integrating by parts at finite density, JHEP 08 (2023), 212 [arXiv:2304.05427].
  * (29) K. G. Chetyrkin and F. V. Tkachov, Integration by parts: The algorithm to calculate $\beta$-functions in 4 loops, Nucl. Phys. B 192 (1981), 159.
  * (30) S. Laporta, High-precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A 15 (2000), 5087 [arXiv:hep-ph/0102033].
  * (31) M. Nishimura and Y. Schröder, IBP methods at finite temperature, JHEP 1209 (2012) 051 [arXiv:1207.4042].
  * (32) A. I. Davydychev and J. B. Tausk, Two-loop self-energy diagrams with different masses and the momentum expansion, Nucl. Phys. B 397 (1993) 123.
  * (33) A. I. Davydychev and Y. Schröder, Recursion-free solution for two-loop vacuum integrals with “collinear” masses, JHEP 12 (2022) 047 [arXiv:2210.10593].
  * (34) A. A. Vladimirov, Method for computing renormalization group functions in dimensional renormalization scheme, Theor. Math. Phys. 43 (1980) 417 [Teor. Mat. Fiz. 43 (1980) 210].
  * (35) O. V. Tarasov, Generalized recurrence relations for two loop propagator integrals with arbitrary masses, Nucl. Phys. B 502 (1997) 455 [hep-ph/9703319].
  * (36) Y. Schröder, Loops for hot QCD, Nucl. Phys. Proc. Suppl. 183B (2008) 296 [arXiv:0807.0500].
  * (37) I. Ghisoiu and Y. Schröder, A new three-loop sum-integral of mass dimension two, JHEP 1209 (2012) 016 [arXiv:1207.6214].
  * (38) I. Ghisoiu and Y. Schröder, A new method for taming tensor sum-integrals, JHEP 1211 (2012) 010 [arXiv:1208.0284].
  * (39) M. Laine, P. Schicho and Y. Schröder, A QCD Debye mass in a broad temperature range, Phys. Rev. D 101 (2020) no.2, 023532 [arXiv:1911.09123].
  * (40) M. Laine, M. Meyer and G. Nardini, Thermal phase transition with full 2-loop effective potential, Nucl. Phys. B 920 (2017), 565 [arXiv:1702.07479].