Constraints for the semileptonic $B \to D^{(*)}$ form factors from lattice QCD simulations of two-point correlation functions

In this work we present the first non-perturbative determination of the hadronic susceptibilities that constrain the form factors entering the semileptonic $B \to D^{(*)} \ell \nu_\ell $ transitions due to unitarity and analyticity. The susceptibilities are obtained by evaluating moments of suitable two-point correlation functions obtained on the lattice. Making use of the gauge ensembles produced by the Extended Twisted Mass Collaboration with $N_f = 2+1+1$ dynamical quarks at three values of the lattice spacing ($a \simeq 0.062, 0.082, 0.089$ fm) and with pion masses in the range $\simeq 210 - 450$ MeV, we evaluate the longitudinal and transverse susceptibilities of the vector and axial-vector polarization functions at the physical pion point and in the continuum and infinite volume limits. The ETMC ratio method is adopted to reach the physical $b$-quark mass $m_b^{phys}$. At zero momentum transfer for the $b \to c$ transition we get $\chi_{0^+}(m_b^{phys}) = 7.58\,(59) \cdot 10^{-3}$, $\chi_{1^-}(m_b^{phys}) = 6.72\,(41) \cdot 10^{-4}$ GeV$^{-2}$, $\chi_{0^-}(m_b^{phys}) = 2.58\,(17) \cdot 10^{-2}$ and $\chi_{1^+}(m_b^{phys}) = 4.69\,(30) \cdot 10^{-4}$ GeV$^{-2}$ for the scalar, vector, pseudoscalar and axial susceptibilities, respectively. In the case of the vector and pseudoscalar channels the one-particle contributions due to $B_c^*$- and $B_c$-mesons are evaluated and subtracted to improve the bounds, obtaining: $\chi_{1^-}^{sub}(m_b^{phys}) = 5.84\,(44) \cdot 10^{-4}$ GeV$^{-2}$ and $\chi_{0^-}^{sub}(m_b^{phys}) = 2.19\,(19) \cdot 10^{-2}$.


