Five-particle contributions to the inclusive rare $\bar B \to X_{s(d)} \, \ell^+\ell^-$ decays

We calculate tree-level contributions to the inclusive rare $\bar B \to X_{s(d)} \, \ell^+\ell^-$ decays. At the partonic level they stem from the five-particle process $b \to s(d) \, q \bar q \, \ell^+\ell^-$, with $q \in \{u,d,s\}$. While for $b \to d$ transitions such five-body final states contribute at the same order in the Wolfenstein expansion compared to the three-body partonic decay, they are CKM suppressed in $b \to s$ decays. In the perturbative expansion, we include all leading-order contributions, as well as partial next-to-leading order QCD and QED effects. In the case of the differential branching ratio, we present all results completely analytically in terms of polylogarithmic functions of at most weight three. We also consider the differential forward-backward asymmetry, where all except one interference could be obtained analytically. From a phenomenological point of view the newly calculated contributions are at the percent level or below.


Introduction
Despite tremendous effort, the experiments at the LHC have to date not found any evidence for a direct signal of physics beyond the Standard Model (SM). In the absence of direct new-physics signals, indirect searches via low-energy observables become ever more important. In case the latter are suppressed or forbidden at tree level in the SM, they are capable of probing high scales via virtual effects of new degrees of freedom. Depending on the process, the scales probed in low-energy observables might even outrange those accessible via on-shell production.
As a matter of fact, all currently observed tensions between experimental data and SM predictions -the so-called anomalies -are in low-energy observables located in the (quark or lepton) flavour sector of the SM. Among the most prominent ones, there are P 5 inB → K * + − , R D ( * ) , R K ( * ) , and (g − 2) µ [1]. All except the last one stem from exclusive decays of heavy mesons, and experimental measurements rely mostly on data from the B-factories [2][3][4], LHCb [5][6][7][8] , and partially ATLAS [9] and CMS [10].
One useful way to shed light on the nature of the anomalies is a cross-check via the corresponding observables in the inclusive modes. With the first data-taking at Belle II in sight, there is a unique opportunity to get access to these modes with sufficient precision both on the theoretical and experimental side. On the theory side, obtaining precise predictions amounts to including higher-order perturbative corrections on the one hand, but also multi-particle contributions on the other. In the case of rare and radiative flavour-changing neutral current transitions the focus in the past was largely on b → s decays such asB → X s γ andB → X s + − , which have both reached a highly sophisticated level in terms of inclusion of perturbative corrections (for the latest comprehensive analyses, see [11] and [12], respectively). While inB → X s γ four-body contributions which at the partonic level amount to b → s qq γ (with a light quark q ∈ {u, d, s}) have been calculated to leading [13] and next-to-leading order [14,15], the corresponding five-particle b → s qq + − modes toB → X s + − are yet unknown. Besides b → s decays, also b → d transitions will become relevant in the Belle II era. They are interesting on their own grounds, because contrary to b → s, all three sides of the unitarity triangle for b → d transitions are of the same order O(λ 3 ) in the Wolfenstein expansion parameter λ ≈ |V us |. This means in particular that matrix elements of the effective operators P u 1,2 are not CKM suppressed compared to those of their P c 1,2 siblings. On general grounds, one can therefore expect quite sizeable CP-violating effects in b → d modes.
While there are plenty of papers dealing with the inclusiveB → X s + − decay, there is only little dedicated literature on theB → X d + − counterpart. The next-tonext-to-leading order (NNLO) corrections to the matrix elements have been computed in [16,17], and also the latest phenomenological study dates back fifteen years [16]. In view of prospects on the experimental side, also a theory update is timely. We start this endeavour in the present article by computing the five-particle b → d qq + − process with q ∈ {u, d, s} at tree level. We compute these corrections for two distributions that are differential in the dilepton invariant-mass squared, namely the branching ratio and the forward-backward asymmetry. It turns out that the matrix elements of the dimension-six operators are real-valued due to the absence of strong phases, and hence do not affect CP-violating observables. The corresponding contributions toB → X s + − are obtained by a simple replacement of the CKM factors.
On the technical side we have to deal with the integration of the squared matrix elements over the five-particle massless phase space. While several parametrisations exist [18][19][20] which all turn out to be useful for particular types of interferences (which we specify in subsequent sections) the analytic integration over the kinematic invariants turns out to be the main challenge since we aim at performing all integrations in an analytical manner. While many of the previous studies of processes with five final-state particles [21][22][23][24] use numerical methods, our paper serves as a proof-of-principle that five-particle massless phase-space integrations can be carried out (almost) completely analytically in terms of polylogarithmic functions of at most weight three. Another recent work which evaluates the five-particle master integrals with massless particles analytically in dimensional regularisation was presented in [25]. That paper integrates the integral kernels over the entire phase space, whereas in the present paper we stay differential in the dilepton invariant-mass squared. This entails on the one hand that we do not encounter any infrared divergences in the present calculation. On the other hand, we do not have the freedom to exploit the full symmetry of the phase-space measure.
This article is organised as follows: In section two we specify the operator basis and explain the power-counting by means of which we decide which terms to include in the present calculation. Section three contains details on the phase-space integration, and in section four we present our results for the differential branching ratio and forwardbackward asymmetry. We finally conclude in section five, including a numerical estimate of the impact of the five-particle contributions. In two dedicated appendices, we illustrate how to analytically integrate sample kernels over the five-particle phase space, and collect the functions that appear in our analytical results.

