Unitarity Bounds for Semileptonic Decays in Lattice QCD

In this work we discuss in detail the non-perturbative determination of the momentum dependence of the form factors entering in semileptonic decays using unitarity and analyticity constraints. The method contains several new elements with respect to previous proposals and allows to extract, using suitable two-point functions computed non-perturbatively, the form factors at low momentum transfer $q^2$ from those computed explicitly on the lattice at large $q^2$, without any assumption about their $q^2$-dependence. The approach will be very useful for exclusive semileptonic $B$-meson decays, where the direct calculation of the form factors at low $q^2$ is particularly difficult due to large statistical fluctuations and discretisation effects. As a testing ground we apply our approach to the semileptonic $D \to K \ell \nu_\ell$ decay, where we can compare the results of the unitarity approach to the explicit direct lattice calculation of the form factors in the full $q^2$-range. We show that the method is very effective and that it allows to compute the form factors with rather good precision.


I. INTRODUCTION
In this work we present an extended study of two-and three-point lattice correlation functions which are used, together with dispersive techniques [1]- [9], to constrain the lattice predictions for the form factors (FFs) relevant to exclusive semileptonic decays. The form factors are then obtained in a substantially non-perturbative way and without any specific assumption on their momentum dependence. To achieve this goal a complex strategy and several improvements with respect to the original proposal of Ref. [10] are introduced and applied to several calculations. This strategy will be described in detail in the following.
The measurement of the weak charged-current b → c and b → u transitions, more specifically the semileptonic B → D ( * ) ν decays, received an increasing attention in the recent past. The first reason is the precise determination of two of the fundamental parameters of the Standard Model, namely the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements |V cb | and |V ub | [11]. The second reason is the apparent tension between inclusive [12,13] and exclusive determinations of |V cb |, which may be related to other lepton anomalies which have been observed experimentally, see for example [14] and references therein. Among the others let us mention the deviations from lepton-flavour universality (LFU) in the measurements of R D ( * ) [15], the ratios of the branching fractions B → D ( * ) τ ν over B → D ( * ) ν, = e, µ, made by Belle, BABAR and LHCb [16][17][18][19][20][21][22][23][24]. These deviations may be interpreted as a hint of the presence of New Physics (NP) [25][26][27][28][29][30][31][32][33][34]. Anomalous values of the ratios of the branching fractions BR(B → K ( * ) µ + µ − ) to BR(B → K ( * ) e + e − ) have also been interpreted as further signals of a possible violation of LFU, see [35] and references therein. The tension between inclusive and exclusive determinations of |V ub | is even larger and without a satisfactory explanation so far [14,36,37]. From the experimental point of view new data and a better understanding of the experimental systematics could still change the present scenario.
An improvement of the theory, mainly if not uniquely, for exclusive decays is expected from progress in lattice QCD calculations of the relevant form factors. Indeed, most of the reliable information about the form factors relevant in semileptonic B decays is given by first principles calculations made in lattice QCD. One is however limited by the cutoff effects induced by the presence of a quark as heavy as the b-quark in calculations done at a finite lattice spacing a.
Actually, most of the numerical simulations with heavy quarks are performed at a larger than about 0.05 fm, so that for the physical b-quark mass we have m b a 1 and an extrapolation in m b from unphysical values is necessary. In this context discretisation errors affect the value of the form factors at zero recoil and make it difficult to study the momentum dependence of the form factors at large recoil, namely at small q 2 , where q = p B − p D ( * ) is the lepton pair momentum in the decay 1 . For this reason, for example, results for the B → D form factors from lattice QCD [38,39], together with their uncertainties and correlations, are only available in the range 9.3 GeV 2 q 2 11.7 GeV 2 , much smaller than the physical range, 0 q 2 11.7 GeV 2 .
More recently the first, preliminary results for the momentum-dependence of the form factors in B → D * ν decays appeared, but the kinematical region is always restricted at small recoil and not all the form factors are determined with good accuracy [40,41]. There exist also calculations of the form factors relevant for B s → D ( * ) s decays, as well as recent updates of these quantities [42][43][44].
In order to supply the lack of information from explicit calculations of the form factors in the full kinematical range, both the experimental analyses (in order to account for efficiencies and response functions) and the theoretical studies have to assume some parameterisation of the form factors. It is well possible then that the extraction of |V cb | from experiments is biased by the theoretical model adopted in the fits of the data. In the years most of the analyses used two popular parameterisations called Boyd-Grinstein-Lebed (BGL) [4][5][6] or Caprini-Lellouch-Neubert (CLN) [7,8] after the name of the authors.
An important step in constraining in a model independent way the q 2 -dependence of the semileptonic FFs, extracted on the lattice, along the full kinematical range was proposed by L. Lellouch in a pioneering work [10]. This proposal uses the dispersive techniques mentioned before [4][5][6], see also the original proposal in Refs. [1][2][3]9], applied to lattice data and introduces a formalism to take into account the errors of the lattice results. In spite of the use of the form factors derived from first principles in lattice QCD, the proposal of Ref. [10] relies, for the unitarity constraints, on the perturbative calculation of the two-point current correlation functions. To our knowledge, in spite of the large use of the dispersive techniques discussed above, no one has systematically used the lattice two-point correlators computed non perturbatively in numerical simulations to constraint the FFs in semileptonic decays, not even in the original work [10].
In this work we present an extended study of the two-and three-point lattice correlation functions which are used, together with the dispersive techniques, to constrain the lattice predictions for the form factors. This approach will require a detailed understanding of dispersion relations on the lattice, the combination of perturbative and non perturbative calculations, the renormalisation of lattice T-products, the control of lattice artefacts and the peculiar treatment of the statistical 1 In some cases the dependence of the form factors on q 2 in the whole allowed kinematical region is supplemented by using also the results of QCD sum rules calculations at small q 2 .
errors to be discussed in the following. At the end we will be able to constrain FFs in a substantially non-perturbative way.
With respect to the proposal by L. Lellouch [10] and other previous studies, the main novelties in this work are as follows: 1. The non perturbative determination of the relevant two-point current correlation functions on the lattice which are then used to implement the dispersive bounds; 2. The possibility of implementing the constraints from the two-point correlation function computed non-perturbatively also in regions where the perturbative calculations used in previous analyses are non reliable and could not be used; 3. The reduction of lattice artefacts in the two-point correlation functions using fixed-order perturbation theory on the lattice and in the continuum; 4. A quite simpler treatment of the lattice uncertainties with respect to the method proposed in Ref. [10];