Largo Bruno Pontecorvo 3, I-56127 Pisa, Italy
In this work we present the first non-perturbative determination of the hadronic suscep- Consequently, efforts have been devoted to both the experimental and the theoretical investigations of inclusive and exclusive semileptonic decays of hadrons.
In the last few years the exclusive semileptonic B → D ( * ) ν decays have received significant attention. As well known, the crucial ingredients for the extraction of |V cb | are the hadronic form factors entering the exclusive B-meson decays. Reliable information on the latter is given by first-principle calculations performed using QCD simulations on the lattice. However, since the bottom quark is so heavy, its mass in lattice units is still larger than one on the currently available lattice setups. Thus, cutoff effects (as well as large statistical fluctuations) are a limiting factor, which makes very difficult to determine on the lattice the dependence of the form factors on the squared 4-momentum transfer q 2 , where q = p B − p D ( * ) is the lepton-pair 4-momentum, in the full kinematical range m 2 ≤ q 2 ≤ (m B − m D ( * ) ) 2 . As a matter of fact, the results for the B → D ν form factors from lattice QCD [22,23] are available only in a limited range of values of q 2 , namely 9.3 GeV 2 q 2 11.7 GeV 2 , much smaller than the full physical range, 0 q 2 11.7 GeV 2 . Also preliminary results for the momentum dependence of the form factors entering the B → D * ν decays are still restricted at small recoil, i.e. at large values of q 2 [24,25].
In order to supply the lack of information at large recoil, i.e. at small values of q 2 , both experimental and theoretical analyses have to parameterise the form factors in order to describe their momentum dependence in the full kinematical range. In this way, the extraction of |V cb | from experiments may be biased by the theoretical model adopted for fitting the data. In the years most of the analyses used two popular parameterisations, called Boyd-Grinstein-Lebed (BGL) [26][27][28] or Caprini-Lellouch-Neubert (CLN) [29,30]. Recent experimental results [31,32] for the momentum and angular distributions have allowed studies of the role played by the BGL and CLN parameterisations of the semileptonic form factors [33][34][35][36][37][38]. It turned out that the sensitivity to the specific parametrization employed can be solved only by the precise knowledge of the form factors at non-zero recoil.
A possible way to constrain in a model-independent way the q 2 -dependence of the form factors extracted from lattice QCD along the full kinematical range has been proposed recently in Ref. [39]. Starting from the pioneering work of Ref. [40] the dispersive matrix (DM) method of Ref. [39] makes use of unitarity and analyticity bounds [26][27][28] applied to lattice data and introduces a formalism to take into account the uncertainties of the lattice results. The DM method contains several new elements with respect to the original proposal of Ref. [40] and allows to extract, using suitable two-point functions computed non-perturbatively, the form factors at low values of q 2 from those computed explicitly on the lattice at large q 2 , without any assumption about their q 2 -dependence.
The DM method was tested in the case of the semileptonic D → K ν decay, by comparing its results with the explicit direct lattice calculation of the form factors available in the full q 2 -range. It was shown [39] that the DM method is very effective and allows to compute the form factors in a model-independent way and with rather good precision in the low-q 2 region not accessible directly to lattice calculations.
Our aim is to apply the DM method of Ref. [39] to the case of exclusive semileptonic B-meson decays. The DM method requires the non-perturbative evaluation of moments of suitable two-point correlation functions and the subtraction of non-perturbative singleparticle contributions. These are essential ingredients of the whole strategy and they represent per se a complex step, since it involves the determination of single-particle contributions to the two-point functions and, for heavy quarks, the control of large discretization effects. In this work we address the determination of the subtracted susceptibilities in the case of the b → c transitions, while the application of the DM method to the extraction of |V cb | from the experimental data on the exclusive B → D ( * ) ν decays will be the subject of a separate work [41].
The two-point correlation functions are evaluated making use of the gauge ensembles produced by the Extended Twisted Mass Collaboration (ETMC) with N f = 2 + 1 + 1 dynamical quarks at three values of the lattice spacing (a 0.062, 0.082, 0.089 fm) and with pion masses in the range 210 − 450 MeV [42,43]. Our results are extrapolated to the continuum limit and to the physical pion point. In order to reduce cutoff effects we take advantage of the calculation of the two-point correlation functions carried out in perturbation theory both in the continuum and on the lattice in Ref. [39]. For reaching the b-quark point we employ the ETMC ratio method of Ref. [44], that was already applied to determine the physical b-quark mass, m phys b , on the same gauge ensembles in Ref. [45].
In terms of the scalar (0 + ), vector (1 − ), pseudoscalar (0 − ) and axial (1 + ) susceptibilities, which at zero momentum transfer are moments of the appropriate two-point correlation functions, we get the following results for the b → c transition In Refs. [33,35,36] the susceptibilities have been estimated using perturbation theory (PT) at next-to-next-to-leading order (NNLO), obtaining: While the vector and pseudoscalar susceptibilities are only 4% and 7% lower than our findings (2) and (3), the scalar and axial susceptibilities are 20% lower than our results (1) and (4), respectively. In all cases the differences are within ∼ 2.5 standard deviations.
In order to improve the bounds of the DM method, in the case of the vector and pseudoscalar channels the one-particle contributions due to B * c -and B c -mesons are evaluated and subtracted from Eqs. (2) and (3), respectively, obtaining Eqs. (1), (4), (5) and (6)  The plan of the paper is as follows. In Section II we give the basic definitions of the the quantities relevant in this work, namely the longitudinal and transverse susceptibilities, which are moments of suitable two-point correlation functions. In Section III we describe the lattice calculations of the latter ones using the gauge ensembles produced by the ETMC with N f = 2 + 1 + 1 dynamical quarks (see Appendix A) for the h → c transition, where h is a quark heavier than the charm. We illustrate the presence of contact terms in the longitudinal susceptibilities and how to avoid them by means of Ward Identities (WIs). We apply the perturbative calculations of Ref. [39] to understand the main features of the contact terms and also to get a beneficial reduction of the cutoff effects on the susceptibilities. In Section IV we apply the ETMC ratio method of Ref. [44] to reach the physical b-quark point. We first construct suitable ratios of the susceptibilities, evaluated at two subsequent values of the heavy-quark mass m h , and perform the extrapolation to the continuum limit and to the physical pion point. Then, the resulting ratios are smoothly interpolated at the physical b-quark mass, determined in Ref. [45], taking advantage of the fact that the ratios are defined in such a way to guarantee that the their values are exactly known in the heavy-quark limit. In Section V we evaluate the contributions of the B * c -and B c -mesons to the vector and pseudoscalar susceptibilities, respectively, obtaining in this way a more stringent bound for the DM method of Ref. [39]. Finally, Section VI is devoted to our conclusions and outlooks for future developments.

