Numerical Multi-Loop Calculations via Finite Integrals and One-Mass EW-QCD Drell-Yan Master Integrals

TL;DR

This paper demonstrates that using finite integrals as a basis significantly improves the efficiency of numerical multi-loop Feynman integral evaluations, enabling calculations that were previously challenging or infeasible.

Contribution

It introduces a basis of finite integrals for numerical evaluation, leading to substantial performance improvements in multi-loop integral calculations.

Findings

Performance improvements with SecDec3 for two-loop integrals.

Order of magnitude speedup in three-loop form factor evaluations.

Numerical accessibility of complex topologies previously difficult to compute.

Abstract

We study a recently-proposed approach to the numerical evaluation of multi-loop Feynman integrals using available sector decomposition programs. As our main example, we consider the two-loop integrals for the $\alpha \alpha_s$ corrections to Drell-Yan lepton production with up to one massive vector boson in physical kinematics. As a reference, we evaluate these planar and non-planar integrals by the method of differential equations through to weight five. Choosing a basis of finite integrals for the numerical evaluation with SecDec3 leads to tremendous performance improvements and renders the otherwise problematic seven-line topologies numerically accessible. As another example, basis integrals for massless QCD three loop form factors are evaluated with FIESTA4. Here, employing a basis of finite integrals results in an overall speedup of more than an order of magnitude.

Tables

Table 1. Table 1: The four integral families used for the two-loop mixed EW-QCD Drell-Yan master integrals with either no massive internal lines or exactly one massive internal line.

Family $A$	Family $B$	Family $C$	Family $D$
$k_1^2$	$k_2^2$	$k_1^2$	$k_2^2$
$k_2^2$	${\left(k_2-p_3\right)}^2$	$k_2^2$	${\left(k_2-p_3\right)}^2$
${\left(k_1-k_2\right)}^2$	${\left(k_2-p_1-p_2\right)}^2$	${\left(k_1-k_2\right)}^2$	${\left(k_2-p_1-p_2\right)}^2$
${\left(k_1-p_1\right)}^2$	${\left(k_1-p_1\right)}^2$	${\left(k_1-p_1\right)}^2$	${\left(k_1-p_1\right)}^2$
${\left(k_2-p_1\right)}^2$	${\left(k_1-k_2\right)}^2$	${\left(k_2-p_1\right)}^2$	${\left(k_1-k_2\right)}^2$
${\left(k_1-p_1-p_2\right)}^2$	${\left(k_1-k_2+p_2\right)}^2$	${\left(k_1-p_1-p_2\right)}^2$	${\left(k_1-k_2+p_2\right)}^2$
${\left(k_2-p_1-p_2\right)}^2$	$k_1^2$	${\left(k_2-p_1-p_2\right)}^2-m^2$	$k_1^2-m^2$
${\left(k_1-p_3\right)}^2$	${\left(k_1-p_3\right)}^2$	${\left(k_1-p_3\right)}^2$	${\left(k_1-p_3\right)}^2-m^2$
${\left(k_2-p_3\right)}^2$	${\left(k_1-k_2+p_1+p_2\right)}^2$	${\left(k_2-p_3\right)}^2$	${\left(k_1-k_2+p_1+p_2\right)}^2$

Table 2. Table 2: Numerical performance of finite and conventional integral bases for two-loop mixed EW-QCD corrections to Drell-Yan lepton production in the physical region using SecDec 3 . For each integral in the above, the run time is given in seconds and the fractional difference from the analytical solution is given for the expansion coefficient which first gives rise to weight four multiple polylogarithms. In the final row of the table, total run times and worst-case relative accuracies are recorded for both integral bases.

$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}\left(s\right)$	$4s$	$3.17\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$4s$	$3.42\times {10}^{-3}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$8.39\times {10}^{-7}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$1.40\times {10}^{-6}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$11s$	$3.77\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$48s$	$2.11\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$1.33\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$12s$	$4.98\times {10}^{-5}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}\left(s\right)$	$5s$	$1.40\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$21s$	$3.31\times {10}^{-3}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$7.20\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$5s$	$1.11\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$13s$	$5.76\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$48s$	$4.56\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(t,s)$	$31s$	$3.31\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(t,s)$	$37s$	$2.74\times {10}^{-3}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}(s,t)$	$8s$	$6.70\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$8s$	$1.27\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$212s$	$3.44\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$17503s$	$3.59\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$197s$	$1.84\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$1337751s$	$.263$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}\left(s\right)$	$8s$	$4.18\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$1261s$	$1.63\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$6s$	$7.71\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$19s$	$6.84\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(t\right)$	$3s$	$1.03\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$5s$	$4.60\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$4s$	$4.60\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$5s$	$5.39\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$7s$	$5.06\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$370s$	$2.20\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$9s$	$1.89\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$143s$	$9.62\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$10s$	$2.85\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$82s$	$1.76\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$46s$	$1.21\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$101s$	$4.29\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$72s$	$4.90\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$332s$	$1.01\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(t\right)$	$6s$	$9.44\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$21s$	$5.18\times {10}^{-6}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$8s$	$2.25\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$1004s$	$4.82\times {10}^{-4}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$1.55\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$5s$	$2.29\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$2.28\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$3.86\times {10}^{-6}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$9s$	$2.88\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$441s$	$7.73\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$11s$	$3.66\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$10s$	$4.12\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$4.60\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$7s$	$3.45\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(t\right)$	$4s$	$1.38\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$6s$	$8.49\times {10}^{-6}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(t\right)$	$9s$	$2.76\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$16s$	$2.63\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$18s$	$1.04\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$41s$	$6.47\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$18s$	$1.23\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$12s$	$4.70\times {10}^{-4}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$9s$	$2.13\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(t\right)$	$3s$	$2.13\times {10}^{-5}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$31s$	$7.51\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$17s$	$7.51\times {10}^{-6}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$35s$	$3.28\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$19s$	$3.28\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$36s$	$3.94\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$113s$	$7.78\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$28s$	$8.37\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$101s$	$2.13\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$37s$	$5.84\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$26s$	$5.35\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$18s$	$2.92\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$319s$	$5.92\times {10}^{-2}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$50s$	$6.76\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$13s$	$9.49\times {10}^{-4}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}(s,t)$	$35s$	$7.64\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$20605s$	$9.87\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$1609s$	$4.39\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$564s$	$2.04\times {10}^{-2}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$202s$	$7.31\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$96s$	$2.35\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$201s$	$2.34\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$384s$	$8.12\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$150s$	$4.83\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$56538s$	$1.67\times {10}^{-2}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$280s$	$1.00\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$214135s$	$8.29\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}(s,t)$	$294s$	$1.21\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}(s,t)$	$3484378s$	$30.9$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$91s$	$3.76\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$87s$	$3.76\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$17s$	$5.15\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$20s$	$1.95\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}\left(s\right)$	$119s$	$2.32\times {10}^{-3}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}\left(s\right)$	$118s$	$2.12\times {10}^{-3}$
Total/Max:	$3995s$	$5.84\times {10}^{-3}$	Total/Max:	$5136862s$	$30.9$