5.
A new approach to a realistic estimate of the systematic errors present at small values of q 2 , namely at large momenta of the final meson, both at finite lattice spacing and in the continuum limit, based on the results of Refs. [45,46].
We stress an important feature of the method that will be presented in this work. Consider a set of lattice data for a form factor evaluated at a series of values q 2 j of the squared 4-momentum transfer (j = 1, ..., N ). The data are distributed according to a multivariate distribution with given uncertainties and correlations. The (non-perturbative) unitarity bounds act as a filter by selecting only those combinations of the data that satisfy unitarity and analyticity. Since the twopoint correlation functions provide an extra information, in general a new multivariate distribution is obtained, which may even correspond to a more precise (and differently correlated) data set.
Using the new distribution the method of Ref. [10] reproduces exactly each of the data point when q 2 → q 2 j . It behaves like a fitting procedure passing exactly for the given data set. This is at variance with what may happen adopting the BGL or CLN parameterisations. Indeed, in these cases there is no guarantee that the parameterization reproduces exactly the data set and, therefore, the impact of the unitarity filter may be different.
We have applied the method described in this work to the analysis of the lattice data of the semileptonic D → K ν decays obtained in Ref. [47]. We use this process as a training ground for the dispersive approach to show that starting from a limited set of data at large q 2 it is possible to determine quite precisely the form factors in a model independent way in the full kinematical range, obtaining a remarkable agreement with the direct calculations from Ref. [47]. This finding opens the possibility to obtain non-perturbatively the form factors entering the semileptonic B decays in their full kinematical range. The extension of the non-perturbative dispersive approach described here to B decays requires a remarkable effort in the reduction of discretisation effects due to the large masses and momenta involved and it will be the subject of a forthcoming publication.
The plan of the paper is as follows. The definition of the relevant FFs entering in semileptonic decays is introduced in Section II, where we also give their expressions in the Heavy Quark Effective Theory (HQET) [8]. The definition of the hadronic tensors from the two-point Green functions of suitable bilinear currents and their expression in terms of the form factors is given in Section III, where we also recall the dispersive bounds that can be extracted from these two-point Green functions. The material in Sections II-III can be found in several papers. It is, however, useful to have it collected before its use in the present analysis. In particular, we recall the basic features of the dispersive matrix approach of Ref. [10] and in Appendix A we also provide new analytical expressions for the numerical evaluation of the unitarity bands of the form factors. The Euclidean lattice correlation functions corresponding to the Minkowskian Green functions that we used in our numerical simulation are presented in Section IV. In Section V we discuss the treatment of the statistical and systematic errors with the dispersive bounds and in the presence of kinematical constraints among different form factors. This Section will give us the opportunity of discussing briefly the treatment of the statistical and systematic errors using the Bayesian approach of Refs. [45,46], which is different from the treatment of the errors suggested in the original proposal by L. Lellouch in Ref. [10]. In Section VI we give the results of the calculation of the two-point correlation functions in perturbation theory both in the continuum and on the lattice, while in Section VII 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 [48,49] (see Appendix B). The perturbative calculations are used to reduce discretisation errors and to improve the extrapolation of the numerical non-perturbative results to the continuum limit. In Section VIII the dispersive matrix approach is applied to the D → K form factors obtained in Ref. [47] both in the continuum limit and at finite lattice spacing (using the same gauge ensembles adopted for the evaluation of the two-point correlation functions). Finally, Section IX contains our conclusion and outlooks for future developments.

II. FORM FACTORS IN SEMILEPTONIC DECAYS
In this Section we introduce the relevant FFs for semileptonic decays and convert them in the form used in the HQET, which is particularly suitable for the expansion in inverse powers of the heavy quark mass. For definiteness, since this is our final goal, the formulae will refer to B → D and B → D * decays. With trivial modifications, the same formalism can be applied to D → K and D → K * semileptonic decays, which is the case study for which we present a complete numerical analysis here, as well as also to K → π and other semileptonic decays.
We use the following form factor classification: where the weak currents are given by f i = f i (q 2 ) is the generic hadronic FF that depends only on the squared 4-momentum transfer which is the only non-trivial Lorentz-invariant quantity, and D(B) is a polarization 4-vector. Although the initial state is always aB ( * ) meson containing a b quark, rather than a B ( * ) meson, to simplify the notation in the following we will omit the bar. The normalisation of the states differs from the Feynman one by a factor equal to the square root of the meson mass, |p M Feynman = √ M M |p M , because this is more convenient in the framework of the HQET.
One can define the recoil variable w, given by the scalar product of the meson four-velocities, and a new dimensionless variable r ≡ m D ( * ) /m B ( * ) . In Ref. [8] Caprini, Lellouch and Neubert introduced the quantities h i , which are linear combinations of the FFs f i previously defined, expressed as functions of the recoil variable w rather than of q 2 (the index i labels a generic form factor). The h i describe the decompositions (1)- (2) in terms of the meson 4-velocities instead of the meson 4-momenta, according to the following classification: • scalar FFs • vector FFs • pseudoscalar FFs • axial FFs The last classification of the FFs allows to separate different values of spin and parity quantum numbers from each other. For simplicity, all scalar, pseudoscalar, vector and axial quantities have been presented indicating initial and final states as superscripts, so that we can easily remember to which process each form factor refers to.

III. TWO-POINT CORRELATION FUNCTIONS
The bounds on the different FFs are derived from the two-point functions of suitable currents.
The starting point is the Fourier transform of the T-product of two hadronic currents, which generalizes the definition of the hadronic vacuum polarization (HVP) tensor. Assuming where q is the current 4-momentum transfer and p n the 4-momentum of the intermediate n-particle state. A similar result can be derived for the case x 0 < 0. The completeness sum runs over all possible intermediate hadronic states and in particular we will focus our attention onto either a single-particle B ( * ) c -meson state or two-particle states composed by a B ( * ) -meson and a D ( * ) -meson.
Note that pseudoscalar (scalar) particles will be characterized only by their four-momentum p µ while the vector (axial) mesons will be distinguished also by their polarization four-vector µ . The link between the two-particle states appearing in the completeness sum (9) and the classification of the FFs introduced in Eqs. (1)-(2) is given by the substitution which can be simply realised by inverting the sign of p D ( * ) and by analytic continuation of the FFs Thus, the amplitudes entering in semileptonic decays are strongly correlated to the T-product in Eq. (9), which is the reason why the latter is so important to constrain the form factors.
The Fourier transform of the T-product defines the following HVP tensors: where V µ , A µ are defined in Eq. (3) and the subscripts 0 ± ,1 ± represent spin-parity quantum numbers of the intermediate states. Note that inserting a completeness sum between the vector or axial four-currents we are able to relate these expressions to Eq. (9).
The quantities Π 0 ± ,Π 1 ∓ are called polarization functions. In particular, the term proportional to Π 0 + (Π 0 − ) represents the longitudinal part of the HVP tensor with vector (axial) four-currents, while the term proportional to Π 1 − (Π 1 + ) is the transverse contribution to the HVP tensor with vector (axial) four-currents.