Effective theory and power counting
The calculation of transition amplitudes in heavy-quark decays is most conveniently performed in an effective field theory where the top quark, the heavy gauge bosons and the Higgs field have been integrated out [26]. In what follows we present the formulas relevant for b → d transitions. The corresponding formulas for b → s are obtained by obvious adjustments of quark fields and CKM factors.

The effective theory
The Lagrangian of the effective theory reads b P u 1,2 , P 3,...,6 b P u 1,2 , P 3,...,6 Figure 1: Contributions of the operators P u 1,2 , P 3,...,6 to the five-body decay, where the lepton pair can be emitted and absorbed from any of the crosses. Insertion (b) is only present for penguin operators P 3,...,6 with q = d. where and q = u, d, s, c, b and l runs over the three charged lepton flavors. The five-particle contributions toB → X d + − amount at the partonic level to b → d + − qq, where by definition the final state is free of charm quarks, i.e. q ∈ {u, d, s}. At leading order in the perturbative expansion, one has to calculate interferences of the operators P u 1,2 , P 3,...,6 at tree level, which contribute at O(α 2 e ) as illustrated in Fig. 1. Here the lepton pair can be emitted and absorbed in the 16 different ways indicated by the crosses. For the penguin operators P 3,...,6 for q = d in addition also the insertion in Fig. 1b contributes.
From a technical point of view, one has to evaluate the squared matrix elements according to the interferences shown in Fig. 1, and subsequently integrate them over the  Figure 2: Contributions of the P u 1,2 , P 3,...,6 interference with P 9 via (a) emission of the quark pair via a gluon (QCD) (b) emission of the quark pair via a photon (QED). Again the photon or gluon can be emitted and absorbed from the different fermion lines as indicated by the crosses.
five-particle massless phase space. Since we stay differential in the dilepton invariantmass squared, we do not encounter any ultraviolet or infrared divergences and hence can perform the entire calculation in four-dimensional space-time. As a result, the interferences between the operators P u 1,2 , P 3,...,6 can be computed completely analytically. When calculating the matrix elements, the following relations turn out to be useful, They hold in four dimensions and shorten the occurring traces over Dirac structures considerably.