Table 3. Table 3: Numerical performance of finite and conventional integral bases for massless three-loop form factors with FIESTA 4 . In addition to the run times, the fractional difference from the exact value is shown. In the final row of the table, total run times and average relative accuracies are recorded.

$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$6s$	$2.46\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$7s$	$3.84\times {10}^{-6}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$49s$	$8.88\times {10}^{-7}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$57s$	$1.36\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$35s$	$4.31\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$33s$	$2.68\times {10}^{-6}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$17s$	$4.09\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$15s$	$6.97\times {10}^{-6}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$16s$	$1.28\times {10}^{-4}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$22s$	$1.45\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$121s$	$5.42\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$92s$	$5.33\times {10}^{-6}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$93s$	$1.37\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$62s$	$1.96\times {10}^{-5}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$39s$	$2.25\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$87s$	$1.82\times {10}^{-5}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$34s$	$2.90\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$23s$	$4.55\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$30s$	$2.02\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$43s$	$5.05\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$19s$	$1.86\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$16s$	$1.42\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$69s$	$8.83\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$46s$	$1.57\times {10}^{-3}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$88s$	$1.93\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$117s$	$1.35\times {10}^{-3}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$74s$	$6.62\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$46s$	$1.40\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$26s$	$4.29\times {10}^{-7}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$17s$	$4.89\times {10}^{-6}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$20s$	$6.27\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$20s$	$6.27\times {10}^{-6}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$212s$	$1.50\times {10}^{-5}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$407s$	$5.66\times {10}^{-4}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$72s$	$3.11\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$75s$	$3.11\times {10}^{-6}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$74s$	$9.04\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$126s$	$2.18\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$128s$	$5.12\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$39094s$	$9.91\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$192s$	$2.68\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$19025s$	$9.38\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$127s$	$2.26\times {10}^{-6}$	$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$19586s$	$1.07\times {10}^{-4}$
Total/Average:	$1541s$	$1.76\times {10}^{-5}$	Total/Average:	$79016s$	$2.33\times {10}^{-4}$

Table 4. Table 4: Numerical performance of finite master integrals with FIESTA 4 for higher orders in the ϵ italic-ϵ \epsilon expansion. On the left, the weight six results of Table LABEL:tab:ff3L are reproduced. On the right, results for weight eight are shown.

$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$6s$	$2.46\times {10}^{-5}$	$9s$	$6.01\times {10}^{-5}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$49s$	$8.88\times {10}^{-7}$	$123s$	$7.72\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$35s$	$4.31\times {10}^{-5}$	$82s$	$9.55\times {10}^{-5}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$17s$	$4.09\times {10}^{-5}$	$46s$	$1.37\times {10}^{-4}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$16s$	$1.28\times {10}^{-4}$	$38s$	$3.03\times {10}^{-4}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$121s$	$5.42\times {10}^{-6}$	$313s$	$4.17\times {10}^{-5}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$93s$	$1.37\times {10}^{-5}$	$263s$	$4.47\times {10}^{-5}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$39s$	$2.25\times {10}^{-6}$	$103s$	$3.96\times {10}^{-6}$
$\stackrel{\left(10-2ϵ\right)}{\text{<img>}}$	$34s$	$2.90\times {10}^{-5}$	$90s$	$9.54\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$30s$	$2.02\times {10}^{-5}$	$80s$	$3.82\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$19s$	$1.86\times {10}^{-5}$	$45s$	$4.43\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$69s$	$8.83\times {10}^{-6}$	$202s$	$1.29\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$88s$	$1.93\times {10}^{-6}$	$274s$	$9.25\times {10}^{-6}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$74s$	$6.62\times {10}^{-6}$	$220s$	$4.79\times {10}^{-6}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$26s$	$4.29\times {10}^{-7}$	$82s$	$4.91\times {10}^{-6}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$20s$	$6.27\times {10}^{-6}$	$70s$	$1.72\times {10}^{-5}$
$\stackrel{\left(8-2ϵ\right)}{\text{<img>}}$	$212s$	$1.50\times {10}^{-5}$	$655s$	$2.36\times {10}^{-5}$
$\stackrel{\left(4-2ϵ\right)}{\text{<img>}}$	$72s$	$3.11\times {10}^{-6}$	$213s$	$1.57\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$74s$	$9.04\times {10}^{-6}$	$234s$	$1.29\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$128s$	$5.12\times {10}^{-6}$	$491s$	$2.22\times {10}^{-5}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$192s$	$2.68\times {10}^{-6}$	$761s$	$5.84\times {10}^{-6}$
$\stackrel{\left(6-2ϵ\right)}{\text{<img>}}$	$127s$	$2.26\times {10}^{-6}$	$485s$	$8.45\times {10}^{-6}$
Totals/Averages:	$1541s$	$1.76\times {10}^{-5}$	$4879s$	$4.90\times {10}^{-5}$