A. Dispersion relations and analytic expressions for exclusive B decays
The imaginary parts of the longitudinal and transverse polarization functions introduced in Eqs. (12) are related to their derivatives with respect to q 2 by the dispersion relations In what follows we will denote by χ a generic susceptibility. From a dimensional point of view note that the longitudinal (scalar/pseudoscalar) susceptibilities χ 0 ± are dimensionless, while the transverse (vector/axial) ones have dimension [E] −2 , where E is an energy.
The two-particle contribution to the polarization functions can be expressed in terms of the FFs defined in the previous Section: • scalar channel • vector channel • pseudoscalar channel • axial channel where w(q 2 ) is given by Eq. (4) and the quantity β i is defined as Note that the above expressions hold for a single two-particle B (  to ImΠ J P ,2p .
We will also need the expression for the one-particle contribution of mesons to the imaginary part of the polarization functions. In the case of the vector B * c and of the pseudoscalar B c mesons the one-particle current matrix elements are defined as These one-particle contributions can be opportunely subtracted from the corresponding susceptibilities through the dispersion relations. Thus, we can define where in general we can consider N pseudoscalar and M vector poles (see, e.g., Ref. [50] for reference values of the masses and the decay constants of such poles). The generalization to possible one-particle states below the annihilation threshold in the scalar and axial-vector channels is straightforward.

B. Dispersive bounds
First of all, we will describe the general ideas behind the dispersive method of Refs.
[1]- [9] thanks to which one can obtain bounds for a generic form factor f (q 2 ). We define and we use the analytic continuation of the amplitudes from the kinematical decay region, where m 2 ≤ q 2 ≤ t − , to the single-meson or pair production region, where m 2 Bc ≤ q 2 or t + ≤ q 2 . By means of the dispersion relations we can use the two-point correlation functions computed in QCD to obtain the constraints on f (q 2 ) and then use analyticity to translate these constraints into the physical FFs relevant to semileptonic decays.
The dispersion relations that we will consider in this work have already been introduced in Eqs. (13). We also recall that the derivatives of the various polarization functions can be determined in perturbative QCD only for values of q 2 that are far from the region of production of resonance states, namely where m b ∼ 4.2 GeV and m c ∼ 1.3 GeV are approximate values of the bottom and charm quark masses, respectively. A possible choice is thus also the value q 2 = 0, which has been widely used in the past, particularly in all the calculations that used the perturbative expression of the susceptibilities χ(q 2 ). 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 form factors.
By inserting a complete set of states with the same quantum numbers of a generic current J we have 2 where dµ(n) is the measure of the phase space for the set of states n. As the completeness sum is semi-positive definite, we can restrict our attention to a subset of hadronic states and thus produce a strict inequality. This consideration allows us to rewrite the dispersion relations for χ(q 2 ) as where f (t) is the generic form factor and W (t) is a computable function that depends on the particular form factor under consideration and is related to phase space factors. These have been given explicitly in Eqs. (14)- (17) for all the possible bilinears.
We can now use analyticity to turn the result (25) into a constraint for the semileptonic region. To achieve this goal, it is necessary that the integrand is analytic below the pair-production threshold t < t + . To this end we define which is real for t s < t + , zero for t = t s and a complex number on the unitary circle for t ≥ t + .
We can remove the poles of the integrand of Eq.(25) by multiplying it by appropriate powers of the z(t, t s )s, as determined by the positions t s of the sub-threshold poles. Each pole has a distinct value of t s , and the product z(t, t s1 ) k 1 z(t, t s2 ) k 2 . . . removes all of them. Hence, we re-express Eq. (25) as where t 0 is an arbitrary point that we will define below and we have introduced the Blaschke factor P (t). The latter is a product of many quantities of the form (26) at the position of the sub-threshold poles (i.e., t s ≤ t + ), and the outer functionΦ(t, t 0 ), which is defined for the vector/axial channel In the last expressionP (t) represents a product of the z(t, t s )'s and z(t, t s )'s that remove the sub-threshold singularities and cuts in the kinematical part W (t). Regarding the choice of t 0 in Eq. (27), there are several possibilities which may be more or less convenient depending on the quantity at hand. We have followed the common lore and, in the remaining of the paper, we used t 0 = t − , so that the allowed kinematic interval [0, t − ] corresponds to the range [z max , 0]. The physical values of z max for different semileptonic decays are given in Table I. In a lattice calculation z max will depend on the values of the initial and final meson masses at which the simulations are performed.
A way to account for the bounds imposed by the susceptibilities on the form factors is provided by the two popular parameterisations BGL [4][5][6] or CLN [7,8]. Coming back to Eq. (27), the quantityΦ(t, t 0 )P (t)f (t) can be expanded in a set of orthonormal functions, proportional to powers of z(t, t 0 ). The consequence of this strategy is that the form factor f (t) in the semileptonic region  can be expressed as where, because of Eq. (27), the coefficients a n have to satisfy the unitarity condition In the case of B → D ( * ) decays, since z max < 0.07, the series present in Eq. (29) are usually truncated after the first two or three terms, introducing only small uncertainties in the theoretical predictions.
In this work we do not adopt the BGL or CLN approach, since there is another alternative formulation [10], which is very convenient for translating the information given by the susceptibility χ(q 2 ) into a bound on the form factors. Indeed, we make the transformation where the kinematical functions φ(z, q 2 ) for the different form factors entering B → D ( * ) decays will be specified in Eq. (40) below. At this point, by introducing an inner product defined as whereḡ(z) is the complex conjugate of the function g(z), the inequality (33) can simply be written where we have also used the positivity of the inner product.
C. Extending the method of the dispersive bounds Following Refs. [9]- [10], we define the function g t (z) as wherez(t) is the complex conjugate of the variable z(t) defined in Eq. (32) and z is the integration variable of Eq. (34). It is then straightforward to show that .
Let us introduce the matrix g tn |φf g tn |g t g tn |g t 1 · · · g tn |g tn In a numerical simulations of lattice QCD the values t 1 , · · · , t n will correspond to the squared 4-momenta at which the FFs have been computed non-perturbatively and that will be used as inputs for constraining the FF in regions non accessible to the calculation.
Note that the first matrix element in (38) is the quantity directly related to the susceptibility χ(q 2 ) through the dispersion relations, see Eq. (35). To be more specific, in the case of B → D decays, in terms of the longitudinal and transverse susceptibilities χ 0 + (q 2 ) and χ 1 − (q 2 ) we have that: where φ 0,+ are kinematical functions where n I is an isospin Clebsch-Gordan factor and When analyticity does not hold, i.e. when a form factor has, for instance, N poles at t = t P 1 , t P 2 , · · · ..., t P N , it is sufficient to modify the kinematical function φ according to and the previous definitions will continue to be valid.
The positivity of the inner products (37) guarantees that the determinant of the matrix (38) is This condition can be rephrased in the second order inequality with ,··· } is the minor obtained by deleting the rows i 1 , i 2 , · · · and the columns j 1 , j 2 , · · · . Calling ∆ the discriminant of the inequality (44), one can show that so that at the end the relevant quantities will only be α, β, ∆ 1 , ∆ 2 . Note that α and ∆ 2 are tindependent, i.e. they are given numbers once the susceptibility χ(q 2 ) and the lattice QCD inputs are chosen. On the contrary, β and ∆ 1 are t-dependent. Moreover, only the quantities β and ∆ 2 depend on the chosen value of q 2 .
At this point, since ∆ 1 ≥ 0 by construction, the inequality (44) will have an acceptable solution only when ∆ 2 ≥ 0. If this condition is satisfied, by expressing the scalar product g t |φf according to Eq. (37) we obtain the following unitarity constraints on the form factor f (t) where Thus, by using a direct lattice measurement of the form factors at the points t 1 , t 2 , . . . , t n and the two-point functions of the suitable currents we can constrain the form factors in regions of momenta which for several reasons may not be accessible to lattice simulations. The application to the case of the semileptonic D → K decays will be presented in Section VIII.
The positivity of the inner product implies that α and ∆ 1 are strictly positive, in such a way that, when also ∆ 2 (q 2 ) ≥ 0, the bounds in Eq. (48) are well defined. We stress that the unitarity filter ∆ 2 (q 2 ) ≥ 0 is t-independent, which implies that, when it is not satisfied, no prediction for f (t) is possible at any value of t.
We point out an interesting feature of the dispersive approach based on the matrix (38). When the momentum transfer t coincides with one of the data points, i.e. when t → t j , the determinant ). In other words the form factor f (t), obtained from the dispersive matrix method, reproduces exactly the given set of data points. This is at variance with what may happen using the BGL or the CLN parameterisations, i.e. when the number of powers of z included in Eq. (29) is truncated below the number of data points. In this case there is no guarantee that the parameterization reproduces exactly the set of input data.
Explicit analytical expressions for f lo(up) (t, q 2 ), which are very useful for their direct numerical evaluation, are given in Appendix A.

IV. EUCLIDEAN TWO-POINT FUNCTIONS
In this Section we discuss in details the approach that has been followed in order to constrain the values of the FFs from the two-point correlation functions computed non perturbatively by numerical QCD simulation on the lattice. Many of the definitions and formulae introduced in this Section have been used on the lattice to compute the HVP function of two electromagnetic currents [51] and its isospin-breaking corrections [52] contributing to the muon g − 2.

