NNLO QCD corrections to the $q^2$ spectrum of inclusive semileptonic $B$-meson decays

We calculate the next-to-next-to-leading order QCD corrections to the leptonic invariant mass ($q^2$) spectrum of semileptonic $b \to c$ inclusive decays, taking into account the mass of the charm quark and the charged lepton in the final state. We obtain analytic results in terms of generalized polylogarithms and present numerical studies of the $\mathcal{O}(\alpha_s^2)$ corrections to the $q^2$ spectrum of $b \to c \ell \bar \nu_\ell$ decays, for $\ell =e, \mu$ and $\tau$, in the kinetic scheme. Our computation can be used to incorporate the recent measurements of $q^2$ moments by Belle and Belle II into global fits of inclusive semileptonic $B$-decays.


Introduction
The Heavy Quark Expansion (HQE) has become a pillar in the theoretical description of inclusive decays of heavy hadrons, allowing the derivation of precise predictions with reliable estimates of the uncertainties.One of its main applications is the study of the inclusive decay B → X c ℓν ℓ .Thanks to its relatively large rate and clean experimental signature, studies of inclusive semileptonic decays have lead to precise determinations of the magnitude of the Cabibbo-Kobayashi-Maskawa matrix element V cb .
The HQE allows to describe both the total decay rate and various kinematic distributions as a double series expansions in the strong coupling constant α s and Λ QCD /m b .Various measurements of moments of the charged lepton energy and the hadronic invariant mass in B → X c ℓν ℓ decays have been performed by BABAR [1,2], BELLE [3,4], CLEO [5], CDF [6] and DELPHI [7].The comparison of experimental measurements with the predictions calculated within the HQE has lead to determinations of |V cb | with a 1.1% accuracy [8][9][10] and the so called "HQE parameters", the non-perturbative matrix elements, such as µ 2 π , µ 2 G , ρ 3 D and ρ LS .
The HQE parameters are important theoretical inputs not only for the extraction of |V cb |, but also for the extraction of |V ub | from B → X u ℓν ℓ decays since the moments of the shape functions can be related to the HQE parameters extracted in b → c decays.Moreover predictions for other kinds of processes require precise knowledge of the HQE parameters, as for instance the B-meson lifetimes [11,12] and the rare decay B → X d,s ℓ l [13,14].
An alternative method for the determination of |V cb | has been proposed in Ref. [15] and is based on the measurement of the leptonic invariant mass (q 2 ) spectrum and the branching ratio as a function of a lower cut on q 2 .These observables are invariant under reparametrization, a symmetry within the HQE reflecting Lorentz invariance of the underlying QCD.Reparametrization invariant (RPI) quantities depend on a reduced set of HQE parameters [15,16].The smaller set of parameters necessary in a global fit of these observables (eight instead of 13 up to 1/m 4 b ) enabled to extract |V cb | and the HQE parameters up to 1/m 4  b in a completely data-driven way [17], based on recent measurements of q 2 moments by Belle [18] and Belle-II [19].The new measurements of the q 2 moments have also been included in a global fit of |V cb | [10], together with the lepton energy and M X moments, finding that the new data are compatible with the other measurements and slightly decreasing the uncertainty on the HQE parameters and on |V cb |.Recently, the RPI operator basis has been extended up to order 1/m 5 b [20].Given the precision achieved by experimental measurements, which show a percent or even sub-percent relative accuracy for certain observables, a good control of perturbative and non-perturbative effects in the HQE is mandatory, also in light of the rather large value of the strong coupling constant α s (m b ) ≃ 0.22.At the partonic level, next-to-leading order (NLO) corrections are available from Ref. [21][22][23].The O(α s ) triple differential distribution up to the power-suppressed terms µ 2 π and µ 2 G were presented in Ref. [24][25][26], while NLO corrections proportional to ρ D are available only for the total rate and the q 2 spectrum [27].
The O(α 2 s ) corrections to the free quark decay b → X c ℓν ℓ are also required for a consistent theoretical description and to match the experimental accuracy.The next-to-next-to-leading order (NNLO) corrections to the hadronic invariant mass and charged-lepton energy moments have been calculated in [28,29].No result is available for the q 2 spectrum for arbitrary values of q 2 , only for q 2 = 0 [30], q 2 = (m b − m c ) 2 [31,32] and q 2 = m 2 c [33].While the calculations in these three special points allowed the authors of Ref. [33] to estimate the non-BLM corrections at O(α 2 s ) with a relative 30% uncertainty, their result is unsuited to calculate higher moments of the q 2 spectrum with sufficient precision. 1Analytic expressions up to O(α 3 s ) for the q 2 moments without threshold selection cuts have been presented in Ref. [34], while Ref. [10] presented an evaluation of the α 2 s β 0 corrections, utilizing the BLM correction to the triple differential rate from Ref. [23].
The goal of this paper is to present the complete NNLO QCD corrections to the q 2 spectrum of b → X c ℓν ℓ .At variance with the numerical approach used in [28,29], based on sector decomposition, recent developments in analytic approaches to multi-loop computations allow us to calculate the differential rate w.r.t.q 2 in an analytic form and write it in terms of generalized polylogarithms (GPLs) [35,36].Our results can be used to calculate the NNLO corrections to the q 2 moments with arbitrary cuts on q 2 .The inclusion of the results presented in this paper into global fits will allow to better assess the theoretical uncertainty in the prediction.
The paper is organized as follows.In Section 2 we present the details of the calculation, in particular we discuss how we obtain analytic results for the three-loop master integrals at NNLO.Section 3 presents our numerical results for the differential rate and the moments in the on-shell scheme and the kinetic scheme.We will discuss also the impact on the decay B → X c τ ντ .We conclude in Section 4.