II. EUCLIDEAN TWO-POINT CORRELATION FUNCTIONS
In this Section we briefly recall the definition of the longitudinal and transverse susceptibilities, which are the quantities relevant in this work. More details can be found in Ref. [39].
Let's start with the vacuum polarization tensors Π related to the product of vector and axial-vector currents of the the flavor changing b → c weak transition [28,30].
Performing a Wick rotation from Minkowskian coordinates to Euclidean ones, x ≡ (t, x), the vacuum polarization tensors can be written as [39] Π where γ E µ are (Hermitean) Euclidean Dirac matrices (i.e., γ E 0 = γ M 0 , γ E = −i γ M and γ E 5 = γ M 5 ) and Q ≡ (Q 0 , Q) is an Euclidean 4-momentum (i.e., Q 2 ≡ Q 2 0 + | Q| 2 ). In Eqs. (7)-(8) the quantities Π 0 ± and Π 1 ∓ are called the polarization functions and the subscirpt identifies the spin-parity of the various channels. In particular, the term proportional to Π 0 + (Π 0 − ) represents the longitudinal part of the polarization tensor with vector (axial) four-currents, while the term proportional to Π 1 − (Π 1 + ) is the transverse contribution to the polarization tensor with vector (axial) four-currents. In what follows, for sake of simplicity, we will omit the explicit superscript E in the definition of the Euclidean γ-matrices.
Choosing the four-momentum Q in the temporal direction Q = (Q, 0) one has where the Euclidean correlators C (0 + ,1 − ,0 − ,1 + ) (t) are given by The quantities relevant in this work are the susceptibilities χ(Q 2 ), which are either first or second derivatives of the polarization functions Π(Q 2 ), namely where j 0 (x) = sin(x)/x and j 1 (x) = [sin(x)/x − cos(x)]/x are spherical Bessel functions. Note that the longitudinal derivatives (14) and (16) are dimensionless, while the transverse ones (15) and (17) have the dimension of [E] −2 , where E is an energy.
Eqs. (14)- (17) have been obtained in the Euclidean region Q 2 ≥ 0, but, as shown in Ref. [39], they can be easily generalized also to the case Q 2 < 0. In the Euclidean region Q 2 ≥ 0 a good convergence of the perturbative calculation of the above derivatives is expected to occur far from the kinematical regions where resonances can contribute. In the case of the b → c weak transition this means down to Q 2 = 0 [28,30]. Thus, the value Q 2 = 0 has been generally employed in the evaluation of the dispersive bounds on heavy-to-heavy [28,30,33,35,36] and heavy-to-light [40,46] semileptonic form factors.
On the contrary, with a non-perturbative determination of the two-point correlation functions we can use the most convenient value of Q 2 at disposal, namely the value which will allow the most stringent bounds on the semileptonic form factors [39]. In this work we will limit ourselves to the usual choice Q 2 = 0, which will allow the comparison with the perturbative results at NNLO obtained in Refs. [33,35,36], and we will leave the investigation of the choice Q 2 = 0 to a future work.
At Q 2 = 0 the derivatives of the longitudinal and transverse polarization functions correspond to the second and fourth moments of the longitudinal and transverse Euclidean correlators, i.e.
In Ref. [39] it has been shown that the Ward Identities (WIs), which should be satisfied by the vector and axial-vector quark currents, allow to express the longitudinal susceptibilities (18) and (20) as the fourth moments of the scalar and pseudoscalar correlation functions, namely one has where