A. Basic definitions
We compute the correlation functions at the Euclidean four-momentum Q ≡ (Q 0 , Q), given in terms of the Minkowskian momentum q ≡ (q 0 , q) by the relations Q 0 = iq 0 and Q = q. With this choice Q 2 = −q 2 . Furthermore, we perform a Wick rotation on the coordinates, so that we pass from the Minkoskian coordinates x M = (τ, x) to the Euclidean ones x = (t, x), with t = i τ . The vector and axial HVP tensors take the form where we have introduced the currents defined in terms of Hermitian, Euclidean Dirac matrices satisfying the anticommutation relations In the following we will omit the explicit subscript E in the definition of the Euclidean currents and γ-matrices.

B. Correlators and derivatives of the polarization functions on the lattice
A convenient choice is to work with the momentum Q = (Q 0 , 0) so that where the explicit expressions of the various correlators computed at the time distance t are By recalling the definition of the spherical Bessel functions and given that we get In this work, in view of a comparison with the results obtained by using the perturbative calculation of the susceptibilities, we will take Q 2 = 0. In this case j 0 (0) = 1 and lim x→0 j 1 (x)/x = 1/3, so that the derivatives of the longitudinal and transverse polarization functions are equal to the second and the fourth moments of the longitudinal and transverse Euclidean correlators, respectively. We point out again that, by using the two-point correlation functions determined non-perturbatively, we may constrain the form factors also at Q 2 = 0.

C. Ward Identities
Some relations, which will be particularly useful in the analysis of the numerical results, can be derived using the Ward Identities (WIs) that the vector and axial vector quark currents satisfy where m b and m c are the (bare) masses of the bottom and charm quarks respectively. Hence, by defining two further (Euclidean) polarization functions connected to the scalar and the pseudoscalar currents, namely the WIs imply that from which we obtain Moreover, by performing a double derivative with respect to Q 2 we get At this point, we can define the scalar and pseudoscalar analogues of Eqs. (53): where the scalar and pseudoscalar Euclidean correlators are defined as The WIs offer thus the possibility to express the derivatives of vector longitudinal and axial longitudinal polarization functions in a different way, namely Thus, when setting Q 2 = 0, we can compute the derivatives of the longitudinal vector and axialvector polarization functions directly through the fourth moments of the scalar and pseudoscalar correlators, respectively. The advantage of the above procedure based on the WIs will be clarified in Section VII.
Eqs. (56) and (64) where In general, on the lattice, the WIs are modified by discretisation terms [53] that vanish in the continuum limit, namely when the lattice spacing a → 0. For example, with improved fermion actions like Wilson improved [54,55], Maximal Twisted [56][57][58] and Domain Wall Fermions [59], we have where Z V and Z A are appropriate renormalization constants for the lattice vector and axial-vector currents [53].
Further discretisation effects may enter in the T-products computed on the lattice, that we use in our analysis. We propose to reduce discretisation errors by using a combination of nonperturbative and perturbative subtractions which were found very effective in the past (see later Section VII). Although the treatment of the statistical errors in the presence of kinematical constraints was already discussed in Ref. [10], in this work we introduce a different method which looks to us simpler to implement. This method, together with the skeptical Bayesian approach of Refs. [45,46], allows us also to treat the effects induced by discretisation and other systematic effects.

A. Generation of the bootstrap events
The machinery described in Section III allows us to compute the lower/upper bounds of f 0(+) (t), once we have chosen our set of input data, i.e. {χ 0 + (1 − ) , f 0(+) (t 1 ), · · · , f 0(+) (t n )}. Thus, the input data set is made of 2n + 2 quantities: the n values of the scalar form factor f 0 , the n values of the vector form factor f + and the two susceptibilities χ 0 + and χ 1 − . For sake of simplicity we are considering the same number of data points for both the scalar and the vector form factors evaluated at the same series of values t i (i = 1, · · · , n).
The crucial question is, however, how to propagate the uncertainties related to these quantities into the evaluation of the FFs f 0(+) (t) at a generic value of t. To answer this question, we propose a method different from the one described in Ref. [10]. We start by building up a multivariate Gaussian distribution with mean values and covariance matrix given respectively by The first option is certainly to be preferred to reduce the statistical noise by taking properly into account all the correlations of the data.
For each jackknife/bootstrap event we consider the (n + 1) × (n + 1) matrices M 0 and M + (see Eq. (38)) corresponding to the scalar and vector form factors, respectively. Since ∆ 0 1 = ∆ + 1 is non-negative by definition, the positivity condition (46) implies that both ∆ 0 2 and ∆ + 2 should be positive. Thus, we compute ∆ 0(+) 2 and verify their signs. If either ∆ 0 2 or ∆ + 2 results to be negative, then the event is eliminated from the sample. From the physical point of view, this step can be read as a consistency check between all the input data, namely the susceptibilities and the FFs for that particular bootstrap. At the end of the procedure, we will be left with N boot ≤ N boot events.

B. Implementation of the constraints
At t = 0 the FFs f 0 and f + are subject to the constraint In order to satisfy this condition, in the subset of the N boot events satisfying the unitarity filters ∆ 0(+) 2 ≥ 0, we select only the N * boot ≤ N boot events for which the dispersive bands for f 0 and f + overlap each other at t = 0. This corresponds to impose the conditions where f lo,up (t, q 2 ) were defined in Eq. (48) for a generic form factor f . Omitting for simplicity the argument q 2 at which the susceptibilities χ 0(+) are calculated, the conditions (69) can then be rephrased as As already said, the above condition select N * boot ≤ N boot events. Following Ref. [10] for each of the N * boot events we define We now consider the form factor f (0) to be uniformly distributed in the range given by Eq. (72) and we add it to the input data set as a new point at t n+1 = 0. To be more precise, for each of We then consider two modified (n + 2) × (n + 2) matrices, M 0 C and M + C , that have one more row and one more column with respect to matrices M 0 and M + and contain the common form factor f (t n+1 = 0), namely matrices of the form For any point t at which we want to predict the allowed dispersive band of the form factor f (t) (which can be either f 0 (t) or f + (t)) without directly computing it in our simulation, we compute the matrix M C and using Eq. (48) we get f lo (t) and f up (t). This can be done for each of the N 0 events. Let us indicate the result of the k-th extraction by f k lo (t) and f k up (t), respectively. Then, for each of the N * boot events the lower and upper bounds f lo (t) and f up (t) can be defined as At this point we can generate the bounds of the form factor f (t). To achieve this goal, we combine all the N * boot results f i lo,up (t) (i = 1, · · · , N * boot ) to generate the corresponding histograms and fit them with a Gaussian Ansatz, as it is shown in Fig. 1 in an illustrative case. From these fits we extract the average values f lo(up) (t), the standard deviations σ lo(up) (t) and the corresponding correlation factor ρ lo,up (t) = ρ up,lo (t), namely It is understood that the above procedure is applied for both the scalar f 0 (t) and the vector An important aspect of the full procedure is a good determination of the correlations among the matrix elements in Eq. (73), namely among the corresponding form factors. This is automatically achieved by generating the jackknife/bootstrap events with the procedure A) described in subsection V A, whereas an accurate determination of the correlation matrix of the form factors must be provided when one wants to use data produced by other groups.

C. Combination of the lower and upper bounds for each FF
After the steps described before, for any choice for t we obtain from the bootstrap events (pseudogaussian) distributions for f 0,lo (t), f 0,up (t), f +,lo (t) and f +,up (t) as well as the corresponding mean values, standard deviations and correlations. We combine them according to the following procedure.
Let us consider a single bootstrap event in which f L is the lower bound and f U is the upper one for a generic FF at the given value of t (for sake of simplicity we omit the t-dependence for a while). We associate to the FF f a flat distribution between f L and f U where f = f 0(+) , f U = f 0(+),up , f L = f 0(+),lo and Θ is the Heaviside step function. The mean value and the variance associated to the distribution (78) are respectively given by It is however necessary to average over the whole set of bootstrap events. Since the lower and the upper bounds of a generic FF are strongly correlated, we adopt a multivariate Gaussian distribution to describe them, i.e.
where f lo(up) represents the mean of the lower(upper) bound over all the bootstrap events, given by Eq. (75), and C is the covariance matrix with σ lo(up) and ρ lo,up = ρ up,lo being given by Eqs. (76) and (77), respectively. In Eq. (81) the normalization has been chosen so that Using the product of the distributions (78) and (81) we can compute the final values of the form factor f (t) and its variance σ 2 f (t) as