Power counting and higher orders
In the presence of electroweak corrections, the perturbative expansion in inclusiveB → X s(d) + − is consistently done in terms of α s = α s /(4π) and κ = α e /α s , as was pointed out in [27]. Following [27,28], we only keep contributions up to O( α 3 s κ 3 ) to the differential decay rate. As discussed above, the leading five-body contributions toB → X d + − arise from the P u 1,2 , P 3,...,6 operators and are O( α 2 s κ 2 ). In order to specify the set of interferences that we take into account beyond O( α 2 s κ 2 ), we also consider possible suppressions by the Wilson coefficients. The numerical values for the Wilson coefficients can be found in Ref. [27]. Through renormalisation group running, the Wilson coefficients C 1,...,8 (µ b ) at a scale µ b ∼ m b start to receive contributions at O( α 0 s κ 0 ), while C 9 and C 10 are given by [27] at µ b = 5 GeV. The Wilson coefficients C 9,10 will be counted as starting from O( α s κ) due to the numerical smallness of their respective first coefficient.
We are aware of the fact that there are one-loop insertions from P 1,...,6 that start contributing at O( α 3 s κ 2 ), i.e. at the same order as the tree-level interferences involving the P 7,...,10 operators. Including these loops would require performing our calculation in D dimensions. Since the aim of the present paper is to perform a first exploration of five-body contributions, inclusion of these loops is beyond the scope of the current analysis. To set the stage for a future calculation of five-body loop corrections, we find it nevertheless beneficial to take into account all tree-level contributions up to O( α 3 s κ 3 ) and to give a numerical estimate of the impact of the higher-order terms.
As discussed, the P 7,...,10 operators contribute via a QCD or QED emission of the qq-pair and interference with P u 1,2 , P 3,...,6 (see Figs. 2a and 2b). While we take all QCD emissions into account, we treat the QED emissions as follows. In case of the branching ratio, we include the full QED emission of the operators P 7, 8,9 . Here it turns out that the interferences where the qq pair is emitted from the lepton pair is zero for the P 7 and P 9 operators by symmetry of the phase-space measure. However, for the operator P 10 , which involves the axialvector current, these particular QED emissions do contribute and turn out to be infrared divergent. We illustrate them in Fig. 3a. To cancel the divergent part of this contribution would require the inclusion of fewer-particle cuts involving loops such as the two-loop contribution of the three-body decay depicted in Fig. 3b. As mentioned above, we presently only work at tree level and therefore omit the full QED emission of P 10 and leave it for future work. For the forward-backward asymmetry, the situation is reversed and actually the P u 1,2 , P 3,...,6 interferences with P 7,9 via QED emission are infrared divergent. By the same argument, we discard these contributions, but keep the infrared finite P u 1,2 , P 3,...,6 interferences with P 8,10 via QED emission. Finally, for P 7 and P 8 we use the scheme independent effective Wilson coefficients which effectively also includes universal penguin-loop corrections. However, not taken into account are the finite, non-universal parts, induced for instance by the cc pair in the penguin loop, as well as those contributions where both the gluon and the photon are emitted from the penguin loop (see [15] for details). Figure 3: Illustration of the P u 1,2 , P 3,...,6 interference with P 10 where the qq-pair is emitted from the charged leptons. The cancellation of infrared divergences in the five-particle cuts (a) requires fewer-particle cuts involving loops as exemplified in (b).
To summarize, we only consider the tree-level contributions of the P u 1,2 , P 3,...,6 operators and their interference with P 7,..., 10 . Infrared divergent parts that require fewer-particle cuts involving loops to render them finite are left for future study. We emphasize that this set of diagrams represents a gauge invariant subset.