III. LONGITUDINAL AND TRANSVERSE SUSCEPTIBILITIES
The gauge ensembles used in this work have been generated by ETMC with N f = 2 + 1 + 1 dynamical quarks, which include in the sea, besides two light mass-degenerate quarks (m u = m d = m ud ), also the strange and the charm quarks with masses close to their physical values [42,43]. The ensembles are the same adopted to determine the up, down, strange and charm quark masses in Ref. [47] and the bottom quark mass in Ref. [45]. Details are given in Appendix A.
Using the ETMC gauge ensembles of Table VII we have evaluated the following twopoint correlation functions where q 1 and q 2 are the two valence quarks with bare masses aµ 1 and aµ 2 given in Table VII, while the multiplicative factor Z γ (γ = {V, A, S, P }) is an appropriate renormalization constant (RC), which will be specified in a while. Indeed, we consider either opposite or equal values for the Wilson parameters r 1 and r 2 of the two valence quarks, namely either the case r 1 = −r 2 or the case r 1 = r 2 . Since our twisted-mass setup is at its maximal twist, in the case . The value of λ, which is the same as the one adopted in Ref. [45], is such that Using the gauge ensemble B25.32 as a representative case, our results for the vector and axial longitudinal susceptibilities are shown in Fig. 1 and for the transverse ones in Fig. 2 at either opposite or equal values of the valence-quark Wilson parameters, which will be denoted hereafter by (r, −r) and (r, r).
The following comments are in order.
• Differences among results corresponding to the two r-combinations are expected to occur because of (twisted-mass) discretization effects, but it is clear that such differences are much larger in the longitudinal channels (particularly for χ 0 + ) with respect to the transverse cases.
• For both r-combinations the transverse susceptibilities increase as m h increases, while the opposite behavior occurs for the longitudinal ones with the exception of • Because of charge conservation, the susceptibility χ 0 + evaluated for m h = m c should vanish in the continuum limit. This seems to hold for the (r, r) combination, while it is strongly violated in the (r, −r) case.
The above observations point toward the presence of extra contributions coming from possible contact terms related to the product of two currents, which appear in all the correlators (26)- (31). The issue of contact terms, which may affect the evaluation of the correlators for any lattice formulation of QCD, has been throughly investigated for our ETMC setups in Refs. [48,49] in the case of the HVP contribution to the muon (g − 2), which, as known, involves the product of two electromagnetic currents (i.e. it corresponds to the degenerate case m h = m c ). The main outcome is that contact terms may not

expressions (22) and (23) derived using WIs
where C S (t) and C P (t) are given by Eqs. (30) and (31) The comparison with Fig. 1 indicates very clearly that, without the effects of the contact terms, the differences among the results corresponding to the two r-combinations are much smaller. More precisely, the contact terms heavily affect the vector longitudinal susceptibility χ 0 + for the combination (r, −r), but they seem to be almost absent for the combination (r, r). Moreover, without contact terms both the vector and axial longitudinal susceptibilities increase as m h increases, as it happens for the transverse ones.