D. Systematic effects and the skeptical approach
In various analyses of the lattice data that we performed to test our method we have encountered the following phenomenon. For some of the bootstrap events no solution could be found, either because ∆ 2 < 0 or because there is no overlap between the allowed regions for f 0 (0) and f + (0).
This may obviously happen for a statistical fluctuation of the sample at hand. When, however, the fraction of rejected bootstrap events is large, say much larger than 50%, one may argue that systematic effects come into play (e.g. lattice artefacts), that have not been corrected for, which may jeopardise the unitary relations. In cases like these, however, it is still possible to extract the relevant information, which in our case is the allowed interval for the form factors, by using different statistical approaches like the skeptical one discussed in Refs. [45,46].
We stress that the above issue is not important in the case of the D → K semileptonic decays for the analyses of the lattice data either extrapolated to the physical point and to the continuum limit (see Section VIII A) or obtained at finite lattice spacing and unphysical pion masses (see Section VIII B). In these cases we see that practically all the bootstraps generated for the FFs survive after both the unitarity and the kinematical constraints. We anticipate that the skeptical approach, instead, turns out to be useful in the study of the b → c transitions, where we expect much larger discretisation effects which may manifest as apparent violation of the unitarity constraints.
Focusing our attention for instance onto the B → D case, the skeptical approach is useful when the matrix procedure explained in this work is applied to the FFs presented by MILC [39]. Indeed only the ∼ 15% of the initially generated events survive to the simultaneous application of the unitarity and kinematical constraints. On the contrary, the introduction of the skeptical method described below allows the largest part of the generated bootstraps events to pass both filters. This fact will be properly illustrated in a forthcoming paper about the application of the matrix method to the b → c transitions.
In the reminder of this Section we explain the general idea of the skeptical approach.
Probability theory helps us in building up a model in which the values and uncertainties of the physical quantities about which we are in doubt are allowed to vary from the nominal ones.
Obviously, the model is not unique, as not unique are the probability distributions that can be used. The simplest choice consists in enlarging the reported standard deviations σ i of the measured points, by assuming the true standard deviations, σ t i are related to the σ i by a factor r i , one for each of the measured points, σ t i = r i σ i , whereas the average valuesf i is the same. All the points i are treated democratically and fairly, i.e. our prior belief of each r i has expected value equal to one; its prior distribution does not depend on the point; we are skeptical, and hence each r i has a priori a wide range of possibilities described by a probability distribution, with a prior 100% standard uncertainty on r i , i.e. σ(r i )/E[r i ] = 1.
The simplest model is to introduce a Gamma probability distribution for the variable r ≥ 0 The parameters A and B are fixed by imposing that this distribution has both mean value and variance equal to 1. A simple calculation shows that this request corresponds to the choice A = B = 1, i.e. P skept (r) = e −r . At this point, we slightly modify the procedure described in Sections V B-V C, namely we build up a multivariate Gaussian distribution whose covariance matrix now is where the σ's are the uncertainties of the measured points in the numerical simulation and ρ ij the corresponding correlation matrix. In lattice simulations performed at finite lattice spacing one can attempt to obtain the physical results either by extrapolating the lattice quantities to the continuum or by reducing the discretisation effects by a subtraction procedure based on perturbation theory. A combination of the two strategies is indeed the most effective one. Some of the lattice artefacts can also be eliminated non perturbatively using, for example, the WIs of the theory [53]. In this Section we shall deal with the perturbative approach which can be implemented in one-loop (or higher-loops) order by computing for a given quantity, say the polarization function or its derivatives, the corresponding Feynman diagrams, at finite lattice spacing. In the analysis of lattice data we will also discuss the extrapolation of the results to the continuum theory which can be combined with the perturbative subtraction of lattice artefacts discussed in this Section.
Recalling the definitions of Π µν V,A given in Section IV, we analyse the structure of a twisted fermions loop, the graphical representation of which, at lowest order, corresponds to the first we are interested in the expansion of the generic polarization function Π( where one can show that C 1 = 0 with maximally TMF.
In the following we discuss as an example the case of the vector current. In this case the coefficients of the expansion can be derived from the generic form of the lattice polarization tensor using the symmetries of the lattice action [60] Π µν where Π V con (Q 2 ) is the continuum polarization tensor, S 5 and S 6 are defined by the expansion of the lattice action close to the continuum limit S ef f = S 4 + a S 5 + a 2 S 6 + a 2 S 7 + . . . , where S k = d 4 x L k with the terms L k containing linear combinations of fields with mass dimension k, and . . . 0 stands for vacuum expectation value of some combination of operators, either local, e.g. S 6 0 , or non local, S 2 5 0 . We will indicate explicitly the dependence of Π µν (Q 2 ) or of the χ's on the quark masses only when it will be necessary for the discussion of the results.
Luckily enough all the divergent or mass dependent lattice artefacts in the first line of Eq. (89) disappear when we apply the derivative with respect to Q 2 to obtain the susceptibilities χ's, see Eqs. (13). There are however terms of O(a 0 ) which remain and that, in some cases, can even make the longitudinal polarization function different form zero even with degenerate quark masses µ 1 = µ 2 . Besides these terms there are discretisation errors of O(a 2 ) or higher which remain. The strategy to reduce their effect is that widely used in literature, see for example Ref. [61]. In our case it is even simpler since the quantities that we consider, namely the χ's, are finite in perturbation theory. Let us call χ LAT (Q 2 , a) the generic susceptibility computed non perturbatively on the lattice, χ(Q 2 , a) the corresponding susceptibility computed in lattice perturbation theory,χ(Q 2 ) the expression resulting from χ(Q 2 , a) by neglecting all contributions which vanish for a → 0, lim a→0 χ(Q 2 , a) →χ(Q 2 ), χ con (Q 2 ) the susceptibility computed perturbatively in the continuum theory. We introduce the following quantities where ∆ 1χ (Q 2 , a) represents the discretisation errors that we want to subtract, and ∆ 2χ (Q 2 , a) the finite terms which are different in the continuum with respect to the lattice case. For local operators a typical example of these kind of corrections is represented by the current renormalisation . In order to extract the subtracted susceptibility we construct then the combination Thus, up to a certain order in perturbation theory and up to non perturbative effects, χ s (Q 2 , a) reduces to the continuum result without discretisation errors. If one uses the one-loop perturbative calculations, the discretisation error then reduces to O(α s a 2 ).
This procedure can easily extended to higher orders in α s . In higher orders in perturbation theory, for the transverse susceptibilities, it will be also necessary to renormalize the quark masses.
For two-point correlation functions involving scalar or pseudoscalar densities, it will be also necessary to renormalise these bilinear operators in some chosen renormalisation scheme. We assume that we use the same renormalisation scheme, e.g. MS or the on-shell scheme for the quark masses.
In this respect we should also write χ LAT (Q 2 , a, µ) andχ(Q 2 , µ) since these quantities will contain the contribution of the counter terms necessary to renormalise the strong coupling constant or the quark masses.
After the subtraction procedure we expect then to have discretisation errors of O(α n s a 2 ), where n depends on the order at which we have computed perturbatively the current-current correlation functions, or non-perturbative discretisation errors of O(a 2 Λ 2 QCD ). Note that χ s (Q 2 , a) is the quantity to be used at finite lattice spacing for extracting the bounds on the form factors f (t) according to the procedure explained in Section III. A detailed discussion of the counter terms up to two-loop perturbation theory will be given in a future publication where the subtraction procedure will be extended to O(α s a 2 ).

A. An instructive example: the vector current polarization tensor at one loop
In order to illustrate the procedure followed to reduce the discretisation errors, as an instructive example we discuss in details the one-loop perturbative calculation of the vector current polarization tensor and of the corresponding susceptibilities.
At lowest order in perturbation theory, by calling k the internal momentum, we may easily compute the correlator of two local vector currents on the lattice where the integration interval represents the first Brillouin zone. Here G i=1,2 indicates the tree-level Wilson twisted-mass propagator, namely where we have defined on the lattice In order to make the calculation it is convenient to define the dimensionless quantities and express Eq.(94) as Taking into account the change of the integration variables, we have that where P µν V (Qa) is a dimensionless quantity which can only depend on dimensionless quantities (Qa, m 1 a, m 2 a, . . . ). At this point we may obtain the χ's by applying the appropriate derivatives with respect to Q µ to the expression given in Eq. (98). Note that any derivative with respect to Q µ we make to obtain the χ's implies the appearance of a factor a in front of the r.h.s. of Eq. When we want to obtain the continuum expression (at this order we do not need to define the renormalisation scheme since everything is finite) it is enough to take the limit a → 0 in the integrand (93) and apply to P µν V (Qa) the derivatives with respect to Q 0 .
It is useful to start by discussing the case of the susceptibilities at Q = 0. In this case, in the continuum, one obtains where N c is the number of colours and the quantities on the l.h.s. being dimensionless can only depend on the ratio u ≡ m 1 /m 2 . In what follows m 2 will always denote the heavier of the two valence quarks in the decaying meson, namely the b quark for B → D ( * ) decays, the charm for D → K ( * ) decays and the strange for kaon semileptonic decays. Note that in the limit m 1 → m 2 (i.e. u → 1) the longitudinal susceptibility χ 0 + (Q 2 ) vanishes because the currents are conserved in this limit. Also on the lattice, as a → 0, the χ's can only depend on u and thus, in perturbation theory we expect where the quantities δχ ( ) i can be eliminated in perturbation theory following the scheme described in Eq. (92). The lattice susceptibilities m 2 2 χ 1 − (Q 2 = 0, a) LAT and χ 0 + (Q 2 = 0, a) LAT are obtained by applying the appropriate derivatives with respect to Q 0 to the expression in Eq. (98) and putting Q 0 = 0. The four dimensional integral can be performed numerically without difficulties.
Since we are able to compute the polarization tensor non perturbatively, in principle we are also able to enforce the unitarity constraints on the FFs at Q 2 = 0. Thus in order to reduce the lattice artefacts we also need the lattice and continuum perturbative calculation in this case.
At one loop, in the continuum, for Q 2 = 0 we may construct two dimensionless quantities, namely u as before and z = Q 2 /m 2 2 . In terms of these variables we obtain [8] where Λ ± = (1 ± u) 2 − z 2 . Also in this case χ 0 + (Q 2 = z 2 m 2 2 ) vanishes for m 1 → m 2 (u → 1). The corresponding lattice quantities are easily obtained as before by applying the appropriate derivatives with respect to Q 0 to the expression in Eq. (98) and putting Q 0 = Q 2 = 0.
On the lattice there is another equivalent way to compute the correlation function in perturbation theory, which is more suited for perturbative calculations at higher orders (more loops). Let us continue the discussion with the example of the vector current. We define the free propagator of the quark (and in higher orders also of the gluon propagator) by numerical Fourier transform of the momentum-space propagator given in Eq. (97) where x ≡ (t, x), y ≡ (t 0 , y) and we have elected one of the lattice directions to our Euclidean time.
At one loop we define It is straightforward to show the relations where the quantities C 0 + (t) and C 1 − (t) are those defined in Eq. (53). We can then obtain the susceptibilities using the expressions in Eq. (56).
The difference between this latter way of computing the χ's, that we will call the x-space method, and the calculation of the χ's from the derivatives applied to Eq. (98), that we denote as the Q-space method, is the following. Eq. (98) refers to a discretised, infinite volume lattice. The Fourier transform in Eq. (103), however, must be done, in practice, on a finite lattice of volume L 3 and time extent T Thus, in order to compare the results with the two different methods we must extrapolate the results obtained from G µν V (t) on a finite lattice to the infinite volume limit. We have done the calculation in the two ways for different values of Q and found very good agreement between them.
A further advantage of the x-space method is that we can perform the calculation on the same finite volume used for the non perturbative calculation of the two-point function.
In order to further reduce lattice artefacts we may extend the subtraction procedure of Eq.   Table VIII, we have evaluated the following two-point correlation where q 1 and q 2 are the two valence quarks with bare masses aµ 1 and aµ 2 given in Table VIII, while the multiplicative factors Z Γ (Γ = {V, A, S, P }) are the appropriate RC of the bilinear currents, which will be specified in a while. 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 r 1 = −r 2 we have  Table IX. The results for the vector and axial longitudinal susceptibilities are shown in Fig. 3 and those for the transverse ones in Fig. 4 at either opposite or equal values of the valence-quark Wilson parameters, which will be denoted hereafter by (r, −r) and (r, r). Differences among the re- sults corresponding to the two r-combinations are expected to occur because of (twisted-mass) discretization effects, but it can clearly be seen that such differences are much larger in the longi- should vanish in the continuum limit. This is strongly violated for the (r, −r) combination, as shown in Table II Table VIII) correspond to nearly the same pion mass and differ in the values of the lattice spacing.
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 (107).
The issue of contact terms, which may affect the evaluation of the correlators for any lattice formulation of QCD, has been throughoutly investigated for our ETMC setup in Refs. [52,60] in the case of the HVP contribution to the muon (g − 2), which as known involves the product of two electromagnetic currents (i.e. the degenerate case m 1 = m 2 ). The presence of contact terms is also evident from the explicit calculation in lattice perturbation theory of Section VI. A strategy to subtract non perturbatively the largest contamination due to these contact terms will be discussed in the following.
The main outcome is that contact terms may not vanish in the continuum limit due to the mixing of the product of two currents with terms proportional to second derivatives of the Dirac delta function. Therefore, a quick inspection of Eqs. (107) reveals that the longitudinal susceptibilities χ 0 + (0) and χ 0 − (0) are affected by contact terms (being second moments), while the transverse ones χ 1 − (0) and χ 1 + (0) are not (being fourth moments). A way to avoid contact terms is to replace the longitudinal susceptibilities in Eqs. (107) with the corresponding expressions (66) derived using WI identities at q 2 = 0, namely where C S (t) and C P (t) are given by Eqs. (107). In this way the longitudinal susceptibilities are evaluated using the fourth moments of the scalar and pseudoscalar correlators and they become free from contact terms. Note that in the degenerate case m 1 = m 2 the susceptibility χ 0 + (0) given by Eq. (108) vanishes at any finite value of the lattice spacing.
The results for the contact-free longitudinal susceptibilities χ 0 + (0) and χ 0 − (0) are shown in Even if Eqs. (108) are essentially free from contact terms, it is, however, very interesting to investigate the impact of the contact terms on the longitudinal susceptibilities in Eqs. (107). In particular, from the observation that χ 0 + (0) is much smaller for the combination (r, r), see Table II, we deduce that in this case the contact terms are much smaller than for the other combination (r, −r).
In general a sizeable reduction of the contact terms can be achieved by using the subtraction procedure developed in Section VI. The perturbative contribution, evaluated at order O(α 0 s ), can already explain the large impact of the contact terms for the combination (r, −r). An even more effective cancellation of the contact terms is obtained by using the non perturbative subtraction proposed below. A detailed, general presentation of the numerical implementation of the subtraction procedure will not be given here. We will discuss it in details in a forthcoming study of the b → c transition, for which also discretisation effects are much larger than in the case of the c → s( ) decays considered in the present work. We find that the subtraction turns out to be beneficial also for reducing the discretization effects, that for b → c are much larger, for both r-combinations.
Here we anticipate that the perturbative contribution to the susceptibility χ 0 + (0), evaluated at order O(α 0 s ) in the degenerate case m 1 = m 2 = m phys c , turns out to vanish for the combination (r, r), while in the case (r, −r) it is equal to ≈ 0.016 for the three ensembles A30.32, B25.32 and D20.48, see Appendix B, which represents ≈ 60% of the corresponding non-perturbative value (≈ 0.026) shown in Table II.
An alternative, rather effective, way to get rid of the contact terms in the susceptibility χ 0 + (0) is to subtract the contact terms evaluated non perturbatively at m 1 = m 2 , more precisely by using the formula obtaining in this way that χ 0 + (0; m 1 = m 2 ) = 0 as in the case of the WI-based formula (108).
In Fig. 6 the results obtained using the non-perturbative subtraction in Eq.
where, for sake of simplicity, χ j (0) stands for χ j (0; m phys ud , 0) and we have not written explicitly dependence of the coefficients A 1 and D 1 on the specific channel j (j = {0 + , 1 − , 0 − , 1 + }). The qual- ity of the various fits is always very good (χ 2 /(d.o.f.) 0.6) and our findings for the extrapolated quantities χ j (0) are collected in Table III.
Averaging over all the eight branches [66] of our bootstrap analysis, see also Appendix B, and over the results corresponding separately to the two r-combinations our final results for the longitudinal and transverse susceptibilities relevant for the c → s transition are χ 0 + (0) = 9.29 (64) · 10 −3 , to improve the dispersive bounds on the semileptonic D → K(K * ) ν decays. In this Section we describe such a subtraction.
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. (112) performed in the temporal region t = [t min , t max ], where the effective mass M ef f At large time distances the quality of the plateaux is good for j = 1 − and j = 0 − , while it is definitely poor in the cases j = 0 + and j = 1 + . The latter ones are likely to be plagued by the effects of parity breaking (mixing with j = 0 − and j = 1 − ) present in our lattice formulation.
We stress that we do not use M ef f j (t) to extract the ground-state masses M j , but we perform the where, for sake of simplicity, M j stands for M j (m phys ud , 0). Our results for M j are collected in Table V.     Table V).
We can now proceed to the evaluation of the ground-state contributions to the susceptibilities, χ (gs) j (0), using the results of the exponential fit (112), and to their subtraction from the global susceptibilities χ j (0), namely By repeating the analysis made in Section VII A for the global susceptibilities χ j (0) in the case of the subtracted ones defined in Eq. (116), we obtain for the χ As a check of our subtraction procedure, we consider the ground-state contribution to χ  [47]. In the next Sections we consider two analyses. In the first one we adopt the results of Ref. [47] already extrapolated to the physical point and to the continuum limit. In the second analysis we make use directly of the results obtained at finite values of the lattice spacing and for unphysical pion masses.

