Two-loop master integrals for the leading QCD corrections to the Higgs coupling to a $W$ pair and to the triple gauge couplings $ZWW$ and $\gamma^*WW$

We compute the two-loop master integrals required for the leading QCD corrections to the interaction vertex of a massive neutral boson $X^0$, e.g. $H,Z$ or $\gamma^{*}$, with a pair of $W$ bosons, mediated by a $SU(2)_L$ quark doublet composed of one massive and one massless flavor. All the external legs are allowed to have arbitrary invariant masses. The Magnus exponential is employed to identify a set of master integrals that, around $d=4$ space-time dimensions, obey a canonical system of differential equations. The canonical master integrals are given as a Taylor series in $\epsilon = (4-d)/2$, up to order four, with coefficients written as combination of Goncharov polylogarithms, respectively up to weight four. In the context of the Standard Model, our results are relevant for the mixed EW-QCD corrections to the Higgs decay to a $W$ pair, as well as to the production channels obtained by crossing, and to the triple gauge boson vertices $ZWW$ and $\gamma^*WW$.


Introduction
The calculation of Feynman integrals with massive particles both as external legs and as internal lines is not an easy task. Indeed, the combination of external kinematic invariants and internal masses may give rise to physical and non-physical singularities which require the use of special functions, with non-trivial arguments, embedding them. Identifying the rise of such functions within the parametric representation of Feynman integrals is very challenging. Therefore, rather than by direct integration, multivariate Feynman integrals may be more simply determined by solving differential equations (DEs) [1][2][3].
In general, Feynman integrals in dimensional regularization obey relations that can be used, on the one side, to identify a basis of independent integrals, dubbed master integrals (MIs), and, on the other side, to write special equations satisfied by the MIs themselves. MIs are found to obey systems of linear, first-order partial DEs in the kinematic variables. Solving these equations, provided that their values or behavior at special points is known, becomes a method to completely determine MIs, hence to compute Feynman integrals alternatively to their direct integration (as reviewed in [4,5]).
The entries of the matrix associated to the system of DEs depend, in general, on the kinematic invariants and on the space-time dimension d. Although it is a mathematically interesting problem, finding the expression of the MIs for arbitrary values of d is not always possible and, according to the physical context, it may be sufficient to know the MIs around a critical dimension d c , with d = d c + and → 0. In a perturbative approach, the solution of the system around d = d c may admit a representation in terms of iterated integrals (as reviewed in [6,7]), where the matrix associated to the system constitutes the integration kernel. Therefore, the structure of such a matrix has a direct impact on the form of the solutions, namely on the functions required to classify them: simplifying the matrix means simplifying the solutions.
The idea of finding MIs that obey canonical systems of DEs, i.e. systems with an associated matrix where the dependence on the space-time dimensions is decoupled from the kinematics [8,9], has led to a substantial improvement of the system-solving strategy [10][11][12][13][14][15][16][17][18][19], and to the availability of many novel results. In the case of Feynman integrals that depend on several scales, we have shown that the Magnus exponential [11] is an efficient tool to derive MIs obeying canonical systems starting from a basis of MIs that obey systems of DEs whose matrix has a linear dependence on the space-time dimension [20][21][22].
In the Standard Model, the coupling between two W bosons and one neutral boson X 0 = H, Z, γ * is present in the tree-level Lagrangian 1 . At one-loop, the X 0 W + W − interaction receives electro-weak (EW) corrections, either via bosonic-or via fermionic-loop. Strong (QCD) corrections must proceed through a closed quark-loop so that they can first occur at the two-loop level. In this article, we present the calculation of the two-loop three-point integrals required for the determination of the leading QCD corrections to the interaction vertex between a neutral boson X 0 with arbitrary mass and a pair of W bosons of arbitrary squared four-momenta (X 0 W + W − ), mediated by a fermion loop of a SU (2) L quark doublet, with one massive and one massless flavors. In what follows, we refer to the massive flavor as to the top (m t = m), and to the massless one as to the bottom (m b = 0). Representative Feynman graphs for the considered integrals are shown in figure 1. The MIs for the case in which only massless quarks propagate in the loops has been previously studied in [23][24][25].
Our results represent the full set of MIs needed to compute the O(αα s ) corrections to the Higgs decay into a pair of W bosons, and to the triple gauge boson processes Z * W W and γ * W W , with leptonic final states, at e + e − colliders. As for the latter process with 1 Several motivated extensions of the Standard Model feature an extended Higgs sector with Yukawa couplings to the SU (2) fermion doublets. In particular, one or more neutral pseudoscalar bosons might be part of the spectrum, together with other neutral scalar bosons. While we do not refer explicitly to this possibility, our results would also be applicable to the case X 0 = S 0 , A 0 , where we schematically indicate with S 0 the scalars and with A 0 the pseudoscalars. semi-leptonic or hadronic final states, our MIs would only be a subset of the needed MIs. They are also a subset of the MIs needed for the computation of the two-loop mixed EW-QCD corrections to the Higgs production cross section either in the W W -fusion channel or in association with a W boson, and to W W production in higher multiplicity processes. The same MIs would also be needed for the computation of a class of diagrams entering the NNLO EW corrections. Except for the first-generation quarks (that are approximately degenerate), the fermionic one-loop diagrams always involve an SU (2) L doublet with a (nearly or exactly) massless flavor. Is is then clear that the corrections due to photon exchange between the fermionic lines share the same topologies as the ones of the leading QCD corrections.
We distinguish two sets of integrals, according to the flavor that couples to the X 0 boson, i.e. either the massive or the massless one, for which integration-by-parts reduction returns 24 and 23 MIs, respectively. The calculation of the MIs proceeds according to the following strategy. We identify a set of MIs that obey systems of DEs whose matrix has a linear dependence on d = 4 − 2 , and, by means of the Magnus exponential, we derive a canonical set of master integrals. The matrices associated to the canonical systems admit a logarithmic-differential form (d log) with rational arguments, therefore, the canonical MIs can be cast in Taylor series around d = 4 with coefficients written as combinations of Goncharov polylogarithms (GPLs) [26][27][28][29]. Boundary conditions are imposed by requiring the regularity of the solutions at special kinematics points, and by using simpler integrals as independent input. The analytic expressions of the MIs have been numerically evaluated with the help of GiNaC [30] and successfully tested against the values provided by the public computer code SecDec [31]. The package Reduze [32] has been used throughout the calculations.
The paper is organized as follows. In section 2 we fix our notation and conventions. In section 3 we discuss the general features of the systems of DEs satisfied by the master integrals and its general solution in terms of iterated integrals. In section 4 we describe the computation of the two-loop MIs for X 0 W + W − in Euclidean kinematics and in section 5 we discuss the analytic continuation of our result. Conclusions are given in section 6. In appendix A we recall the main properties of iterated integrals. In appendix B we list the coefficients of the linear combinations of MIs that satisfy a canonical system of DEs and finally, in appendix C, we give the expressions of the d log-form of the matrices associated to such systems. The analytic expressions of the canonical MIs up to O( 5 ) are attached to the arXiv version of the manuscript as ancillary files.

