Evaluation of multi-loop multi-scale Feynman integrals for precision physics

Modern particle physics is increasingly becoming a precision science that relies on advanced theoretical predictions for the analysis and interpretation of experimental results. The planned physics program at the LHC and future colliders will require three-loop electroweak and mixed electroweak-QCD corrections to single-particle production and decay processes and two-loop electroweak corrections to pair production processes, all of which are beyond the reach of existing analytical and numerical techniques in their current form. This article presents a new semi-numerical approach based on differential equations with boundary terms specified at Euclidean kinematic points. These Euclidean boundary terms can be computed numerically with high accuracy using sector decomposition or other numerical methods. They are then mapped to the physical kinematic configuration with a series solution of the differential equation system. The method is able to deliver 8 or more digits precision, and it has a built-in mechanism for checking the accuracy of the obtained results. Its efficacy is illustrated with examples for three-loop self-energy and vertex integrals and two-loop box integrals.