Phase-space integration
Having specified all contributions that we include in the present work, one has to calculate the squared matrix element of the process and subsequently integrate it over the massless five-particle phase space (PS). Due to the fact that we stay differential in the rescaled invariant-mass squared of the final-state lepton pair, the latter invariant acts as a regulator of infrared divergences in all interferences that we include. As a consequence, our calculation, and in particular the PS integration, can be done entirely in D = 4 space-time dimensions. The completely differential five-particle PS measure reads It turns out that our squared matrix element only depends on m b , scalar products (p i · p j ), 1 ≤ i < j ≤ 5, and linear factors µνρσ p µ i p ν j p ρ k p σ l . Moreover, the expression for cos θ that we use to project onto the forward-backward asymmetry, can be written in terms of scalar products as well, see eq. (35). The fact that d 3 p i sums over all threemomentum configurations of particle i ensures that, after PS integration, all terms that involve a single -tensor vanish, which also makes all functions of the final result being real-valued. Hence, from a practical point of view, we use PS parametrisations in terms of rescaled invariant masses s ij = 2p i · p j /m 2 b = (p i + p j ) 2 /m 2 b , which are, by construction, dimensionless.
In our four-dimensional setup, there remain eight independent PS variables: The ten different s ij , which are subject to the two constraints (a) that their sum equals to unity, which is a direct consequence of momentum conservation, and (b) that the Gramdeterminant det [2(p i · p j )], i, j = 1, . . . , 5 vanishes. The latter condition holds since in D = 4 dimensions, four light-like vectors are sufficient to span Minkowski space, and hence the five momenta p µ 1 , . . . , p µ 5 must be linearly dependent. In practice, we therefore have to perform seven integrations, keeping in mind that we stay differential in s 45 . To this end, we adopt two different PS parametrisations, one by Kumar (K) [19], and the other one by Heinrich (H) [20]. We give details on both parametrisations in turn below.

Integration according to Kumar
Following Kumar's parametrisation [19], we have where the Lorentz-invariant Mandelstam-like variables are defined by The integration limits are given by Here is the Källén function and λ(1, t 2 , t 2 ) λ(1, s 3 , s 3 ) , Using eq. (11), we can now express the rescaled invariant masses s ij in terms of the integration variables s 1,2,3 , u 1,2,3 and t 2,3 . For all except one invariant, this is straightforward, Equation (19) incorporates already the constraint that the sum of all ten s ij must equal to unity. The remaining invariant s 34 now gets extracted from the condition that the Gram determinant vanishes. Plugging eq. (19) into det [2(p i · p j )] = 0 leaves us with a quadratic equation for s 34 , whose two solutions read where s r/s 34 refer to the rational and square-root part of the solution, respectively. The correct implementation of this two-fold solution in the PS integration then amounts to We note in passing that our results are at variance with some of the earlier works that dealt with the five-particle PS. The authors of [21] obtained an expression for s 34 which captures only the s r 34 part of the full solution. The reason can be traced back to the method proposed in Appendix D of [19], where the delicate point is that only the integration 1 -rather than p µ 4 itself -is a linear combination of p µ b , (p b −p 123 ) µ and (p b −p 23 ) µ , although the conditions represented by the four delta functions in the integration truly have to be satisfied. Solving the four equations corresponding to the four delta functions, we find indeed two valid solutions for p µ 4 , whose sum is a linear combination of p µ b , (p b − p 123 ) µ and (p b − p 23 ) µ . Therefore, what was obtained in [21] is actually the average of the two solutions for s 34 .

Integration according to Heinrich
The parametrisation of Heinrich [20] assumes a form in which the PS measure completely factorises in the eight independent variables, which are labelled t 2 , . . . , t 4 , t 6 , . . . , t 10 and whose integrations all run from 0 to 1. In order to make the present paper self-contained, let us repeat the essential features from [20] in the following.
The derivation starts from the PS parametrisation in D dimensions, where only the constraint s 12 + . . . + s 45 = 1 holds, but not det [2(p i · p j )] = 0. The relation between the rescaled invariant masses and the nine variables t 2 , . . . , t 10 reads and y ± 5 are the solutions of det [2(p i · p j )] = 0. In the transition to D = 4 dimensions the latter constraint is then implemented by letting t 5 assume the values t 5 = 0 or t 5 = 1 only, which renders s 24 = y ± 5 and the number of independent variables is reduced to eight. The four-dimensional PS integral including a matrix element then reads Staying differential in one of the s ij is then simply taken into account by a corresponding delta-function factor, e.g. δ(t 2 t 4 t 6 t 7 − s 14 ), and hence the number of integrations that actually have to be performed is seven, as in the case of the parametrisation (K).