Notation and conventions
In this paper we will consider the two-loop three-point functions of a X 0 boson with momentum q, and two W bosons with momenta p 1 and p 2 , where s = q 2 = (p 1 + p 2 ) 2 and p 2 1 = p 2 2 = 0 . Figure 1: Representative two-loop Feynman diagrams contributing to the X 0 W + W − interaction, where X 0 = H, Z, γ * . Similar diagrams where t and b quarks are exchanged are also taken into account. The diagrams have been generated using FeynArts [33].
The calculation involves the evaluation of Feynman integrals in d = 4 − 2 dimensions of the type In our conventions, the integration measure is defined as where m 2 is the mass of the top quark circulating in the loops, µ the 't Hooft scale of dimensional regularization and S = (4π) Γ(1 + ) . (2.5)

System of differential equations for master integrals
In this section we briefly discuss the general features of the systems of DEs obeyed by the MIs and the properties of the corresponding solutions. The details of the calculations of the MIs for X 0 W + W − are described in section 4. The two-loop Feynman diagrams contributing to X 0 W + W − can be reduced to four parent topologies, which are depicted in figure 2. The integrals belonging to these topologies depend on the three external invariants as well as on the top mass m 2 . These four dimensionful parameters can be combined into three independent dimensionless variables, x = (u, z,z) for topologies (a)-(b) and x = (v, z,z) for topologies (c)-(d), whose explicit definition will be later specified. The MIs satisfy a linear system of partial DEs in these variables, which, if we organize the MIs into a vector F, can be combined into a matrix equation for the total differential of F, In general, the matrix-valued differential form K = K a dx a (a = 1, 2, 3) depends both on the kinematic variables and on the spacetime dimension d = 4 − 2 . Since the left-hand Figure 2: Two-loop topologies for X 0 W + W − interactions. Thin lines represent massless propagators and thick lines stand for massive ones. The dashed external line corresponds to the off-shell leg with squared momentum equal to s whereas the red and blue lines represent the two external vector bosons with off-shell momenta p 2 1 and p 2 2 respectively. side of the system in eq. (3.2) is a total differential by construction, it is easy to show that K satisfies the (matrix) integrability condition Starting from a basis of MIs associated to a matrix K with a linear dependence on , one can use the Magnus exponential [11,20] and apply the procedure outlined in section 2 of [21] in order to perform a basis transformation and obtain a canonical set of MIs [8] I enjoying two remarkable features: first, the canonical basis I obeys a system of DEs where the dependence on is factorized from the kinematics and, in addition, the kinematic matrices can be organized into a logarithmic differential form, referred to as canonical d log-form. Thus, the canonical basis I satisfies a system of equations of the form is the d log matrix written in terms of the differentials d log η i , whose arguments η i = η i ( x) solely enclose the kinematic dependence and form the so called alphabet of the problem. The coefficient matrices M i have rational-number entries. Due to the -factorization, the integrability condition for eq. (3.4) splits into