A. Extraction of the form factors in the continuum
The values of f 0 (q 2 i ) and f + (q 2 i ) with the corresponding uncertainties, extrapolated at several values of q 2 i to the physical point and to the continuum limit in Ref. [47], are given in the second and third columns of Table VI.
As already stated in the introduction one of the main reasons to implement the dispersive bounds is that in the case of B decays (B → D ( * ) and B → π), in most of present lattice simulations, one is able to access only the kinematical region close to q 2 max where the final meson has a small momentum. The dispersive bounds instead, following the method discussed in Sections III and IV, allow us to determine with good accuracy also unaccessible points, close to the minimal q 2 , i.e. q 2 ∼ m 2 ∼ 0, without assuming any specific functional form for the q 2 -dependence of the FFs. Bearing in mind possible differences and further difficulties, in order to check the validity of this strategy, we have used only the lattice QCD points for D → K decays at high-q 2 (in boldface in the Table VI). We have implemented our unitarity method and then compared our results with the direct lattice computations in the region at smaller q 2 (see the first five rows of Table VI).
Using the average values and errors given in Ref. [47] we have first applied the procedure B) given in the Subsection V A. Using 20000 bootstrap events we observe that only a fraction of the events equal to ∼ 15% is rejected by the ∆ 0(+) 2 ≥ 0 condition (see Eq. (46)), while almost all the survived ones respect the kinematical constraint f + (0) = f 0 (0) at q 2 = 0. In Fig. 10 we show our results. The orange band represents the extrapolated vector FF f + (q 2 ), while the cyan one corresponds to the scalar FF f 0 (q 2 ). The three red points at the largest values of q 2 , for each form factor, are the bolded values in Table VI used as inputs for our method, while the green and blue ones are the unbolded values shown for comparison.
Our unitarity results shown in Table VI and Fig. 10 exhibit an excellent agreement with the lattice QCD values from Ref. [47] and also a quite similar precision in the whole kinematical range.
In particular, our value of the form factors at q 2 = 0 is which agrees very nicely with the corresponding result 0.765(31) from Ref. [47] , having, we stress, a comparable error even if only three points at large q 2 have been used as input.