Comparison and concluding remarks
Each of the two parametrisations has its virtues and its drawbacks. At a first glance, the parametrisation (H) seems superior to that of (K): The integration limits are independent of each other and the symmetry under renaming of all final-state particles is manifest, which leaves us with a freedom to choose a labelling that makes the matrix element look as simple as possible. In the parametrisation (K), we use the freedom to rename {1, 2, 3} (corresponding to d, q,q) and {4, 5} (corresponding to the lepton pair). From a practical point of view, the parametrisation (H) is easier when only massless propagators are present in the squared matrix element, whereas we prefer parametrisation (K) whenever massive propagators occur. In appendix A we illustrate with sample kernels how to perform the analytic integration in the (H) and (K) case, respectively. Besides the analytic integration, we perform numerical checks in all cases. To this end we use a Monte-Carlo integration strategy with the Vegas [29] algorithm from the Cuba library [30]. We run the numerical integration for several values of s 45 , typically with a million integration points. This gives us a relative precision of better than one percent, and in many cases even at the per-mill level. In the region of high s 45 , the performance is a bit worse since the results rapidly approach zero due to the lack of PS volume.

Results
In the following, we present our results for the branching ratio and the forward-backward asymmetry, both as functions ofŝ ≡ q 2 /m 2 b = (p 4 + p 5 ) 2 /m 2 b = s 3 = s 45 . For the former observable, all results are obtained completely analytically in terms of transcendental functions of at most polylogarithmic weight three, while for the latter one function required a fit.

Branching ratio
The differential five-particle decay width can be expressed as Following [27], we normalise the decay width to the branching ratio of theB → X c eν decay, which is further expressed in terms of the perturbative expansion of theB → X u eν decay (including power-corrections) and the ratio [31,32] We use C = 0.574 ± 0.019 [33]. Consequently, the expression of the branching ratio reads where Φ u is defined by [12,27] This normalisation strategy is introduced to reduce the uncertainties arising from m b,pole , CKM matrix elements, and phase-space factors involving m c . We leave the details of the error analysis to a later work in which we perform a complete phenomenological study of the inclusiveB → X d + − decay [34]. The indices i and j in (28) run over all the operators in eq. (1). For the Wilson coefficients C i (µ), we take the values at the low scale µ = µ b = 5 GeV. The analytical forms for C i (µ b ) can be found in [27]. As explained in section 2 we use C eff 7 and C eff 8 given by (7). The CKM ratios are given as −ξ * u , for i = 1, 2, j = 3, ..., 10; −ξ u , for i = 3, ..., 10, j = 1, 2; 1, for i, j = 3, ..., 10, . Finally, the tree-level contributions to the differential branching ratio from different operators are summarized in the 10 × 10 matrix where The explicit forms of all the U ( ) i (ŝ) functions can be found in Appendix B. The color factors are defined by (33) with N c = 3, t f = 1/2, C F = 4/3. As discussed in section 2, the lower right 4 × 4 block, the last row and the last column of F(ŝ) are absent. For convenience, the 10 × 10 matrix F(ŝ) is provided in electronic form in the file FA.txt that is attached to the arXiv submission of this article.

Forward-backward asymmetry
We define the forward-backward asymmetry A FB by Here, z = cos θ, with θ the angle between the + and b-quark three-momenta in the dilepton rest frame. It turns out that z can be expressed in terms of the Lorentz invariant momenta products as with s ib ≡ 2p i · p b /m 2 b and s 123 = (p 1 + p 2 + p 3 ) 2 /m 2 b . As anticipated in section 3 the structure of the squared amplitude in all terms that we keep is such that projecting the double differential decay width onto the forward-backward asymmetry by means of the function 3/2 z is equivalent to the projection using sign(z), which we checked by explicit computation.
Using the same normalisation strategy as for the branching ratio (28), the forwardbackward asymmetry is further expressed as with the tree-level contributions from different operators summarized in with The functions U ( ) i (ŝ) are again relegated to Appendix B, and the quantity A(ŝ) is also provided in electronic form in the file FA.txt.  Table 1: Numerical estimates for the branching ratio and the forward-backward asymmetry of b → d + − qq in the low-q 2 region.