General solution
At any point x, the general solution of the canonical system of DEs (3.4) can be expressed in terms of Chen's iterated integrals [34] as the path-ordered exponential where I( , x 0 ) is a constant vector depending on only and γ is a piecewise-smooth path connecting x 0 to x, (3.8) In the limit x → x 0 the integration path γ shrinks to a point and I( , x) → I( , x 0 ). In this perspective, the integration constants I( , x 0 ) which, together with dA, completely specify the solution, can be thought as the initial values of the MIs, which then evolve to arbitrary points x under the action of the path-ordered exponential. It can be proven that, whenever γ does not cross any singularity or branch cuts of dA (but at its endpoints), the path-ordered exponential is independent of the explicit choice of the path. By choosing a proper normalization, we can assume all canonical MIs to be finite in the → 0 limit, in such a way that I( x) admits a Taylor expansion in , and, according to eq. (3.7), the n-th order coefficient is given by where we introduced the weight-k operator which iterates k ordered integrations of the matrix-valued 1-form dA along the path γ. According to eq. (3.5), each entry of ∆ (k) is a linear combination of Chen's iterated integrals of the type It should be remarked that, as explicitly indicated in (3.12), individual Chen's iterated integrals are, in general, functionals of the path and that only the full combinations appearing as entries of ∆ (k) are independent of the particular choice of γ. The most relevant properties of Chen's iterated integrals are summarized in appendix A.
4 Two-loop master integrals for X 0 W + W − In this section we present the solution of the system of DEs for the MIs associated to the four topologies (a)-(d). Since topologies (a)-(b) and (c)-(d) belong to two distinct integral families, we discuss their evaluation separately.