Equations

$$q(p_1)+\bar{q}(p_2)$$

$$q^{\prime }(p_1)+\bar{q}(p_2)$$

$$s=(p_1+p_2{)}^2,t=(p_1-p_3{)}^2,u=(p_1-p_4{)}^2\text{mbox}withs+t+u=0.$$

$$s>0,-s<t<0,-s<u<0$$

$$\displaystyle f_{1}^{\mathrm{A:38}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=54.2025pt,width=54.2025pt,keepaspectratio]{f1andf2}}(s)\quad f_{2}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{A}:38}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=54.2025pt,width=54.2025pt,keepaspectratio]{f1andf2}}(t)\quad f_{3}^{\mathrm{A:99}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f3andf4}}(s)$$

$$\displaystyle f_{4}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{A}:99}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f3andf4}}(t)\quad f_{5}^{\mathrm{A:53}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f5andf6}}(s)\quad f_{6}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{A}:53}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f5andf6}}(t)$$

$$\displaystyle f_{7}^{\mathrm{A:174}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f7}}(s,t)\quad f_{8}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{A}:174}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f8}}(t,s)\quad f_{9}^{\mathrm{A:182}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f9}}(s,t)$$

$$\displaystyle f_{10}^{\mathrm{A:247}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.22 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f10andf11}}(s,t)\quad f_{11}^{\mathrm{A:247}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.22 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f10andf11}}\left[(k_{2}-p_{3})^{2}\right](s,t)$$

$$\displaystyle f_{12}^{\mathrm{B:125}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f12}}(s)$$

$$\displaystyle f_{13}^{\mathrm{C:97}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=86.72267pt,width=86.72267pt,keepaspectratio]{f13andf14}}(s)\quad f_{14}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{C}:97}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=86.72267pt,width=86.72267pt,keepaspectratio]{f13andf14}}(t)\quad f_{15}^{\mathrm{C:76}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=21.68231pt,width=21.68231pt,keepaspectratio]{f15}}\quad$$

$$\displaystyle f_{16}^{\mathrm{C:69}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=54.2025pt,width=54.2025pt,keepaspectratio]{f16}}(s)\quad f_{17}^{\mathrm{C:69}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=54.2025pt,width=54.2025pt,keepaspectratio]{f17}}(s)\quad f_{18}^{\mathrm{C:99}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f18}}(s)$$

$$\displaystyle f_{19}^{\mathrm{C:102}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.6 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f19}}(s)\quad f_{20}^{\mathrm{C:78}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.1 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f20}}(s)\quad f_{21}^{\mathrm{C:212}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.1 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f21andf22}}(t)$$

$$\displaystyle f_{22}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{C}:212}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.1 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f21andf22}}(s)\quad f_{23}^{\mathrm{C:204}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.1 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f23andf25}}(t)\quad f_{24}^{\mathrm{C:204}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.1 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f24andf26}}(t)$$

$$\displaystyle f_{25}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{C}:204}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.1 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f23andf25}}(s)\quad f_{26}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{C}:204}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.1 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f24andf26}}(s)\quad f_{27}^{\mathrm{C:472}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f27}}(t)$$

$$\displaystyle f_{28}^{\mathrm{C:372}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f28}}(t)\quad f_{29}^{\mathrm{C:244}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f29andf30}}(t)\quad f_{30}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{C}:244}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f29andf30}}(s)\quad$$

$$\displaystyle f_{31}^{\mathrm{C:110}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f31}}(s)\quad f_{32}^{\mathrm{C:220}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f32andf33}}(t)\quad f_{33}^{\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{C}:220}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f32andf33}}(s)$$

$$\displaystyle f_{34}^{\mathrm{C:117}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f34}}(s)\quad f_{35}^{\mathrm{C:117}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f35}}(s)\quad f_{36}^{\mathrm{C:214}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f36}}(s,t)$$

$$\displaystyle f_{37}^{\mathrm{C:341}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f37}}(s,t)\quad f_{38}^{\mathrm{C:341}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f38}}(s,t)\quad f_{39}^{\mathrm{C:213}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f39}}(s,t)$$

$$\displaystyle f_{40}^{\mathrm{C:213}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.2 cm},height=58.53632pt,width=58.53632pt,keepaspectratio]{f40}}(s,t)\quad f_{41}^{\mathrm{C:374}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.45 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f41}}(s,t)\quad f_{42}^{\mathrm{C:246}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.45 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f42}}(s,t)$$

$$\displaystyle f_{43}^{\mathrm{C:245}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.45 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f43}}(s,t)\quad\quad f_{44}^{\mathrm{C:245}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.45 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f44}}(s,t)$$

$$\displaystyle f_{45}^{\mathrm{C:247}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.22 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f45andf46}}(s,t)\quad f_{46}^{\mathrm{C:247}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.22 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f45andf46}}\left[(k_{2}-p_{3})^{2}\right](s,t)$$