B. Perturbative subtractions of contact terms and discretization effects
Even if the use of the definitions (34)- (35) are free from contact terms, it is very interesting to understand better the behaviour of the effects of the contact terms on the longitudinal susceptibilities (18) and (20). In particular, we want to improve our understanding of the almost total absence of contact terms in the longitudinal vector susceptibility χ 0 + for the combination (r, r) and of its huge quantitative impact for the other combination (r, −r).
To this end, following the procedure discussed in Section VI of Ref. [39], we have calculated the susceptibilities in the free theory on the lattice, i.e. at order O(α 0 s ), using twisted-mass fermions with masses given in lattice units by am h (n) = λ n−1 am phys c for n = 1, 2, ... . These results, which hereafter will be denoted by χ f ree j , contain a physical contribution related to the LO term of PT, i.e. at order O(α 0 s ), and the corrections due to contact terms and discretization effects present in the free theory for our lattice setup, i.e. at all orders O(α 0 s a m ) with m ≥ 0. The former ones, which will be denoted by χ LO j , are known from Ref. [28], namely where u ≡ m c /m h → λ 1−n . Thus, we subtract from our non-perturbative susceptibilities The corresponding results are shown in Fig. 4 for the longitudinal vector and axial susceptibilities and in Fig. 5 for the transverse ones.
The perturbative calculation explains nicely that in the case of the longitudinal vector susceptibility χ 0 + the contact terms are almost absent for the combination (r, r), while they show a huge quantitative effect for the other combination (r, −r). In the latter case the perturbative estimate is not enough to cancel out the contact term for degenerate quark masses m h = m phys c , where χ 0 + is expected to vanish due to current conservation, and orders higher than O(α 0 s ) are clearly required. By comparing the results shown in Fig. 4 with those in Fig. 1 a reduction of the difference between the two r-combinations is visible also in the case of the longitudinal axial susceptibility χ 0 − .
In the case of the transverse susceptibilities the contact terms are absent, but the subtraction of the discretization effects contained in the difference [χ f ree j − χ LO j ] turns out to be beneficial. By comparing the results shown in Fig. 5 with those in Fig. 2 the difference between the combinations (r, r) and (r, −r) is significantly reduced.
An alternative, rather effective, way to get rid of the contact terms in the susceptibility χ 0 + is to subtract the contact terms evaluated non-perturbatively at m h = m c , more precisely by using the formula obtaining in this way that χ 0 + (m h = m c ) = 0 as in the case of the WI-based formula (34).
In Fig. 6 the results 1 obtained using Eq. (41) are compared with those based on the WI formula (34). A reassuring qualitative agreement is obtained, taking into account that discretization effects are different in the two procedures. Note that smaller differences between the two r-combinations occur also in the case of the WI-based formula. Since the subtraction procedure given in Eq. (41) is not applicable to the axial longi-
In Eq. (42) the factor ρ j (m h ) is introduced to guarantee that in the heavy-quark limit m h → ∞ (i.e., n → ∞) one has R j → 1. Using the PT results of Ref. condition is satisfied by  Thanks to the large correlation between the numerator and the denominator in Eq. (42) the statistical uncertainty of the ETMC ratios R j (n; a 2 , m ud ) is much smaller than those of the separate susceptibilities and it may reach the permille level, as it is illustrated in Fig. 8 in the case j = 1 − and n = 5. Moreover, the light-quark mass dependence where r 0 0.47 fm is the value of the Sommer parameter (determined for our lattice setup in Ref. [47]) and R j (n) stands for R j (n; 0, m phys ud ). For sake of simplicity, we have dropped in the notation of the coefficients A 1 and D 1 their dependence on the specific channel j (as well as on the specific value of the heavy-quark mass m h ).
The results obtained with the fitting function (46) are shown in Fig. 8 as an illustrative example. The quality of the fitting procedure may be quite good in several cases, as shown in the left panel of Fig. 8 where the value of χ 2 /(d.o.f.) is significantly less than 1, but it may be also quite poor, as shown in the right panel of Fig. 8 where the value of χ 2 /(d.o.f.) is significantly larger than 1. In the latter cases discretization effects beyond the order O(a 2 ) seem to be required. Moreover, in Eq. (46) the coefficient R j (n) represents the value of the ETMC ratio extrapolated to the physical pion point and to the continuum limit. However, the susceptibilities corresponding to the two combinations and we include a (gaussian) prior, D prior , on the two parameters D 1 and D 2 by increasing the χ 2 variable with the additive term (D 1 /D prior ) 2 + (D 2 /D 2 prior ) 2 . We remind that, for sake of simplicity, we have dropped in the parameters A 1 , D 1 and D 2 their explicit dependence on the specific channel j as well as on the specific value of the heavy-quark mass m h . As for the dependence of the discretization terms on m h , we have found that the fitting values of D 1 and D 2 generally increase as m h increases and they are roughly where ρ j is defined by Eqs. (43)- (44) and 2 In principle, Dprior can be chosen to be different for different values of the heavy-quark mass m h and for different channels j. In practice we have adopted a unique value of Dprior common to all heavy-quark masses and channels. with χ LO j given by Eqs. (36)- (39) and In Eqs. The important feature of the ETMC ratio method is that the extrapolation to the physical b-quark point of the ratios R j can be carried out taking advantage of the fact that by construction lim n→∞ R j (n) = 1 .
However, in order to make smoother the extrapolation to the physical b-quark point (which we remind corresponds to n = 11 for our choice of λ) we consider the double ratios R j (n)/R P T j (n). In this way the nontrivial heavy-quark mass dependence of R P T j is taken into account and the deviations from unity are expected to be further reduced, as it is shown in Figs. 13-14.
Since the PT contribution R P T j is given up to order O(α s ) and the condensate terms start at order 1/(m pole h ) 4 [28], we fit the lattice data for the double ratios adopting the following Ansatz which contains six parameters (A 0 , A 1 , A 2 , A 3 , A 4 and A s 4 ) to be determined by a χ 2minimization procedure. We remind that, for sake of simplicity, we have dropped in the notation of all the parameters their dependence on the specific channel j. Using the double ratios the susceptibilities χ j (m phys b ) can be expressed as and In Eqs. (56)-(57) the products over n include the lattice data up to n = 9 and then the results of the fitting function (55) for n = 10 and n = 11. We remind that the case j = 0 + is treated differently (i.e., n ≥ 3 in Eq.  Table I and compared with the NLO PT predictions based on Eq. (49). It can be seen that non-perturbative effects (as well as PT ones at order O(α 2 s ) and beyond) may affect the vector susceptibility ratios at the level of ∼ 10 − 20%, while the axial ones at the level of few percent only. The values of the susceptibilities extrapolated to the physical pion point (m ud = m phys ud ) and to the continuum limit are collected in Table II and compared with the corresponding PT predictions based on Eq. (49). It can be seen that non-perturbative effects (as well as PT ones at order O(α 2 s ) and beyond) affect significantly all the sus-  Tables I and II the longitudinal and transverse susceptibilities can be calculated at the physical b-quark point using Eqs. (56)- (57). Our findings are collected in Table III  Appendix A) our final results are In Refs. [33,35,36]  While the transverse vector and longitudinal axial susceptibilities are only 4% and 7% lower than our findings (59) and (60), the longitudinal vector and the transverse axial susceptibilities are 20% lower than our results (58) and (61), respectively. In all cases the differences are within ∼ 2.5 standard deviations.