Details of the calculation
Let us now discuss the details of the calculation.We consider the semileptonic decay of a bottom quark mediated by the weak interaction where X c generically denotes a state containing a charm quark, plus additional gluons and/or quarks.The mass of the charged lepton ℓ is denoted by m ℓ while the neutrino is considered massless.The masses of the bottom and charm quark are m b and m c , respectively, and we introduce their ratio ρ = m c /m b .In the following we study the spectrum of the leptonic invariant mass q 2 = p 2 L with p L = p ℓ + p ν .We begin by writing the differential rate w.r.t.q 2 as where q2 = q 2 /m 2 b and F i stands for the differential decay rate at leading, next-to-leading and next-to-next-to-leading order, respectively.The functions F 0 and F 1 are known since a long time [21].For the power corrections up to 1/m 3 b the NLO corrections have been presented in Refs.[27,37,38].The quark masses m b and m c are renormalized in the onshell scheme.The strong coupling constant α s = α (5) s (µ s ) is renormalized in the MS scheme with five active flavours with µ s being the renormalization scale.
In order to calculate the NNLO QCD corrections to the q 2 spectrum, we follow the method of Refs.[27,37].The idea is to consider the differential rate in the presence of a constraint on q 2 .The phase-space decomposition suitable for b → X c ℓν ℓ is carried out by assuming a sequence of two two-body decays.First, the bottom quark decays into an off-shell W -boson and the charm quark, then the virtual W -boson decays into the lepton and the neutrino: where the integration is performed in a d-dimensional space, with d = 4 − 2ϵ.W µν and L µν are the hadronic and leptonic tensors.The element of n-body phase space is given by The constraint δ(p 2 L − q 2 ) enforces the dilepton system to have an invariant mass equal to q 2 .Since the hadronic tensor depends only on p L and p b , we can integrate the leptonic tensor with respect to the phase-space of charged lepton and neutrino [38]: Higher order terms in ϵ are not necessary since the leptonic tensor does not enter into the renormalization, i.e.L µν is always contracted with the renormalized hadronic tensor.After inserting Eq. ( 5) into the differential rate formula, we represent the δ function as the imaginary part of a propagator: In other words, we treat δ(p 2 L − q 2 ) as an on-shell condition for a "fake" particle.Next we apply the reverse unitarity method [39] and map the calculation of the various interference terms integrated over the final state phase-space into the evaluation of "cuts" of forward b → b scattering amplitudes.In this way, the differential rate can be obtained from the the imaginary part of a b → b two-point amplitude, where the two leptons with constrained invariant mass are replaced by a fake (colorless) spin-1 particle with mass equal to q 2 (see Fig. 1).Note that with this method, we need to consider loop diagrams with one loop ) is mapped into the cuts of forward b → b scattering amplitude where the two leptons are replaced by a fake spin-1 massive particle.less compared to the original diagrams, i.e. to calculate the q 2 spectrum at LO, NLO and NNLO we have to consider one-, two-and three-loop diagrams.One the other hand, the Feynman integrals now depend on two dimensionless variables: ρ and q2 .
Let us now discuss the technical details of the calculation.We generate with qgraf [40] one, two and three loop diagrams (as the ones on the r.h.s. of Fig. 1) and use Tapir [41] to map each diagram to a predefined integral family.We use the program exp [42] to rewrite the output to FORM [43] notation.In this way we express the three-loop b → b amplitude as linear combinations of scalar Feynman integrals with nine indices, where eight correspond to the exponents of propagators and the remaining one to the exponent of an irreducible numerator.In total we have 21 integral families at three loops.
Before performing the IBP reduction, we find with a basis of master integrals such that the denominators in the reduction tables completely factorize into polynomials depending either on ρ and q2 or d (see Ref. [44,45]).To construct this basis, we first reduce a set of seed integrals up to two dots and one scalar product for every integral family individually with the help of Kira [46,47] and Fermat [48].As initial basis we simply take the default master integrals.These reduction tables then serve as input to search for a good basis for every family with the help ImproveMasters.mdeveloped in Ref. [44].
The IBP reduction of the integrals appearing in the amplitude is then performed with Kira.First, we reduce the integrals for every family individually to the good basis of this family.Then we employ symmetries between the families to reduce the number of master integrals.Afterwards we identify the master integrals which have an imaginary part while setting to zero those which are purely real, e.g.tadpole integrals.We find one, six and 98 master integrals at one, two and three loops.
We solve the master integrals in an analytic way by using the method of differential equations [49,50].We establish a set of differential equations by differentiating the 98 master integrals with respect to ρ and q2 and reducing the resulting integrals again to master integrals with Kira.We obtain a system of the form where J = (J 1 , . . ., J 98 ) is the set of master integrals, ϵ is the dimensional regularization parameter while M ρ and M q 2 are 98 × 98 matrices, rational in ρ, q2 and ϵ.For the decay rate, we need only the imaginary parts of the master integrals, i.e. the sum of all possible cuts.Before discussing their solution, it is useful to divide them into three classes according to their cuts: (i) master integrals with cuts only through one charm quark propagator (and the fake particle with mass q 2 ).A sample integral is shown in Fig. 2a.
(ii) integrals where the cuts go through three charm propagators (see Fig. 2b).
(iii) integrals with both kind of cuts, i.e. one and three charm cuts (see Fig. 2c).
We would like to bring the system of differential equations in canonical form (or ϵ-form) [51] and express the master integrals in terms of GPLs.Class II contains integrals beyond GPLs such as the sunrise diagram in Fig. 2b with two unequal masses.However, as observed already for the analytic calculation of the total rate at NNLO presented in Ref. [52], the contribution given by cuts through three charm quarks can be neglected for realistic values of the ratio m c /m b .It corresponds to the decay channel b → cccℓν ℓ , which is rare and present only if 0 < ρ < 1/3.For m c /m b ≃ 0.25, the branching ratio is very small, O(10 −7 ), because of the phase-space suppression and totally negligible compared to the current experimental accuracy.
Our strategy therefore is to solve the differential equations by considering only the cuts with one charm quark and neglecting the cuts with three charm quarks.To this end, we can effectively remove the masters in class II, while for class III we pick up only the one-charm cut contributions in the boundary conditions.
After such simplifications, the number of master integrals reduces to 87 and we proceed to find a basis transformation T(ρ, q2 , ϵ) such that the masters in the new basis J = T J satisfy a set of differential equations in canonical form.We change variables from ρ and q2 to To find a rational transformation T, we use Libra [53], which implements the Lee's algorithm [54], in combination with Fermatica [55] to speed up the matrix transformations via the interface to the CAS program Fermat.
As a first step, we find suitable transformations acting on the diagonal blocks that put them into ϵ-form.At NLO, we observe that it is possible to bring the system of differential equations in canonical form when written in terms of the variables x ± , however at NNLO this is not possible anymore.The eigenvalues of the residue matrices of some blocks are of the form aϵ ± 1/2, with a an integer number.Balanced transformations can only raise or lower the eigenvalues by an integer, so in order to bring such blocks into ϵ-form, we need to apply an additional variable change.We find that all such eigenvalues with half-integers can be removed by switching from x + and x − to the variables u and v defined by with 0 < v ≤ u ≤ 1.However, we find it convenient not to perform the variable transformation globally because this would increase the degree of the poles in the residue matrices, making the reduction to ϵ-form computationally more challenging.Instead we take advantage of the Notations mechanism implemented in Libra.The idea is to work with a system still expressed in terms of x + and x − , but with the introduction of the notation In each subblock, all variables x + , x − and u appear.However the dependence of the latter on x + and x − is always taken into account when computing derivatives and it is automatically simplified to first-order polynomial in u.After bringing all diagonal blocks to ϵ-form, Libra automatically reduces the off-diagonal blocks to Fuchsian form and finds a suitable transformation independent ofn x + and x − to factorize out ϵ.In the end, we bring the system of differential equation in canonical form: where A i are 87 × 87 matrices with rational numbers and α i are the letters Since all letters are linear in v, we write the solution of the differential equation in terms of GPLs in v and letters that might depend on u, as well as Harmonic polylogarithms in u.
The boundary conditions to the differential equations are obtained using the auxiliary mass flow method [56,57] as implemented in AMFlow [58].We compute all 87 master integrals in four different kinematic points: These points are sufficient to fix all boundary constants and at least have two additional kinematic points to check the resulting integrals. 2The master integrals are computed with sufficient numerical precision in order to obtain the boundary constants of the differential equations written in terms of transcendental numbers using the PSLQ algorithm [59].