Topologies (a)-(b)
The two topologies (a) and (b) belong to the 7-denominator integral family identified by the set of denominators where k 1 and k 2 are the two loop momenta. The integrals belonging to this family can be reduced to a set of 29 MIs which are conveniently expressed in terms of the dimensionless variables u, z andz defined by The same parametrization for p 2 1 and p 2 2 was used also for the massless triangles considered in [25]. The following set of MIs obeys a system of DEs which is linear in : where the T i are depicted in figure 3. We observe that some of integrals T i are trivially related by p 2 1 ↔ p 2 2 symmetry, so that the actual number of independent integrals is reduced to 23. However, in order to determine the solution of the DEs with the method described in section 3, i.e. by simultaneously integrating the whole system of equations, one has to consider the full set of integrals given in eq. (4.3).
where λ is the Källén function related to the external kinematics, Explicit expressions for the coefficients c i, j are given in appendix B.1. The alphabet of the corresponding d log-form contains the following 10 letters: The coefficient matrices M i are collected in the appendix C.1. It can be easily checked that all letters are real and positive in the region .
If one fixes m 2 > 0, this corresponds to a patch of the Euclidean region, s , p 2 1 , p 2 2 < 0, defined by the following constraints Since the alphabet is rational, the solution can be directly expressed in terms of GPLs with argument depending on the kinematics variables u, z andz. The prescriptions for the analytic continuation to the other patches of the Euclidean region (s, p 2 1 , p 2 2 < 0) and to the physical regions are given in section 5.
Imposing the regularity of our solutions at the unphysical thresholds, z,z = 0 (corresponding to p 2 1 = 0) and z,z = 1 (corresponding to p 2 2 = 0) entails relations between the boundary constants. These relations allow us to derive all boundary constants from five simpler integrals I 1,3,6,7,15 , which are obtained in the following way: • I 1 is a constant to be determined by direct integration and, due to the normalization of the integration measure (2.4), it is simply set to I 1 ( , x) = 1. (4.10) • I 3 can be obtained by direct integration • Besides being regular in the massless kinematic limit z → 1 (p 2 2 → 0), I 6 is reduced, through IBPs, to a two-loop vacuum diagram, Therefore, by using as an input the analytic expression of the two-loop vacuum graph, we can fix the boundary constants by matching the z → 1 limit of the expression of I 6 obtained from the solution of the DE against the -expansion of eq. (4.12), (4.14) • I 7 and I 15 can be directly integrated The results have been numerically checked, in both the Euclidean and the physical regions, with the help of the public computer codes GiNaC and SecDec 3.0, and their analytic expressions are given in electronic form in the ancillary files attached to the arXiv version of the manuscript.

Topologies (c)-(d)
The topologies (c) and (d) belong to the 7-denominator family defined by the set of denominators where k 1 and k 2 are the two loop momenta. The integrals belonging to this integral family can be reduced to a set of 31 MIs which are conveniently expressed in terms of the variables v, z andz, defined by The following set of MIs obeys a system of DEs which is linear in : where the T i are depicted in figure 4. As for the case of topologies (a)-(b), some of integrals T i are related by p 2 1 ↔ p 2 2 , so that the total number of independent integrals is 24. However, as discussed already after eq. (4.4), we work with the complete set of integrals given in eq. (4.19). The Magnus exponential allows us to obtain a set of canonical MIs obeying a system of equations of the form (3.4),  T1   T2  T3  T4  T5  T6   T7  T8  T9  T10  T11  T12   T13  T14  T15  T16  T17  T18   T19  T20  T21  T22  T23  T24   T25  T26  T27  T28  T29  T30 T31 Figure 4: Two-loop MIs T 1,...,31 for the topologies (c)-(d). Graphical conventions are the same as in figure 2. Dots indicate squared propagators.
where ρ ≡ √ −s √ 4m 2 − s and λ is defined as in eq. (4.6). The expression of the coefficients c i, j is given in appendix B.2. The alphabet of the corresponding d log-form contains the following 18 letters The coefficient matrices M i are collected in the appendix C.2. In this case, all the letters are real and positive in the region If one fixes m 2 > 0, this corresponds to a patch of the Euclidean region, s , p 2 1 , p 2 2 < 0, defined by the following constraint The solution of the system of DEs is straightforwardly obtained in terms of Chen's iterated integrals. Moreover, since the alphabet is rational, the solution can be converted in terms of GPLs of argument 1, with kinematic-dependent weights, as we discuss in appendix A. The prescriptions for the analytic continuation to the other patches of the Euclidean region and to the physical regions are given in section 5.
The boundary constants can be fixed by demanding the regularity of the basis (4.19) for vanishing external momenta, s = p 2 1 = p 2 2 = 0. In particular, if we choose as a base-point for the integration The boundary constants of integrals I 1, 5, 10 can be taken from the previous topologies in equations (4.10,4.14), whereas for integrals I 19, 21, 24 the boundary constants are fixed as follows: • The boundary constants for I 19 and I 24 can be determined by imposing regularity at the pseudothresholds v → 1 (s = p 2 1 = p 2 2 = 0) and, respectively, z → 1,z → 1 (both corresponding to p 2 2 = 0), • Finally, the boundary constants for I 21 can be fixed by observing that, from (4.21), we can derive Therefore, in order for F 21 ( , x 0 ) to be regular we must demand All results have been numerically checked, in both the Euclidean and the physical regions, with the help of the computer codes GiNaC and SecDec 3.0, and the analytic expressions of the MIs are given in electronic form in the ancillary files attached to the arXiv version of the manuscript.