Numerical estimate and conclusion
In this section we give the numerical estimates for the branching ratio and the forwardbackward asymmetry of the five-particle decay process b → d + − qq. For each observable we give the integral over bin 1 (q 2 ∈ [1, 3.5] GeV 2 ), bin 2 (q 2 ∈ [3.5, 6] GeV 2 ), and the entire low-q 2 region (q 2 ∈ [1, 6] GeV 2 ). We do not consider the high-dilepton invariantmass region here since there the five-body contributions are negligible because they are suppressed by high powers of (1 −ŝ) due to the lack of PS volume. We use the input parameters from Table 1 of [12], except that for the CKM matrix elements we take the most recent determination of the parameters [35] λ = 0.2251 +0.007 −0.012 , A = 0.825 ± 0.0003,ρ = 0.160 +0.008 −0.007 ,η = 0.350 ± 0.006, which gives Then, we obtain the numerical results for the branching ratio and the forward-backward asymmetry in the three bins, as shown in Table 1. Only the central values are presented, whereas the study of uncertainties is relegated to a forthcoming work focusing on the phenomenological analysis of the inclusiveB → X d + − decay [34]. Keeping only the P u 1,2 , P 3,...,6 interferences among themselves, which start contributing at order O( α 2 s κ 2 ), we find for the low-q 2 integrated branching ratio B(b → d + − qq) [1,6] = 9.60 × 10 −10 . Actually, if we furthermore drop the penguin operators P 3,...,6 we already obtain 9.40 × 10 −10 in the same q 2 region, which demonstrates that the P u 1,2 -P u 1,2 interferences dominate.
Comparing the five-particle branching ratio and the three-particle branching ratio of B → X d + − that we have also estimated, we find that B(b → d + − qq) [1,6] contributes about 1.4% to B(B → X d + − ) [1,6] and B(b → d + − qq) [1,3.5] contributes about 2.5% to B(B → X d + − ) [1,3.5] . Our results are also valid for b → s + − qq, and can be obtained by changing the CKM matrix elements V ud → V us and V td → V ts in Eqs. (28) and (36). The ratios relevant for the calculation are We straightforwardly obtain the corresponding results for b → s + − qq, as listed in Table 2. Compared to the lastest theory prediction for the branching ratio in the low-q 2  Table 2: Numerical estimates for the branching ratio and the forward-backward asymmetry of b → s + − qq in the low-q 2 region.
region [12], the five-particle contribution is at the level of O(0.01%), which shows the expected CKM suppression.
To conclude, we have presented the first study of tree-level five-particle contributions toB → X s(d) + − . While for the b → s transition, such contributions are numerically at the sub-permille level, they play a more pronounced role in b → d where their contribution is estimated to be at the percent level. In light of the anomalies in exclusive b → s systems and the upcoming data-taking at Belle II, a study of inclusiveB → X s(d) + − decays is both timely and relevant since it provides important complementary constraints compared to the widely studied exclusive decays. The present study paves the road for a phenomenological update of both the branching ratio and the forward-backward asymmetry inB → X d + − [34], whose latest analysis dates back 15 years [16]. From a computational point of view, the integrations of the squared matrix elements over the massless five-particle phase-space in four dimensions, while simultaneously staying differential in one of the kinematic invariants, is very challenging. Our paper therefore serves as a proof-of-concept that such integrations can -with one exception -be carried out completely analytically in terms of polylogarithmic functions of at most weight three.