Results
After the evaluation of the master integrals (one-charm cuts), we insert them into the amplitude and perform the wave function and mass renormalization in the on-shell scheme [60][61][62][63][64], while we use MS for the strong coupling constant.
Our main results are the analytic expressions for the functions F 0 , F 1 and F 2 in the differential decay rate in Eq. ( 2).They are written in terms of GPLs depending on u and v, as defined in Eq. ( 9), which can be evaluated numerically to high accuracy, e.g. with GiNaC [65] and PolyLogTools [66].The explicit expressions for F 0 , F 1 and F 2 are given as ancillary files [67].
Moments of the q 2 spectrum are defined by where ).The moments can be expressed as series expansions in α s : From the expressions for F i , we calculate the coefficients in the perturbative expansion via one-dimensional numerical integrations of the functions F i : where which reduces to u min = ρ for q2 cut = 0 .We also define the normalized q 2 moments as and the centralized moments as Moreover, with qn we denote qn = q n /m 2n b .

On-shell scheme
We now present our results for the centralized moments in the on-shell scheme.We set the quark masses to m OS b = 4.6 GeV and m OS c = 1.15 GeV and numerically evaluate the coefficients in the perturbative expansions for Q n , with n = 1, . . . ., 4. Afterwards, we reexpand the ratios in Eqs. ( 19) and (20) in α s up to second order.Our results for the moments without cuts (q 2 cut = 0 GeV which are in good agreement with the NNLO results presented in Ref. [34].For a cut of q 2 cut = 3 GeV

Kinetic scheme
We discuss now the impact of higher-order QCD corrections to q 2 moments once a shortdistance mass scheme is adopted for the quark masses.We concentrate on the kinetic scheme employed in the global fits of Refs.[8-10, 17, 68].In this scheme the on-shell mass of the bottom quark is replaced by the kinetic mass [64,[69][70][71] via the relation while the charm quark mass is converted to the MS scheme.At the same time, we include the contribution from power corrections to the moments up to 1/m 3 b (the relevant expressions can be retrieved from [15,27]).In the kinetic scheme, we redefine the HQE parameters µ 2 π and ρ 3 D in the following way: Note that µ 2 π drops out for centralized q 2 moments, leaving a dependence only on ρ 3 D .The perturbative version of µ 2 π and ρ 3 D up to O(α 3 s ) can be found in the Appendix of Ref. [64].The Wilsonian cutoff µ plays the role of scale separation between the short-and longdistance regimes in QCD.
In order to present our benchmark predictions of the q 2 moments for validation and comparison, we report the series expansion for centralized moments with q 2 cut = 0 GeV 2 and q 2 cut = 4 GeV 2 .We adopt the HQE parameter definitions employed in Refs.[8,9,72] (the so-called perp basis).We use scheme (A) as defined in Ref. [34]: in a first step the expressions for centralized moments are obtained in the on-shell scheme where we retain terms up to O(α 2 s ) at partonic level (1/m 0 b ) while we discard higher QCD corrections for the subleading power corrections.Afterwards we apply the transition to the kinetic scheme.
We set the renormalization scale of the strong coupling constant µ s = m kin b and use α In the following we present the results for the various contribution at leading order in the 1/m b expansion for two different values of q 2 cut .We do not report the contribution from power suppressed terms, however the terms originating from [ρ 2 D (µ)] pert are included.Our results for the moments for q 2 cut = 0 GeV We performed a comparison also with the corrections of order α 2 s β 0 (the so-called BLM corrections [73]) recently presented in Table 1 of Ref. [10] and found good agreement.From the knowledge of the complete NNLO corrections, we observe that in the kinetic scheme the non-BLM contribution to the moments at O(α 2 s ) has in general the opposite sign of the BLM contribution, and is of comparable size.We conclude that the BLM approximation tends to overestimate the NNLO corrections, especially in case one uses m c (3 GeV) as reference mass for the charm quark.
Let us now discuss the size of the NNLO corrections and the impact on the global fits for V cb .In Figures 3 and 4 we show our results for the first four centralized moments as a function of q 2 cut .The predictions are compared with the Belle and Belle II measurements [18,19].At variance with the numerical values given in Eqs. 26 and 27, in the plots we adopt the RPI basis from Refs.[15,16] and the values from the fit in Ref. [17] for the HQE parameters, m b and m c .
The green curves correspond to the LO prediction with power corrections up to 1/m 3 b .The blue curves include QCD NLO correction up to 1/m 3 b , where we also include the α s corrections to µ 2 G and ρ 3 D calculated in Ref. [27].The red curves, compared to the blue ones, additionally include the NNLO corrections at leading order in 1/m b calculated in this article, which are denoted by NNLO' in the plots.The error bands are obtained by varying the renormalization scale in the range m kin b /2 < µ s < 2m kin b and choosing µ s = m kin b as reference scale for the central value.We do not show the parametric uncertainty stemming from the HQE parameters.The lower panel in each plot shows the ratio between the prediction at NNLO and NLO.
In Figure 3, we show the moments obtained with charm mass at a scale of 2 GeV, which is the default choice in the fit in Ref. [17].We observe that the NNLO corrections shift the prediction for q 1 and q 2 by a few percent in the low q 2 cut range.For the third and fourth moment the impact is larger and close to a 10-15% effects.The relative contribution at   [15,16] and values from the fit in Ref. [17].Measurements from Belle [18] and Belle II [19].higher values of q 2 cut becomes larger since the LO central value tends to vanish close to the end point.Note that the use of m c (2 GeV) leads to accidentally small corrections at O(α s ) for all the moments.For q 1 and q 2 one can observe the overlap between the blue and green bands, while the red lines are much more separated.Consequently, scale-variation alone does not provide a reliable uncertainty estimate for the NLO prediction.In fact in this approximation, the scale uncertainty comes only from the variation of α s .Since α s is multiplied by a small number in case one uses m c (2 GeV), a rather small uncertainty is obtained.Improving the prediction from the NLO to the NNLO, we observe that the O(α 2 s ) coefficient is not suppressed anymore and therefore NNLO error bands become larger than the NLO ones.
A better behaviour of the perturbative series is observed utilizing instead m c (3 GeV) =  b [15,16] and values from the fit in Ref. [17].Measurements from Belle [18] and Belle II [19].0.993 GeV.The results are presented in Figure 4. Notice that now the O(α 2 s ) corrections are smaller than the O(α s ) ones, indicating a much better behaved expansion.For the first and second moment the scale uncertainty is reduced from NLO to NNLO and the error bands overlap.Even though the error bands overlap for the third and fourth moment, we still observe a larger uncertainty at NNLO.This can be explained by the fact that the third and fourth moments, especially at high values of q 2 cut , receive sizable contributions form the power suppressed terms, which include perturbative corrections only up to NLO.The uncertainty is not reduced in this case because cancellation of the µ s dependence at partonic level is spoiled by the lower accuracy in the perturbative expansion of the 1/m 2 b and 1/m 3 b terms.Another notable effect observed in Fig. 3 is that after inclusion of the NNLO corrections the curves move to values higher than the experimental data points.This does not indicate a tension between data and theory.In fact, since we use the HQE parameters from a fit accurate only up to NLO at partonic level [17], we naturally expect the NLO curves to show better agreement with the data.Notice that the blue curves include also the NLO corrections at 1/m 2 b and 1/m 3 b , which was not the case in Ref. [17].The major effect in global fits after including the NNLO corrections would be a change of the favoured values of the HQE parameters in order to shift downwards the red curves and accommodate the predictions with the data.In particular, since ρ D has a major impact in the q 2 moments and enters with negative coefficients, a fit with NNLO corrections would prefer higher values for ρ D compared to Ref. [17].
In addition to the analysis of the moments in the RPI basis, we analyzed also the prediction with the fit setup from Ref. [10] and using the perp basis for defining the HQE parameters.We reached similar conclusions for what concerns the use of a charm mass at a scale of 2 GeV or 3 GeV: a charm mass at 3 GeV yields a better behaviour of the pertubative series while a 2 GeV charm mass underestimates uncertainties at NLO.For completeness, we report in Fig. 5 our results for a charm mass m c (2 GeV) = 1.092GeV, the default scheme in Ref. [10].We also mention that since the α 2 s β 0 corrections overestimates the α 2 s corrections, the NNLO predictions shown by the red curves in Fig. 5 lie below the experimental data.A curve showing the q 2 moment prediction with only the α 2 s β 0 corrections would appear above the red curves, more in agreement with data.We conclude that the inclusion of the complete NNLO corrections in the fit of [10] would bring the red curves upwards, towards the experimental data, preferring a lower value for ρ D (also in the perp basis the coefficient of ρ D is negative).

Decay into a massive tau lepton
While the O(α 2 s ) contributions to the q 2 spectrum are most relevant for the decay into light leptons, our expressions also apply to inclusive b → X c τ ντ decays.The first measurement of the ratio was recently performed by the Belle II experiment [74].The current level of experimental precision is severely limited by systematic uncertainties related to the modelling of B → Xτ /ℓν τ /ℓ decays.However, recent progress in the description of B meson decays into excited charm meson states [75] will allow to address this issue in a data-driven way in the future, making the inclusion of O(α 2 s ) corrections relevant.In the on-shell scheme, our results for the integrated b → X c τ ντ decay rate without cut on q 2 agrees with Ref. [29].In the kinetic scheme, with m kin b (1 GeV) = 4.526 GeV, m c (3 GeV) = 0.993 GeV and at leading order in 1/m b , we obtain  The perturbative convergence is similar to the case of the q 2 moments.This prediction can be improved by incorporating 1/m b suppressed terms at LO and NLO [38].
Our results allow to obtain predictions for R(X c ) with a lower cut on q 2 .For q 2 cut = 6 GeV 2 , we obtain R(X c ) q 2 >6 GeV 2 = 0.350 1 − 0.782 The higher order corrections clearly become more relevant if a cut is introduced.Furthermore, the ratio increases with increasing q 2 cut , as terms proportional to m 2 τ /q 2 and phase-space effects become less relevant.Consequently, it could be advantageous to perform measurements of R(X) using a lower cut on q 2 in the future, as it enriches the fraction of B → Xτ ντ decays, rejects backgrounds and cuts away the regions where the B → Xℓν modelling is most problematic.In addition a lower cut on q 2 allows for the improved inclusion of momentum requirements on signal leptons due to detector thresholds and the reduction of uncertainties associated to the modelling of final state radiation [76].Most of the analysis strategy of the Belle measurement of the q 2 moments [18] could thus carry over to a future measurement of R(X), with the exception of a cut on the difference of the missing energy and the missing momentum in a given event. 3

Conclusions
In this article we presented the complete NNLO QCD corrections to the q 2 spectrum of inclusive semileptonic B decays.The differential rate with respect to the leptonic invariant mass q 2 is obtained by calculating the imaginary part of the b → b 2-point function in the presence of a constraint on q 2 , which can be implemented in a convenient way by replacing the charged lepton-neutrino loop with a fake particle with mass q 2 .After reduction to master integrals, we leverage the method of differential equations to calculate the decay rate analytically.To this end, we restricted the calculation to the cuts through only one charm line which allowed us to bring the system of differential equations in canonical form and express the master integrals in terms of GPLs.
In light of the recent measurements of Belle and Belle II of the q 2 moments, we studied the impact of the NNLO corrections on the moments as a function of q 2 cut .We observe that the O(α 2 s ) corrections are sizable especially for the default choice of the charm mass m c (2 GeV) in the global fits of [10,17].The relative ratio to the NLO prediction reaches the 1-5% level, depending on the cut on q 2 , while for the higher moments it is larger and of about 10-20%.For a charm mass at the scale of m c (3 GeV) we observe a better behaviour of the perturbative series, with the expected reduction of theoretical uncertainties when including the O(α 2 s ) corrections.We applied our results also to the decay into a tau lepton, and proposed a measurement of R(X c ) using a lower cut on q 2 to enrich the fraction of B → X c τ ντ events.We provide our results for the q 2 spectrum in electronic form as ancillary files.They can be employed to incorporate the NNLO corrections into global fits of inclusive semileptonic B-decays, in particular to take advantage of the recent measurements of q 2 moments by Belle and Belle II.

δ(p 2 L − q 2 Figure 1 :
Figure 1: Phase space intergration in the presence of the constraint on the leptoinc invariant mass δ(p 2 L − q 2) is mapped into the cuts of forward b → b scattering amplitude where the two leptons are replaced by a fake spin-1 massive particle.

Figure 2 :
Figure 2: Sample of master integrals.Black, red and green lines denoted massive propagators with mass equal to m b , m c and q 2 .Dashed lines are massless propagators.The master integrals can have cuts through only one charm propagator (a), only three charm propagators (b) and both kind of cuts (c).
kin b ) as expansion parameter, i.e. we decouple the bottom quark from the running of α s , and we reexpand the leading 1/m b term in α

Figure 3 :
Figure 3: The first four q 2 moments of B → X c ℓν ℓ as a function of the lower cut q 2 cut .The heavy quark masses are m kin b (1GeV) = 4.562 GeV and m c (2 GeV) = 1.094GeV.For the HQE parameter we adopt the RPI basis up to 1/m 3 b[15,16] and values from the fit in Ref.[17].Measurements from Belle[18] and Belle II[19].

Figure 4 :
Figure 4: The first four q 2 moments of B → X c ℓν ℓ as a function of the lower cut q 2 cut .The heavy quark masses are m kin b (1 GeV) = 4.562 GeV and m c (3 GeV) = 0.993 GeV.For the HQE parameter we adopt the RPI basis up to 1/m 3 b[15,16] and values from the fit in Ref.[17].Measurements from Belle[18] and Belle II[19].

Figure 5 :
Figure 5: The first four q 2 moments of B → X c ℓν ℓ as a function of the lower cut q 2 cut .The heavy quark masses are m kin b (1 GeV) = 4.573 GeV and m c (2 GeV) = 1.092GeV.For the HQE parameter we adopt the perp basis up to 1/m 3b[72] and values from the fit in Ref.[10].The central values are obtained for a renormalization scale µ s = m kin b /2.Measurements from Belle[18] and Belle II[19].