$$\displaystyle f_{47}^{\mathrm{D:341}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.3 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f47}}(s)\quad\quad f_{48}^{\mathrm{D:117}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.7 cm},height=75.88481pt,width=75.88481pt,keepaspectratio]{f48}}(s)$$

$$\displaystyle f_{49}^{\mathrm{D:125}}=\overset{(4-2\epsilon)}{\includegraphics[valign={m},raise={.22 cm},height=65.04034pt,width=65.04034pt,keepaspectratio]{f49}}\left[(k_{1}-k_{2}+p_{1}+p_{2})^{2}\right](s).$$

$$m_1={ϵ}^{2}s{f}_1^{\mathrm{A}:38}{m}_2={ϵ}^{2}t{f}_2^{\text{accentset}\mathrm{A}:38}{m}_3={ϵ}^2{s}^2{f}_3^{\mathrm{A}:99}{m}_4={ϵ}^2{t}^2{f}_4^{\text{accentset}\mathrm{A}:99}{m}_5={ϵ}^{3}s{f}_5^{\mathrm{A}:53}$$

$$m_6={ϵ}^{3}t{f}_6^{\text{accentset}\mathrm{A}:53}{m}_7={ϵ}^{3}st{f}_7^{\mathrm{A}:174}{m}_8={ϵ}^{3}st{f}_8^{\text{accentset}\mathrm{A}:174}{m}_9={ϵ}^4(s+t)f_9^{\mathrm{A}:182}$$

$$m_{10}={ϵ}^4{s}^{2}t{f}_{10}^{\mathrm{A}:247}{m}_{11}={ϵ}^4{s}^2{f}_{11}^{\mathrm{A}:247}{m}_{12}={ϵ}^4{s}^2{f}_{12}^{\mathrm{B}:125}{m}_{13}={ϵ}^{2}s{f}_1^{\mathrm{A}:38}$$

$$m_{14}={ϵ}^{2}t{f}_2^{\text{accentset}\mathrm{A}:38}{m}_{15}={ϵ}^{2}s{f}_{13}^{\mathrm{C}:97}{m}_{16}={ϵ}^{2}t{f}_{14}^{\text{accentset}\mathrm{C}:97}{m}_{17}=ϵ(1-ϵ)m^2{f}_{15}^{\mathrm{C}:76}$$

$$m_{18}={ϵ}^{2}s{f}_{16}^{\mathrm{C}:69}{m}_{19}=2{ϵ}^2(s-m^2)f_{16}^{\mathrm{C}:69}+{ϵ}^2(s-m^2)f_{17}^{\mathrm{C}:69}{m}_{20}={ϵ}^2{s}^2{f}_3^{\mathrm{A}:99}$$

$$m_{21}={ϵ}^2{t}^2{f}_4^{\text{accentset}\mathrm{A}:99}{m}_{22}={ϵ}^{3}s{f}_5^{\mathrm{A}:53}{m}_{23}={ϵ}^{3}t{f}_6^{\text{accentset}\mathrm{A}:53}{m}_{24}={ϵ}^2{s}^2{f}_{18}^{\mathrm{C}:99}$$

$$\displaystyle m_{25}=\frac{\epsilon^{2}s\Big{(}s(1-5\epsilon)-m^{2}(1-3\epsilon)\Big{)}}{4(1-3\epsilon)(s-m^{2})}f_{1}^{\mathrm{A:38}}-\frac{\epsilon(1-\epsilon)m^{4}}{2(s-m^{2})}f_{15}^{\mathrm{C:76}}+\frac{\epsilon^{2}(1-2\epsilon)^{2}}{s-m^{2}}f_{19}^{\mathrm{C:102}}$$

Full text

††preprint: MSUHEP-17-002, MITP/16-101, TCDMATH-17-03

1 Introduction

In the early days of hadron collider physics and Quantum Chromodynamics (QCD), the Drell-Yan process was recognized to be of fundamental importance Drell:1970wh and, over the years, it has received much attention. The electroweak corrections of relative order $\alpha$ were calculated in Baur:2001ze ; Dittmaier:2001ay ; Baur:2004ig . Another milestone of particular interest was set long ago when the corrections of relative order $\alpha_{s}^{2}$ were calculated in references Hamberg:1990np ; Harlander:2002wh ; Anastasiou:2003ds ; Melnikov:2006di . The fact that these corrections were found to be relatively large in size motivates one to consider the subleading terms in the two-loop perturbative expansion of the Standard Model Drell-Yan production cross section as well. The most important class of subleading two-loop contributions are those of relative order $\alpha\alpha_{s}$, what we shall hereafter refer to as mixed Electroweak-Quantum Chromodynamic (EW-QCD) corrections. The mixed Electroweak-Quantum Chromodynamic corrections of relative order $\alpha\alpha_{s}$ have already been studied in certain approximations. First, in reference Kilgore:2011pa , the virtual corrections without propagating $W$ or $Z$ bosons were calculated. More recently, the full set of two-loop corrections were considered in the single massive gauge boson resonance region Dittmaier:2014qza ; Dittmaier:2015rxo . Progress on real emission contributions and the infrared structure has been presented in Kilgore:2013uta ; deFlorian:2015ujt ; Bonciani:2016wya . Recently, the planar two-loop master integrals with up to two same-mass vector bosons in the loops were calculated in the Euclidean region Bonciani:2016ypc . The non-planar master integrals relevant to the problem are vertex integrals, which have been available in the literature for some time Gonsalves:1983nq ; Aglietti:2003yc .