Change of variables and analytic continuation
In this section we discuss in detail the variables used to parametrize the dependence of the MIs on the kinematic invariants. In particular, we elaborate on the prescriptions to analytically continue our results to arbitrary values of s, p 2 1 , p 2 2 . Both topologies (a)-(b) and (c)-(d) feature two independent, kinematic structures: 1. the off-shell external legs are responsible for the presence in the DEs of the square root of the Källén function, λ(s, p 2 1 , p 2 2 ); 2. the presence of massive internal lines can generate square roots in the DEs, as in the case of topologies (c)-(d) where one has also √ −s √ 4m 2 − s.
In the following we separately discuss the variable changes that rationalize the two types of square roots.

Off-shell external legs: the z,z variables
To deal with the square root of the Källén function, we begin by choosing one of the external legs as reference, s, and trading the other squared momenta for dimensionless ratios In the (s, τ 1 , τ 2 ) variables, the square root of the Källén function is proportional to and is rationalized by the following change of variables [25] τ 1 = zz , (5.3) Without loss of generality, we choose the following root of eq. (5.4) Varying the pair (τ 1 , τ 2 ) in the real plane, we identify the following possibilities for z,z where the first five regions were discussed also in [25]. A graphical representation of these eight regions in the ( The variables z,z are complex conjugates in region I, where λ(1, τ 1 , τ 2 ) < 0, and real in all the other regions. In regions I-V one has τ 1,2 > 0, which requires that either s, p 2 1 , p 2 2 < 0 or s, p 2 1 , p 2 2 > 0. The former case defines the Euclidean region. The latter case, for λ(1, τ 1 , τ 2 ) > 0, describes 1 → 2 or 2 → 1 processes involving three timelike particles. Region V is where λ(1, τ 1 , τ 2 ) = 0, so that z =z. Since our expressions are obtained in general for z =z, the limitz → z has to be taken carefully. Regions VI-VIII have at least one of the τ i < 0, which requires either two external legs to be spacelike and the remaining Figure 5: Regions of the (τ 1 , τ 2 )-plane classified in eq (5.8). Region V, which is identified by the condition λ(1, τ 1 , τ 2 ) = 0, corresponds to the blue curve.
one to be timelike, or vice versa. The former configuration, in the 2 → 1 kinematics, describes the vertex entering the production of a timelike particle via the "fusion" of two spacelike particles. In regions other than II, the variables z,z are not in the half of the unit square where all the letters are real, therefore analytic continuation is required. A consistent physical prescription is inherited in regions VI-VIII from the Feynman prescription on the kinematic invariants, and it is naturally extended to the other regions, as we argue below. For the moment we hold s < 0, and we will discuss later the case s > 0. In region VI, s < 0 and p 2 1 , p 2 2 > 0, then so that the vanishing imaginary parts outside the square root in eq (5.7) cancel against each other, and only the one stemming from the square root is left: In region VII, p 2 1 > 0 and s, p 2 2 < 0, then where the approximate equalities are allowed because the factor in the bracket is always positive, and a redefinition of ε is understood. In region VIII, p 2 2 > 0 and s, p 2 1 < 0, then where again the approximate equalities are allowed because the factor in the bracket is always positive, and a redefinition of ε is understood.
We have so far only discussed the case in which s < 0. It is easy to see that, if instead s > 0, the prescription on z,z is the opposite.
In regions I-V there is no physical prescription for the analytic continuation of z,z. Indeed, if s, p 2 1 , p 2 2 > 0, then the vanishing imaginary parts of the Feynman prescription cancel out in the ratios τ 1 and τ 2 : This cancellation affects also region I, where λ(1, τ 1 , τ 2 ) < 0 and z * =z. Indeed, while this condition fixes the relative sign of their imaginary parts, the sign of Im z depends on the choice of the branch of the square root in eq. (5.7), which is not fixed. This last statement holds true also in the Euclidean region. This ambiguity is resolved by the definite iε prescription in regions VI-VIII discussed above. In order to have a smooth analytic continuation in the Euclidean, in region I, one chooses the branch of the square root that gives Im λ(1, τ 1 , τ 2 ) > 0, and in regions III-IV, one assigns vanishing imaginary parts for z,z according to the previous discussion. The opposite prescription should be used if the three external legs are timelike.
Summarizing, according to the sign of s, we choose the following analytic continuation prescriptions for z,z in the whole real (p 2 1 , p 2 2 ) plane

Internal massive lines: the u, v variables
For topologies (a)-(b), the change of variables eq. (5.4) is actually enough to rationalize the DEs completely. In eq. (4.2) we simply rescale s by the internal mass (m 2 > 0), −s/m 2 = u, to deal with a dimensionless variable. If s < 0, u > 0. If s > 0, the Feynman prescription s → s + iε fixes the analytic continuation for u In the case of topologies (c)-(d), the DEs still contain the square roots related to the s-channel threshold at s = 4m 2 . They are rationalized by the usual variable change (see eq. (4.18)), For completeness, we discuss how v varies with s. Holding m 2 > 0, and keeping in mind the Feynman prescription for s > 0, one finds the following cases • For s > 4m 2 , v is on the negative unit interval, and one must replace

Analytic continuation of the master integrals
As discussed in section 4, for all the topologies we start in the patch of the Euclidean region where the alphabet is real and positive (see eqs. (4.8) and (4.23)), and we solve the DEs there. As far as the variables z,z are concerned, the conditions of positivity of the alphabet are the same for all our topologies, i.e. we start from region II (see eq. (5.8)). Regarding the condition on the variables associated to s, i.e. u and v (respectively for topologies (a)-(b) and topologies (c)-(d)), we require It is clear from eq. (5.8) that, if these conditions are satisfied, one does not have access even to the full Euclidean region. Results in the remaining patches of the latter, as well as in the physical regions, are obtained by analytic continuation using the prescriptions described in sections 5.1 and 5.1.
In the present work we performed the analytic continuation numerically, i.e. we assigned to u, v, z,z the vanishing imaginary parts discussed above choosing sufficiently small numerical values. For convenience, we summarize the analytic continuation prescription for the physically interesting cases.
• X 0 → W W : In this region a particle of mass s > 0 decays in two (possibly off-shell) particles with invariant masses p 2 1 > 0 and p 2 Regarding z,z, this corresponds to region II (see eq. (5.8)), therefore no analytic continuation is needed. Furthermore, for topologies (a)-(b) one must replace irrespectively of the value of s, according to eq. (5.18). Instead, for topologies (c) according to eq. (5.21).
• W → W X 0 : This is again a 1 → 2 process involving timelike particles, the only difference being that now The former case corresponds to region IV, the latter to region III (see eq. (5.8)). Therefore, in addition to the analytic continuation in u, v discussed already for the X 0 -decay, one must further use the replacement (5.17) corresponding to region VI, so that z,z inherit the analytic continuation prescription eq. (5.17) from s → s + iε Concerning u, v, the discussion is the same as for X 0 -decay.
corresponding to region VII and VIII respectively, so that z,z inherit from p 2 i → p 2 i +iε the prescription (5. 16) z → z + iε ,z →z − iε .
Since s < 0, no continuation is due on u, v.

Conclusions
In this paper we have computed the two-loop master integrals required for the leading QCD corrections to the interaction vertex between a massive neutral boson X 0 , such as H, Z or γ * , and pair of W bosons, mediated by a SU (2) L quark doublet composed of one massive and one massless flavor. We considered external legs with arbitrary invariant masses. The master integrals were computed by means of the differential equation method.
After identifying a set of master integrals obeying a system of equations which depends linearly on the space-time dimension d, we used the Magnus exponential in order to to find a novel set of integrals that, around d = 4 dimensions, obey a canonical system of differential equations. The canonical master integrals were finally obtained as a Taylor series in = (4 − d)/2, up to order four, with coefficients written as combination of Goncharov polylogarithms, respectively up to weight four.
In the context of the Standard Model, our results are relevant for the mixed EW-QCD corrections to the Higgs decay to a W pair, as well as the production channels obtained by crossing, and to the triple gauge boson vertices ZW W and γ * W W .

A Properties of Chen's iterated integrals
In this appendix we recall the main properties of Chen's iterated integrals [34]. We closely follow the notation of [22]. Chen's iterated integrals are defined by where γ is a piecewise-smooth path connecting x 0 to x, • Invariance under path reparametrization. The integral C [γ] i k ,...,i 1 does not depend on the way one chooses to parametrize the path γ.
• Reverse path formula. If the path γ −1 is the path γ traversed in the opposite direction, then • Recursive structure. From (A.1) and (A.3) it follows that the line integral of one d log is defined, as usual, by and only depends on the endpoints x 0 , x.
It is convenient to introduce the path integral "up to some point along γ": given some path γ and a parameter s ∈ [0, 1], one can define the 1-parameter family of paths which differs from eq. (3.12) by the fact that the outer integration (i.e. the one in dt k ) is performed over [0, s] instead of [0, 1]. Having introduced γ s , we can rewrite (3.12) in a recursive manner: In addition, eq. (A.7) can be used in order to derive the identity • Shuffle algebra. Chen's iterated integrals fulfill shuffle algebra relations: if m = (m M , . . . , m 1 ) and n = (n N , . . . , n 1 ), with M and N natural numbers, one has where the sum runs over all the permutations σ that preserve the relative order of m and n.
• Path composition formula. If α, β : [0, 1] → M are two paths such that α(0) = x 0 , α(1) = β(0), and β(1) = x, then the composed path γ = αβ is obtained by first traversing α and then β. One can prove that the integral over such a composed path satisfies • Integration-by-parts formula. The computation of eq. (A.1) requires, in principle, the evaluation of k nested integrals. Nevertheless, we observe that the innermost integration is always reduced to (A.5), so that one has k − 1 actual integrations to perform. For instance, at weight k = 2, we have and one is left with a single integral to be evaluated, either analytically or numerically.
Moreover, one can show that the integration involving the outermost weight g k can be performed by parts, returning The combined use of eqs. (A.12) and (A.13) allows, for instance, a remarkable simplification in the numerical evaluation of weight k ≥ 3 iterated integrals, since the analytic calculation of the inner-and outermost integrals leaves only k−2 integrations to be performed via numerical methods.
• Conversion to GPLs formula. If all letters appearing in a Chen's iterated integral are rational functions with algebraic roots, then the iterated integral can be converted in terms of GPLs.
Suppose we connect the endpoints x 0 = (a 1 , a 2 , a 3 ) and (A.14) the conversion of the iterated integral to GPLs can be achieved by factorizing all letters such that the dependence on the varied variable x i , becomes linear. Any weight-k integral can then be transformed by where w i are weights which may depend on the constant variables along the path γ i .

B Canonical master integrals for
In this appendix we give the explicit expression of the kinematic coefficients appearing in the definition of canonical master integrals defined by eq. (4.5) and eq. (4.21).

B.1 Topologies (a)-(b)
The coefficients of the canonical MIs listed in eq. (4.5) are given by where the Källén function λ is defined in eq. (4.6).

B.2 Topologies (c)-(d)
The coefficients of the canonical MIs listed in eq. (4.21) are given by where ρ is defined after eq. (4.21) and the Källén function λ is given in eq. (4.6).

C d log-forms
In this appendix we collect the coefficient matrices of the d log-forms for the master integrals (4.5) and (4.21).