B. Alternative analysis
The results obtained in the previous Section suggests that the heavy-quark mass dependence of the ratios R j is mainly dictated by the corresponding NLO PT predictions R P T j . We want to check further this point and we perform a direct fit of the ratios R j using the following Ansatz which contains six parameters (A 1 , A s 1 , A 2 , A s 2 , A 3 and A s 3 ) to be determined by a χ 2minimization procedure. Note that, at variance with Eq. Using the ratios R j the susceptibilities χ j (m phys b ) at the physical b-quark point can be expressed as and After averaging over all the eight branches of our bootstrap analysis we quote the final values obtained for the susceptibilities χ j (m phys b ), namely which nicely agrees (in a reassuring way) with the findings (58)- (61).
Another check is dictated by the fact that the susceptibilities corresponding to the two combinations (r, −r) and (r, r) of the Wilson r-parameters should differ only by discretization effects starting at order O(a 2 ). This suggests to perform a combined extrapolation to the continuum limit (and to the physical pion point) of both combinations.
In terms of the ratios R In this Section we describe the subtraction of the ground-state contribution.
As well known, at large time distances one has where M j is the mass of ground-state meson H j 12 and Z j is the matrix element Z j ≡ | H j 12 |q 1 Γ j q 2 |0 | 2 with Γ j = {γ 0 , γ, γ 0 γ 5 , γγ 5 , 1, γ 5 } for j = {0 + , 1 − , 0 − , 1 + , S, P }. Thus, the ground-state mass M j and the matrix element Z j can be extracted from the exponential fit given in the r.h.s. of Eq. (71) performed in the temporal region exhibits a plateau. The quality of the plateaux for the effective masses M ef f j (t) is shown in Fig. 18 in two illustrative cases. The quality of the plateaux is good for j = 0 + , 1 − , 0 − and j = P , while it is very poor in the cases j = 1 + and j = S. The latter case is likely to be plagued by the effects of parity breaking (mixing with j = P ) present in our lattice formulation. The same is expected to occur for j = 0 + in spite of the good plateaux  Table IV for the various ETMC ensembles.  Since we make use of the WIs for the longitudinal axial-vector susceptibility (see Eq. (23)), the ground-state contributions χ (gs) 1 − (m 1 , m 2 ) and χ (gs) 0 − (m 1 , m 2 ) are explicitly given by M J/ψ (GeV) 3.13 (5)    . Therefore, we adopt the following Ansatz The quality of the fits is illustrated in Fig. 20 As argued in Ref. [30], the ground-state contributions χ (gs) 0 + and χ (gs) 1 − are expected to be small and they can be conservatively ignored. For the same reason we can ignore also the subtraction of the contributions coming from higher excited bound-states.
Thus, after the subtraction of the ground-state contributions given in the last column of   The ETMC ratio method of Ref. [44] has been adopted to reach the physical b-quark The application of the above results to the extraction of the CKM matrix element |V cb | from the exclusive experimental data will be illustrated in a separate work [41].
Finally, we mention that in this work we have limited ourselves to evaluate the susceptibilities at zero momentum transfer. However, our non-perturbative calculations of the relevant two-point correlation functions allow to consider the most convenient value of the 4-momentum transfer leading to the most stringent bounds on the semileptonic form factors. We leave this investigation to a future work. The ETMC setup adopted in this work is based on the Iwasaki action [58] for the gluons and on the Wilson maximally twisted-mass action [59][60][61] for the sea quarks.
Three values of the inverse bare lattice coupling β and different lattice volumes are considered, as it is shown in Table VII (aµ c ) and heavier-than-charm (aµ h ) regions considered for the 15 ETMC gauge ensembles with N f = 2 + 1 + 1 dynamical quarks (see Ref. [47]). N cf g stands for the number of (uncorrelated) gauge configurations used in this work. In the last column the values of the simulated pion mass, M π , are shown in physical units.
In Ref. [47] the analysis is split into eight branches, which differ in: • the continuum extrapolation adopting for the matching of the lattice scale either the Sommer parameter r 0 or the mass of a fictitious P-meson made up of two valence strange(charm)-like quarks; • the chiral extrapolation performed with fitting functions chosen to be either a polynomial expansion or a Chiral Perturbation Theory (ChPT) Ansatz in the light-quark mass; • the choice between the methods M1 and M2, which differ by O(a 2 ) effects, used to determine the mass RC Z m = 1/Z P in the RI -MOM scheme.
In the present analysis we will make use of the input parameters corresponding to each of the eight branches of Ref. [47].  Tables VIII and IX. Besides the RC Z P we need the RCs of other bilinear quark operators, namely Z V , Z A and Z S related respectively to the vector, axial-vector and scalar currents. They have been evaluated in the Appendix of Ref. [47] in the RI -MOM scheme for Z A and Z S , while for Z V we adopt its determination based on the vector Ward-Takahashi identity.
Throughout this work 5 the results (x i ± σ i ) with i = 1, 2, ...N , obtained within N branches, are averaged according to the following general formula (see Ref. [47]) The second term in the r.h.s. of Eq. (A1), coming from the spread among the results of the different branches, corresponds to the systematic error which accounts for the uncertainties due to the chiral extrapolation, the cutoff effects and the choice of the RC  GeV.