In this paper, we consider the two-loop master integrals for mixed EW-QCD corrections to Drell-Yan production with up to a single massive vector boson exchanged and focus on physical kinematics. The integrals may be conveniently computed using the method of differential equations Kotikov:1990kg ; Kotikov:1991hm ; Kotikov:1991pm ; Bern:1992em ; Bern:1993kr ; Remiddi:1997ny ; Gehrmann:1999as . We employ a normal form basis Kotikov:2010gf ; Henn:2013pwa ; Henn:2013tua , integrate the $\epsilon$-expanded differential equations in terms of multiple polylgarithms LDanilevsky , and impose regularity conditions to fix the boundary constants. In some cases this is supplemented by explicit solutions worked out directly from Feynman parameters to all orders in the parameter of dimensional regularization, $\epsilon$. Analytical solutions for all of our $\epsilon$-expanded integrals are included with our arXiv submission through to weight four.111In this work, we make extensive use of the GiNaC-based implementation of the multiple polylogarithms Bauer:2000cp ; Vollinga:2004sn and therefore adopt the notation of reference Vollinga:2004sn for our function definitions. Although the weight four solution is one of the primary goals of our study, we found it interesting to calculate to one order higher in $\epsilon$ than necessary to facilitate checks in the second part of this work.

To cross-check an analytical solution for a multi-loop integral or to even produce the primary result for a phenomenological application, numerical evaluations using sector decomposition programs Binoth:2000ps ; Bogner:2007cr ; Borowka:2015mxa ; Smirnov:2015mct are very valuable. One might hope that one could just pick any integral basis, run one of the available programs for a physical point in phase space, and then obtain a reasonably accurate and precise result in, say, a few days’ time. However, for our mixed EW-QCD integrals we found this not to be the case in practice. Instead, picking a basis which is free of infrared and ultraviolet divergences allows us to obtain reliable and relatively fast results with SecDec 3 Borowka:2015mxa . This method of finite integrals introduced by Erik Panzer and the authors vonManteuffel:2014qoa ; Panzer:2014gra is generic and straightforward to automate. The approach has successfully been employed in the calculation of two-loop integrals for double Higgs production with exact top quark mass dependence Borowka:2016ehy ; Borowka:2016ypz already. Here, we present a detailed numerical study of performance gains stemming from the use of this technique. We quantify the effect also for a very different setup: the evaluation of the massless three-loop form factor integrals with FIESTA 4, which performs particularly well for such Euclidean integrals. Again, we observe drastic improvements with regard to the numerical convergence.

The outline of this article is as follows. In Section 2, we comment on the Feynman diagrams relevant to the mixed EW-QCD corrections and define the integral families which we use to describe our master integrals. In Section 3, we define our normal form basis, present the differential equations they satisfy, and solve the integrals in terms of multiple polylogarithms. In Section 4, we present the results of our numerical analysis of the two-loop mixed EW-QCD Drell-Yan integrals and discuss the crucial performance enhancements obtained by employing a basis of finite integrals. We present our conclusions in Section 5. Finally, in Appendix A, we present a self-contained numerical study of the master integrals for massless three-loop form factors up to contributions of weight eight.

2 Integral Families

We consider neutral-current and charged-current lepton pair production in quark-antiquark annihilation,

[TABLE]

where $p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=p_{4}^{2}=0$. As usual, the Mandelstam invariants are

[TABLE]

For the physical part of the phase space we have

[TABLE]

and selected Feynman diagrams are depicted in Figure 1 above.

We are interested in the two-loop master integrals required for the calculation of the $\mathcal{O}(\alpha_{s}\alpha)$ corrections to these processes. In this paper, we restrict ourselves to the subset of these integrals, which involve at most one massive vector boson, $W$ or $Z$, propagating in the loops. We will use the generic symbol $m$ for the vector boson mass. The relevant loop integrals can be indexed using the four integral families shown in Table 1. Integration by parts reduction Tkachov:1981wb ; Chetyrkin:1981qh ; Laporta:2001dd allows us to rewrite any Feynman integral in these families as a linear combination of master integrals. The reductions for this paper were performed with the implementation Reduze 2 vonManteuffel:2012np ; Studerus:2009ye ; Bauer:2000cp ; fermat and we adapt its conventions for the labeling of sectors, etc. An ideal choice for the master integrals will depend on the application and we will discuss suitable options for different applications in the following sections.

3 Analytical Solution from Differential Equations

We employ the method of differential equations to derive analytical solutions for the master integrals. Altogether, we find that we need to consider forty-nine master integrals for the system of differential equations. Some of these integrals are related by a permutation of external legs. They are covered by crossed integral families, which we will denote by an overline (e.g. we write $\mathrm{\accentset{\rule{2.79996pt}{0.8pt}}{F}}$:$x$ for the crossed version of sector $x$ from family $\mathrm{F}$). Note that, because our four-particle integral topologies are naturally a function of the invariants $s$ and $t$ (and never $u=-s-t$), we may without loss of intelligibility suppress all explicit momentum labels; instead, we use directed external lines and explicit function arguments to help clarify our basis integral definitions when necessary. As usual, we use dots to denote doubled propagators, heavy lines to denote massive propagators, and explicitly write all numerator insertions in square brackets. The starting point for the construction of our basis are the following integrals:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Building upon these definitions, we construct a convenient basis of integrals for the method of differential equations along the lines discussed in reference Gehrmann:2014bfa :222 We would like to point out further recent work on the analysis of master integrals and the construction of a normal form basis Lee:2013hzt ; Argeri:2014qva ; Lee:2014ioa ; Tancredi:2015pta ; Ablinger:2015tua ; Ita:2015tya ; Larsen:2015ped ; Primo:2016ebd ; Meyer:2016slj ; Georgoudis:2016wff ; Prausa:2017ltv ; Gituliar:2017vzm .

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that these definitions repeat nine massless integrals to allow for an independent treatment of the differential equations for the massless integrals, $m_{1},\ldots,m_{12}$, and the one-mass integrals, $m_{13},\ldots,m_{58}$. We employ the integration measure