B. Extraction of the form factors from each ensemble
Since we have access to the original data of Ref. [47], we can redo the analysis of that paper having computed in Section VII A on the same ensembles also the susceptibilities. The goal is now to implement the matrix method directly on the lattice data points bootstrap by bootstrap and make a totally model-independent extraction of the form factors following the procedure A) of Subsection V A.
As in Ref. [47], our analysis is divided in eight branches which differ by the choice of the scaling variable, the fitting procedures and the choice of the method used to determine non perturbatively the values of the mass renormalization constant (see Appendix B). For every branch we generate 100 bootstrap events in order to take into account the statistical uncertainties. We keep separate all the branches until the continuum and physical pion point is reached. There we combine the bootstrap events and also the branches using Eq. (28) of Ref. [66]. As extensively described in Ref. [47], the lattice data obtained using the three-point correlation functions with the insertion of vector and scalar densities are affected by non-negligible hypercubic effects, which break Lorentz symmetry. The latter ones can be sensibly reduced following the strategy of Ref. [47]. An example of this procedure is shown in Fig. 11. To implement this procedure, however, some functional form for subtracting hypercubic effects, and consequently some model dependence, was introduced.
For the D → K decay the lattice data of Ref. [47] (already interpolated to the physical charm and strange quark masses) cover indeed all the kinematical region in q 2 . The idea is to mimic what happens in lattice calculations of B decays where all the lattice data are concentrated at q 2 ∼ q 2 max . Thus, we have chosen to use, for each form factor, only two points at large values of q 2 corresponding to the D-meson at rest, shown as red markers in Fig. 12. The great advantage of The form factors are already interpolated to the physical charm and strange quark masses.
studying the D → K decay is that we can compare our results obtained with the unitarity procedure to the ones obtained from a direct calculation of the form factors. We applied the matrix method described in the previous Sections to the determination of the FFS using 31 bins in q 2 in the range [−0.5GeV 2 , q 2 max ]. The susceptibilities χ 0 + ,1 − are those computed non perturbatively for each ensemble in Section VII. They have been obtained by eliminating the one particle state both for χ 0 + and χ 1 − . Thus, the kinematical functions φ 0(+) have been modified accordingly to Eq. (42) by including respectively the D * s and the D * 0 (2400) poles. Their masses have been calculated on the same ensembles in Section VII B (see Fig. 9).
To illustrate the procedure we show in Fig. 12 the comparison between our predictions for the  Ref. [47] are interpolated to the physical values of the charm and strange quark masses determined in Ref. [66].
allowed bands of the form factors, obtained by using as inputs only the points denoted as red markers at large q 2 , and the rest of the lattice points that are not used as input in our analysis in the case of the ETMC ensembles B25.32 and D30.48 (see Appendix B). The agreement is excellent.
These results suggest that it will be possible to obtain quite precise determinations of the form factors for B decays by combining form factors at large q 2 with the non perturbative calculation of the susceptibilities.
We now combine the results for all the ensembles and perform the extrapolation to the continuum limit and to the physical pion point adopting the following ansatz where being m the renormalized light-quark mass, and B and F the SU (2) low-energy constants entering the chiral Lagrangian at leading order, whose values were determined in Ref. [66]. In the fitting procedure the parameters c 0 , c 1 , c 2 , c 3 and c 4 , which depend on q 2 and on the form factor, are treated as free independent parameters for each bin in q 2 , and the correlations among the different bins are automatically taken into account by generating events within the jackknife/bootstrap procedure. Differently, the parameter A K is the coefficient of the chiral-log. We have checked that it can be safely fixed at the value A K = 1/2 predicted by the hard-pion SU (2) chiral perturbation theory at q 2 = 0 (see Ref. [47]). For each value of q 2 the quality of the fit (121) with a total of 5 free parameters turns out to be quite good being χ 2 /d.o.f. ∼ 1. We stress that no assumption has been made concerning the q 2 -dependence of the parameters appearing in Eq. (121).
At this point, we recombine the bootstrap events and the branches of the analysis to obtain the final results. In Fig. 13 we present the final bands for the vector and scalar form factors, extrapolated to the physical value of the pion mass and to the continuum limit. The bands agree with the results of Ref. [47] and exhibit a good precision. This demonstrates that the dispersive matrix method allows to determine the semileptonic form factors in their whole kinematical range with a quality comparable to the one obtained by the direct calculations, even if only a quite limited number of input lattice data for each FF (and the non-perturbative susceptibilities) are used 3 .
In Table VII we provide explicitly our final results for the vector and scalar form factors, computed at the eight values of q 2 adopted in Ref. [47], including their total uncertainties. The latter  Table VII are consistent within the uncertainties with the lattice data of Ref. [47].
To conclude we recall the advantages of the present method. The first point is that the determination of the form factors at values of q 2 where there isn't a direct lattice calculation does not assume any functional dependence of the FF on the momentum transfer. Indeed, the analysis at each bin in q 2 is independent of the others. In addition, the results obtained by using only two points in q 2 for each FF and the susceptibilities are of comparable precision with the direct lattice calculations in the full physical range of values of q 2 (see Fig. 12). We are confident that this will remain true for B → D ( * ) or B → π decays, where it is much harder, if not impossible, to compute reliably the FFs at small q 2 . Last but not least, in this analysis we used non perturbative susceptibilities. In the future, we will investigate the use of susceptibilities at non zero momentum.   of q 2 adopted in Ref. [47]. The errors correspond to the sum in quadrature of the uncertainties related to statistical, chiral extrapolation and discretization effects (see text). For comparison the results of Refs. [47] are shown in the second and fourth columns.