A Details on the phase-space integration
In this appendix we describe how to integrate two sample kernels analytically over the five-particle PS, using the parametrisation (K) and (H), respectively.

A.1 Integration via parametrisation (K)
As mentioned in the main text, the parametrisation (K) turns out to be quite efficient if massive propagators are present. We therefore illustrate the integration which appears for instance in the O u 2 -O u 2 interference. We substitute the rescaled invariant masses s ij = (p i + p j ) 2 /m 2 b according to section 3.1, and apply in particular eq. (21) for the treatment of s 34 .
First, we perform the t 3 -integration, where it turns out that the dependence of K on t 3 is polynomial. After the substitution the χ-integration runs from 0 to 1 and can be carried out in terms of Γ-functions, yielding and similar if there is a polynomial in t 3 in the numerator.
The obtained integrand has a polynomial dependence on u 3 . Since the limits t ± 2 of the t 2 -integration are independent of u 3 , we can trivially interchange the two integrations, and the u 3 -integration is elementary (i.e. one computes an integral function and plugs in the limits). Surprisingly, after these two integrations, the integrand has simply become − s 2 − s 3 .
Looking at its structure, we see immediately that we can perform the t 2 -integration along the same lines as the t 3 -integration before. However, we experienced also cases where the dependence on t 2 at this stage is more complicated. For instance, a 1/t 2 -dependence can occur, which makes the subsequent substitutions considerably more involved. The next integration to be performed is that over u 2 , which is also elementary, but which now introduces a logarithm. We observe, however, that the dependence of the resulting integrand on the variables s 1 and u 1 is via the combination u 1 + s 1 only! We therefore shift the variable u 1 → u 1 − s 1 and interchange the order of the remaining integrations, which results in the new integration limits s 2 = s 3 . . . 1, u 1 = 2 √ s 2 . . . 1 + s 2 , and s ± 1 = u 1 /2 ± u 2 1 − 4s 2 /2. The integration over s 1 is then trivial, yielding the integrand (s 2 − s 3 ) 4s 2 ln The remaining two integrations over u 1 and s 2 , which we perform in this order, are also elementary if a computer algebra system is used. After simplification, the final result reads where f 1 is defined in eq. (63).

A.2 Integration via parametrisation (H)
The parametrisation (H) is efficient in the case of only massless propagators. We therefore demonstrate how to integrate the sample kernel which also appears in the O u 2 -O u 2 interference. To this end, we first relabel the particle indices in the PS parametrisation eq. (23) (but not in H) according to {1, 2, 3, 4, 5} → {5, 3, 2, 4, 1}. Staying differential in the lepton invariant-mass squared is then implemented via the factor δ(t 2 t 4 t 6 t 7 − s 45 ), and implementing the vanishing of the Gram determinant proceeds via eq. (25).
The first two integrations are over t 8 and t 10 . In both cases, the dependence of the integrand on these variables is of the type t a i (1 − t i ) b , with integer or half-integer a and b. Hence, both integrations can be done in terms of Γ-functions. Afterwards, the integrand is polynomial in t 3 and integration results in The next integration over t 9 is elementary and introduces a logarithm. We observe that in the resulting integrand, the dependence on the variables t 2 and t 4 in the denominator and in the δ-function is via the combination t 2 t 4 only! We therefore substitute t 2 = w/t 4 and integrate over t 4 next, i.e.
The resulting expression is short, 1 8t At this stage, we apply the above trick again twice by first substituting t 7 = x/t 6 (followed by integration over t 6 from x . . . 1) and subsequently setting w = z/x, followed by integrating over x from z . . . 1. Both integrations are elementary with a computer algebra system at hand. What remains is a trivial integration over z due to the factor δ(z − s 45 ). After simplification, the final result reads

B Functions
Here we present the U ( ) functions appearing in (32) and (37). The electric charge factors are Q L = −1, Q d = Q s = −1/3, Q u = +2/3, and the sum over q runs over u, d, s.
The term proportional to Q u in U 8 was obtained from a least-square fit.