[TABLE]

with $\epsilon=(4-d)/2$, which renders our integrals $m_{i}$ functions of the dimensionless variables

[TABLE]

only.

In this basis, we obtain a system of differential equations in normal form Kotikov:2010gf ; Henn:2013pwa ; Henn:2013tua for the vector $\vec{m}=(m_{i})$,

[TABLE]

with the letters $l_{k}(x,y)$ of the symbol alphabet

[TABLE]

and matrices $A^{(k)}$ of rational numbers.

Expanding the masters integrals $m_{i}$ about $\epsilon=0$, the differential equations fully decouple and can be integrated order-by-order in $\epsilon$ in terms of multiple polylogarithms. The alphabet is linear in both variables, which implies that the integration path can be chosen along the $x$ and $y$ axis, resulting in an expression written in terms of Goncharov $G$ functions with either $x$ or $y$ in the final argument. This representation can introduce spurious singularities and therefore may not be optimal for numerical evaluations Gehrmann:2015ora .

Instead, one can employ $\operatorname{Li}$ functions of fewer but possibly more involved arguments. Recently, it has been demonstrated explicitly that $\ln$, $\operatorname{Li}_{2}$, $\operatorname{Li}_{3}$, $\operatorname{Li}_{4}$, and $\operatorname{Li}_{2,2}$ functions are sufficient to reduce a general $G$ function of at most weight four Frellesvig:2016ske . In order to arrive at real-valued functions with a convergent power series representation in a specific region of phase space, it can be useful to allow for $\operatorname{Li}_{2,1}$ functions as well. A representation written in terms of such functions can be obtained either by recasting a solution written in terms of $G$ functions or by directly integrating the differential equations via an ansatz built out of the required functions. We perform these calculations using in-house routines based on the symbol and coproduct calculus Goncharov:2005sla ; Brown:2011ik ; Duhr:2011zq ; Duhr:2012fh .

Note that, for phenomenological applications, we are interested in the solutions of our integrals through to weight four. In the next section, however, we consider alternative bases of master integrals for the purpose of numerical evaluation which partially involve weight five functions at the $\epsilon$ order required for physics applications. In order to validate our numerical analysis of the solution, we therefore find it useful to obtain results through to weight five.

We construct an ansatz for our solutions in terms of the functions

[TABLE]

using the Duhr-Gangl-Rhodes algorithm Duhr:2011zq . We require the symbol of each function to not introduce letters beyond those already present in our differential equation (10) and we construct suitable power products of our letters (11) accordingly. For the arguments of the logarithms $\ln$ we choose

[TABLE]

As admissible arguments for the classical polylogarithms $\operatorname{Li}_{n}$, we obtain the following sixty-six power products of letters

[TABLE]

which have the property that not only the argument $z$ but also $1-z$ factorize over the original symbol alphabet. Extending this set by the argument $1$, forming pairs and selecting those pairs $(z_{1},z_{2})$, for which $1-z_{1}z_{2}$ factorizes over the alphabet, gives us 618 admissible arguments for $\operatorname{Li}_{2,1}$, $\operatorname{Li}_{2,2}$, and $\operatorname{Li}_{3,2}$ functions which we do not list here for the sake of brevity. Similarly, we find 3342 admissible triples $(z_{1},z_{2},z_{3})$ as arguments for our $\operatorname{Li}_{2,2,1}$ functions. While we have no proof that $\operatorname{Li}_{2,2,1}$ functions are strictly required, we were successful in constructing the genuine weight five part of the solution only when including them in addition to $\operatorname{Li}_{3,2}$ and $\operatorname{Li}_{5}$ functions. This observation is independent of further constraints to be described below, which motivate the inclusion of $\operatorname{Li}_{2,1}$ functions.

Not all of these functions are needed and we formulate further objectives for our functional basis. We prefer functions which are real-valued for physical kinematics

[TABLE]

This requirement is the reason why we include $\operatorname{Li}_{2,1}$ functions in our basis. In addition to requiring real-valuedness, we prefer functions which possess a convergent power series representation. Insisting upon such a representation automatically avoids additional manipulations which might otherwise be necessary to carry out numerical evaluations.

Indeed, we find that it is possible to use functions which are real-valued over the entire physical region of the phase space with a single important exception. The letter $l_{4}=1+y$ changes sign at $s=m^{2}$ and the functional form of our results will change in this region if we insist upon separating real and imaginary parts explicitly. This, however, is to be expected because there is a physical singularity at the point $s=m^{2}$, which would be regulated by the width of the massive vector boson inside the loop. Our solution, in fact, contains just a single function which is sensitive to the cut starting at $y=-1$:

[TABLE]

In contrast, no logarithms of the letter $l_{6}=x-y$ appear in our results and all other elements of our functional basis which depend of $l_{6}$ are real-valued for both $x>y$ and $x<y$.