IX. CONCLUSIONS
In this work we have presented an extended study of two-and three-point correlation functions on the lattice, that together with known dispersive techniques [1]- [10] allows to constrain the lattice predictions for the form factors relevant to exclusive semileptonic decays. The constraints on the form factors have been implemented by using two-point functions computed in a non-perturbative way. Contrary to the perturbative calculation of the two-point function, this approach will allow in the future to use the unitarity constraints also at non zero momentum. We also introduced a straightforward, and simple to implement, treatment of the uncertainties. This includes the cases where kinematical constraints between form factors are present.
We have then applied the new method to the analysis of the lattice data of the semileptonic D → K decays obtained in Ref. [47]. We have used this example as a training ground for the method and we have shown that it is possible determine the form factors, in a model-independent way, in the region at low q 2 not accessible directly to lattice calculations, as it is the case of exclusive semileptonic B-meson decays. This was achieved by comparing the results of the method with the direct calculation of the form factors along the full kinematical range and allowed us to test the validity of the approach. The application to exclusive semileptonic B-meson decays will be presented elsewhere.
A simple evaluation by induction shows that where, in the case N = 1 it is understood that N i<j=1 (...) → 1. The matrix of which we want to calculate the determinant is given by Eq. (38) of Section III C, namely it has the form where χ is the susceptibility that bounds the inner product φf |φf and, we remind, φ i f i corresponds to the scalar product φf |g t i for the known values of the form factor f i = f (z i ), whereas φf is the scalar product φf |g t of the form factor f (z(t)) that we want to constrain. In order to use a compact notation let us indicate the values of the conformal variable z and of φ(z)f (z) as z 0 and φ 0 f 0 , respectively, so that in what follow the index i may run from 0 to N .
A simple evaluation by induction, as before, yields det[M] = G (N +1) (z 0 , z 1 , z 2 , ...z N ) The unitarity bounds for the (unknown) form factor f 0 result from the condition which implies 4 where (after some algebraic manipulations) α ≡ G (N ) (z 1 , z 2 , ...z N ) ≥ 0 , (A8) Unitarity is satisfied only when γ ≥ 0, which implies χ ≥ χ. Note that d 0 and φ 0 depend on z 0 , while the quantities d j and φ j with j = 1, 2, ...N do not. Thus, the values of β and γ depend on z 0 , while the value of χ does not depend on z 0 and it depends only on the set of input data.
Consequently, the unitarity condition χ ≥ χ does not depend on z 0 .
• Since in terms of the squared 4-momentum transfer t the variable z 0 is given by the annihilation threshold t = t + corresponds to z 0 = −1, while t → −∞ corresponds to z 0 = 1. From Eq. (A10) it follows that unitarity may have no predictive power (i.e. γ → ∞) both at the annihilation threshold t + and for t → −∞.

Appendix B: Simulation details
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 [48,49].
The ensembles are the same adopted to determine the up, down, strange and charm quark masses in Ref. [66] and the bottom quark mass in Ref. [67].
In the ETMC setup the Iwasaki action [68] for the gluons and the Wilson maximally twistedmass action [56][57][58] for the sea quarks are employed. Three values of the inverse bare lattice coupling β and different lattice volumes are considered, as it is shown in Table VIII, where the number of configurations analyzed (N cf g ) corresponds to a separation of 20 trajectories.
At each lattice spacing different values of the light sea quark mass are considered, and the light valence and sea quark masses are always taken to be degenerate, i.e. m sea ud = m val ud = m ud . In order to avoid the mixing of strange and charm quarks in the valence sector we adopt a nonunitary set up in which the valence strange and charm quarks are regularized as Osterwalder-Seiler fermions [69], while the valence up and down quarks have the same action of the sea. Working at maximal twist such a setup guarantees an automatic O(a)-improvement [58,70]. Quark masses are renormalized through the RC Z m = 1/Z P , computed non-perturbatively using the RI -MOM scheme (see Ref. [66]).
The lattice scale is determined using the experimental value of f π + so that the values of the lat-  regions considered for the 15 ETMC gauge ensembles with N f = 2 + 1 + 1 dynamical quarks (see Ref. [66]).
N cf g stands for the number of (uncorrelated) gauge configurations used in this work.
In Ref. [66] eight branches of the analysis were considered. They 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. [66]. The central values and the errors of the input parameters, evaluated using Besides the RC Z P we need the RC's 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. [66] in the RI -MOM scheme for Z A and Z S , while for Z V we adopt its determination based on the vector WI identity.