Employing our functional basis, we construct an ansatz for the solution and match it against the differential equations (10). We employ regularity conditions and require real-valuedness in the Euclidean domain to fix the integration constants. In addition, we calculated $f_{1}^{\mathrm{A:38}}$, $f_{3}^{\mathrm{A:99}}$, $f_{5}^{\mathrm{A:53}}$, $f_{13}^{\mathrm{C:97}}$, $f_{15}^{\mathrm{C:76}}$, $f_{16}^{\mathrm{C:69}}$, $f_{18}^{\mathrm{C:99}}$, $f_{27}^{\mathrm{C:472}}$, and $f_{29}^{\mathrm{C:244}}$ explicitly from Feynman parameters to all orders in $\epsilon$. Most of these integrals are completely straightforward to evaluate. Some, such as

[TABLE]

turn out to be a bit more non-trivial.

As an additional sanity check, we found it convenient in some cases to look at asymptotic expansions of integrals in particular limits in order to explicitly confirm their scaling behavior. For this purpose, we found the Mathematica package asy.m Jantzen:2012mw to be extremely useful. For the convenience of the reader, we provide our analytical solutions through to weight four in the ancillary files sol-phys-AB.m and sol-phys-CD.m on arXiv.org.

4 Numerical Analysis

In an effort to check an analytical calculation such as the one presented in the previous section, it is desirable to numerically evaluate the loop integrals with an independent method. In other cases, a numerical evaluation might even be the method of choice, e.g. if it is not clear how to evaluate the integrals analytically. A standard method for numerical evaluation is sector decomposition and one might think that it is straightforward to numerically check the integrals $f_{i}$ in (6) using this approach. Unfortunately, this is not the case. Specifically, for double box integrals with additional numerators or dotted propagators we often encountered cases where none of the available software packages is able to provide reliable results. In some cases, issues related to the presence of problematic singularity structures in the integrand can be tracked down and, with dedicated effort, cured Borowka:2013cma . For what concerns the two-loop mixed EW-QCD topologies considered in this paper, the most challenging cases for the code are seven-line integrals which have $\epsilon^{-4}$ poles and doubled propagator denominators.

As discussed in earlier work by the authors and Erik Panzer vonManteuffel:2014qoa ; vonManteuffel:2015gxa , it is possible to employ a basis of finite Feynman integrals defined by allowing for higher-dimensional integrals with potentially higher powers of the propagator denominators to extract all poles in $\epsilon$ analytically before performing any numerical integrations.333Let us remind the reader that, due to the principle of analytical continuation, the existence of an on-shell Euclidean region for all integrals guarantees that our procedure applies equally well to all integrals when physical kinematics is considered. The change of basis is performed with integration by parts reductions, whose computational complexity is non-negligible but still lower than that of the reductions for the amplitude. This is true both for the processes discussed in this work and for all other phenomenologically relevant problems which we have studied.

The finite integrals could, in principle, be evaluated by direct numerical integration. However, their integrands may still contain structures which spoil the convergence of the numerical integrations, such as integrable but large variations near the boundaries. In our experiments with a limited number of Euclidean integrals, we find that an additional sector decomposition step improves the performance of the evaluation. Furthermore, handling branch cuts requires a dedicated treatment for physical kinematics. We therefore find it convenient to employ existing sector decomposition programs for the numerical integration of our finite integrals.

We employ the sector decomposition program SecDec 3 to evaluate our two-loop topologies using either conventional or finite master integrals. As is apparent from Table LABEL:tab:mixedDYnums, working with a basis of finite integrals improves the numerical performance tremendously. Let us now explain in more detail how we used the code to produce Table LABEL:tab:mixedDYnums. Although we made a good effort to use the program in an appropriate way, it is certainly possible that other researchers might manage to produce better results in less time. However, we at least took care to use identical program settings and hardware while carrying out our numerical studies of the conventional and finite integral bases.444We followed the recommendations given in Borowka:2015mxa as much as possible, using, for example, their new sector decomposition strategy G2 and Mathematica 9.0.1 for all of our calculations. Using this version of Mathematica avoided some parallelization issues encountered with versions 8 or 10 in our experiments.

Altogether, we ran for roughly sixty days on the four cores of a i7-4940MX processor. For our finite integral evaluations, just over an hour was required; the excessive run time is due to the fact that, in spite of allowing the VEGAS routine Lepage:1977sw provided by CUBA Hahn:2004fe to sample up to $5\times 10^{8}$ points, numerous conventional basis integrals failed to achieve the very modest goal of $\mathtt{epsrel}=1\times 10^{-2}$ and $\mathtt{epsabs}=1\times 10^{-4}$ at the physical phase space point $p=\{s\rightarrow 17,t\rightarrow-7,m^{2}\rightarrow 6241/1681\}$. Analytical expressions for the finite basis integrals were obtained by using the explicit connection between the finite integral basis and the normal form integral basis. We found good agreement in all cases. At worst, we recorded a fractional difference of $6\times 10^{-3}$ at weight four between our exact results evaluated at $p$ and the numerical results output by the program. The average relative accuracy at weight four was found to be $9\times 10^{-4}$, quite acceptable given our modest accuracy goal.

Due to the fact that there is a substantial loss of precision in rotating from the finite integral basis to the normal form integral basis, we found that our initial $\mathtt{epsrel}=1\times 10^{-2}$ and $\mathtt{epsabs}=1\times 10^{-4}$ run of SecDec 3 did not allow us to directly check the normal form basis to similar part per mille precision. To achieve this, we had to carry out a second, higher precision run of the program with $\mathtt{epsrel}=1\times 10^{-4}$ and $\mathtt{epsabs}=1\times 10^{-8}$, which took roughly eleven days of run time. Let us stress again that it would be hopeless to attempt this exercise using the conventional integral basis.

We find that, in several other cases, significantly more mileage can be squeezed out of publicly available sector decomposition programs simply by working with a well-behaved, finite integral basis. As a second example, we present the evaluation of three-loop form factors in massless QCD in Appendix A. We employ the other publicly available program under active development, FIESTA 4, and observe dramatic gains in both speed and numerical convergence when using a basis of finite integrals.

5 Conclusions and Outlook

In this paper, we considered the two-loop integrals relevant to the $\alpha\alpha_{s}$ corrections to Drell-Yan production with up to a single massive vector boson exchanged. As a reference, we calculated their Laurent expansion through to weight five in terms of multiple polylogarithms using the method of differential equations. Our representation in terms of real-valued functions allows for fast and precise numerical evaluations over the entire physical region of the phase space. We found it challenging to even check our analytical solutions using available sector decomposition programs. Employing a basis of finite integrals systematically improved the situation and rendered all integrals numerically accessible with SecDec 3, both in Euclidean and physical kinematics. Order of magnitude improvements both in program run time and integration error were also found for massless three-loop form factor integrals using FIESTA 4 when using finite instead of conventional master integrals. For the finite integrals, we allowed for shifts to higher numbers of spacetime dimensions and additional powers of denominators; the actual change of basis is performed with integration by parts reductions in a highly-automated way. An implementation of the algorithm to construct a basis of finite integrals is publicly available in the package Reduze 2.1 on HepForge. Numerical evaluations along the lines discussed in this paper are not only interesting for cross-checks of analytical solutions, but may actually serve as the primary method for phenomenological applications in especially complicated cases Borowka:2016ehy ; Borowka:2016ypz .

Acknowledgments

We would like to thank Gudrun Heinrich, Johannes Schlenk, and Alexander Smirnov for help with SecDec and FIESTA. We would also like to thank Erik Panzer for his careful reading of our manuscript. The work of RMS was supported in part by the European Research Council through grants 291144 (EFT4LHC) and 647356 (CutLoops). We are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and support. Our figures were generated using Jaxodraw Binosi:2003yf , based on AxoDraw Vermaseren:1994je .

Appendix A Numerical Analysis for Massless Three-Loop Form Factors

In this appendix, we present the results of our FIESTA 4-based numerical study of the finite and conventional integral bases for the massless three-loop vertex functions (i.e. the bases employed in references vonManteuffel:2015gxa and Gehrmann:2010ue respectively). They show the advantages of using a basis of finite integrals in a setup independent of what was employed in Section 4. Due to the fact that the three-loop form factor master integrals have been calculated to weight eight Gehrmann:2005pd ; Lee:2010ug ; Lee:2010ik ; Henn:2013nsa ; vonManteuffel:2015gxa , we can also study how well our finite basis performs for the evaluation of higher order coefficients in the $\epsilon$ expansion.

By experimenting with the code and discussing some aspects of it with the developer, we settled on the FIESTA 4 settings $\mathtt{ComplexMode=False}$, $\mathtt{STRATEGY=STRATEGY\_X}$, $\mathtt{BucketSize=28}$, $\mathtt{NegativeTermsHandling=None}$, and $\mathtt{CurrentIntegratorSettings=}$

$\mathtt{\{\{}$“$\mathtt{mineval}$”,“$\mathtt{1000}$”$\mathtt{\},\{}$“$\mathtt{maxeval}$”,“$\mathtt{500000}$”$\mathtt{\},\{}$“$\mathtt{nstart}$”,“$\mathtt{50000}$”$\mathtt{\}\}}$ for the trivial phase space point $q^{2}=-1$. Especially crucial is the setting $\mathtt{NegativeTermsHandling=None}$, which seems to substantially improve the performance for Euclidean Feynman integrals in general. We also configured the code to run on all four cores of a i7-4940MX processor. Once again, we used the VEGAS Monte Carlo integration routine provided by CUBA, but, this time, we found that Mathematica 8.0.4 is a better choice than Mathematica 9.0.1 for the front end.

Let us begin by discussing Table LABEL:tab:ff3L, where all twenty-two irreducible integral topologies in both the finite and conventional integral bases are evaluated through to the terms of weight six. Reasonable results are obtained in both integral bases using the program settings described above. However, in going from the conventional basis to our chosen finite basis, the program run time drops by more than a factor of fifty: from just shy of twenty-two hours to about half an hour. There is, however, another crucial difference between the two sets of integrals which plays an even more important role when one goes from three to four loops: for the fixed program settings given above, one finds that the numerical convergence of the VEGAS algorithm is far superior for the finite basis integrals.555Let us point out that, in recent work on massless form factors vonManteuffel:2015gxa , we were able to numerically evaluate a finite non-planar twelve-line master integral through to weight eight to four decimal digit accuracy using FIESTA 4 in about one day on a desktop computer.

In fact, the average relative error of the weight six contributions with respect to the exact results is $2\times 10^{-5}$ for the finite integral basis compared to $2\times 10^{-4}$ for the conventional integral basis, an improvement of an order of magnitude. Finally, we see from Table LABEL:tab:higherff3L that, in going from weight six to weight eight with finite integrals, the run time increases by only a factor of three (with a slight loss of precision for fixed program settings).666This shows that, to obtain acceptably fast and stable numerical evaluations of master integrals at higher orders in $\epsilon$, our integral basis is an effective alternative to the basis suggested in Chetyrkin:2006dh , where master integrals were systematically chosen to have the property that their coefficients remain finite in the $\epsilon\rightarrow 0$ limit. Given the long run times with conventional integrals already at weight six, we refrained from trying to perform a similar weight eight evaluation of them.