Heavy-to-light scalar form factors from Muskhelishvili-Omn\`es dispersion relations

By solving the Muskhelishvili-Omn\`es integral equations, the scalar form factors of the semileptonic heavy meson decays $D\to\pi \bar \ell \nu_\ell$, $D\to \bar{K} \bar \ell \nu_\ell$, $\bar{B}\to \pi \ell \bar\nu_\ell$ and $\bar{B}_s\to K \ell \bar\nu_\ell$ are simultaneously studied. As input, we employ unitarized heavy meson-Goldstone boson chiral coupled-channel amplitudes for the energy regions not far from thresholds, while, at high energies, adequate asymptotic conditions are imposed. The scalar form factors are expressed in terms of Omn\`es matrices multiplied by vector polynomials, which contain some undetermined dispersive subtraction constants. We make use of heavy quark and chiral symmetries to constrain these constants, which are fitted to lattice QCD results both in the charm and the bottom sectors, and in this latter sector to the light-cone sum rule predictions close to $q^2=0$ as well. We find a good simultaneous description of the scalar form factors for the four semileptonic decay reactions. From this combined fit, and taking advantage that scalar and vector form factors are equal at $q^2=0$, we obtain $|V_{cd}|=0.244\pm 0.022$, $|V_{cs}|=0.945\pm 0.041$ and $|V_{ub}|=(4.3\pm 0.7)\times10^{-3}$ for the involved Cabibbo-Kobayashi-Maskawa (CKM) matrix elements. In addition, we predict the following vector form factors at $q^2=0$: $|f_+^{D\to\eta}(0)|=0.01\pm 0.05$, $|f_+^{D_s\to K}(0)|=0.50 \pm 0.08$, $|f_+^{D_s\to\eta}(0)|=0.73\pm 0.03$ and $|f_+^{\bar{B}\to\eta}(0)|=0.82 \pm 0.08$, which might serve as alternatives to determine the CKM elements when experimental measurements of the corresponding differential decay rates become available. Finally, we predict the different form factors above the $q^2-$regions accessible in the semileptonic decays, up to moderate energies amenable to be described using the unitarized coupled-channel chiral approach.

By solving the Muskhelishvili-Omnès integral equations, the scalar form factors of the semileptonic heavy meson decays D → π¯ ν , D →K¯ ν ,B → π ν andBs → K ν are simultaneously studied. As input, we employ unitarized heavy meson-Goldstone boson chiral coupled-channel amplitudes for the energy regions not far from thresholds, while, at high energies, adequate asymptotic conditions are imposed. The scalar form factors are expressed in terms of Omnès matrices multiplied by vector polynomials, which contain some undetermined dispersive subtraction constants. We make use of heavy quark and chiral symmetries to constrain these constants, which are fitted to lattice QCD results both in the charm and the bottom sectors, and in this latter sector to the light-cone sum rule predictions close to q 2 = 0 as well. We find a good simultaneous description of the scalar form factors for the four semileptonic decay reactions. From this combined fit, and taking advantage that scalar and vector form factors are equal at q 2 = 0, we obtain |V cd | = 0.244±0.022, |Vcs| = 0.945±0.041 and |V ub | = (4.3 ± 0.7) × 10 −3 for the involved Cabibbo-Kobayashi-Maskawa (CKM) matrix elements. In addition, we predict the following vector form factors at q 2 = 0: |f D→η + (0)| = 0.01 ± 0.05, |f Ds→K + (0)| = 0.50±0.08, |f Ds→η + (0)| = 0.73±0.03 and |fB →η + (0)| = 0.82±0.08, which might serve as alternatives to determine the CKM elements when experimental measurements of the corresponding differential decay rates become available. Finally, we predict the different form factors above the q 2 −regions accessible in the semileptonic decays, up to moderate energies amenable to be described using the unitarized coupled-channel chiral approach. Exclusive semileptonic decays play a prominent role in the precise determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, which are particularly important to test the standard model (SM)-any violation of the unitarity of the CKM matrix would reveal new physics beyond the SM (see for instance the review on the CKM mixing parameters by the Particle Data Group (PDG) [1]). Experimental and theoretical efforts have been devoted to multitude of inclusive and exclusive semileptonic decays driven by electroweak charge currents. For instance, the K 3 decays and those of the type H → φ¯ ν and H → φ ν (hereafter denoted by H 3 or H → φ), where H ∈ {D,B} is an open heavy-flavor pseudoscalar meson and φ ∈ {π, K,K, η} denotes one of the Goldstone bosons due to the spontaneous breaking of the approximate chiral symmetry of Quantum Chromodynamics (QCD), are important in the extraction of some of the CKM matrix elements. Experimentally, significant progresses have been achieved and absolute decay branching fractions and differential decay rates have been accurately measured [2][3][4][5][6][7][8][9][10][11]. On the theoretical side, determinations of the form factors in the vicinity of q 2 = 0 (with q 2 the invariant mass of the outgoing lepton pair) using light-cone sum rules (LCSR) have significantly improved their precision [12,13], and have reached the level of two-loop accuracy [14]. Meanwhile, improvements have been made by using better actions in lattice QCD (LQCD), which have allowed to extract CKM matrix elements with significantly reduced statistical and systematical uncertainties [15][16][17][18][19][20][21]. As a result of this activity in the past decade, lattice calculations on the scalar form factors in heavy-to-light semileptonic transitions have been also reported by the different groups (see the informative review by the Flavour Lattice Averaging Group (FLAG) [22]).
The extraction of the CKM mixing parameters from K 3 and/or H 3 decays relies on the knowledge of the vector [f + (q 2 )] and scalar [f 0 (q 2 )] hadronic form factors that determine the matrix elements of the charged current between the initial and final hadron states 1 . Various parameterizations, such as the Isgur-Scora-Grinstein-Wise updated model [23] or the series expansion proposed in Ref. [24], are extensively used in LQCD and experimental studies. In this work, we will study the scalar form factors in H 3 decays by using the Muskhelishvili-Omnès (MO) formalism, which is a model independent approach to account for Hφ coupled-channel re-scattering effects. The coupled-channel MO formalism has been extensively applied to the scalar ππ, πK and πη form factors, see, e.g., Refs [25][26][27][28][29]. It builds up an elegant bridge to connect the form factors with the corresponding S-wave scattering amplitudes via dispersion relations. The construction of those equations is rigorous in the sense that the fundamental principles, such as unitarity and analyticity, and the proper QCD asymptotic behaviour are implemented. The first attempts to extend this method to the investigation of the scalar H → φ form factors were made in Refs. [30][31][32][33][34][35][36][37][38], but just for the single-channel case. A similar dispersive MO approach has been also employed to study the semileptonicB → ρlν l [37,[39][40][41] andB s →K * lν l [37] decays and the possible extraction of the CKM element |V ub | from data on the four-bodȳ B → ππlν l andB s →Kπlν l decay-modes.
The study of heavy-light form factors using the MO representation incorporating coupled-channel effects has not been undertaken yet. This is mainly because of the poor knowledge on the Hφ interactions up to very recent years. However, a few intriguing positive-parity charmed mesons, like the D * s0 (2317), have been recently discovered [1], giving support to a new paradigm for heavylight meson spectroscopy [42] that questions their traditional qq constituent quark model interpretation. Hence, the study of the Hφ interactions aiming at understanding the dynamics of these newly observed states has become an interesting subject by itself, see, e.g., Refs. [43][44][45][46][47][48][49][50][51][52][53][54] for phenomenological studies and [55][56][57][58][59][60][61] for LQCD calculations. For the D 3 decays, several LQCD results on the relevant form factors have been recently reported, see, e.g., Refs. [16,17,21]. This situation makes timely the study of the scalar D → φ form factors by means of the MO representation incorporating our current knowledge of Dφ interactions. The extension to the H =B case is straightforward with the help of heavy quark flavour symmetry (HQFS). Based on HQFS, the low energy constants (LECs) involved in the Dφ interactions or D → φ semileptonic form factors are related to their analogues in the bottom sector by specific scaling rules. It is then feasible either to predict quantities in the bottom (charm) sector by making use of the known information in the charm (bottom) case or to check how well HQFS works by testing the scaling rules.
In the present study, we construct the MO representations of the scalar form factors, denoted by f 0 (s), for the semileptonic D → π and D →K transitions, which are related to the unitarized S-wave scattering amplitudes in the Dφ channels with strangeness (S) and isospin (I) quantum numbers (S, I) = (0, 1 2 ) and (S, I) = (1, 0), respectively. These amplitudes are obtained by unitarizing the O(p 2 ) heavy-meson chiral perturbative ones [62], with LECs determined from the lattice calculation [56] of the S-wave scattering lengths in several (S, I) sectors. The scheme provides an accurate description of the Dφ interactions in coupled channels. For instance, as it is shown in Ref. [52], the finite volume energy levels in the (S, I) = (0, 1/2) channel calculated with the unitarized amplitudes, without adjusting any parameter, are in an excellent agreement with those recently reported by the Hadron Spectrum Collaboration [60]. In addition, it is demonstrated in Ref. [42] that these well constrained amplitudes for Goldstone bosons scattering off charm mesons are fully consistent with recent high quality data on the B − → D + π − π − final states provided by the LHCb experiment [63].
The unitarized chiral scattering amplitudes are used in this work as inputs to the dispersive integrals. However, these amplitudes are valid only in the low Goldstoneboson energy region. Hence, asymptotic behaviors at high energies for the phase shifts and inelasticities are imposed in the solution of the MO integral equations. The Omnès matrices obtained in this way incorporate the strong final state interactions, and the scalar form factors are calculated by multiplying the former by polynomials. The (a priori unknown) coefficients of the polynomials are expressed in terms of the LECs appearing at next-toleading order (NLO) in the chiral expansion of the form factors [64,65]. The scheme employed in the charm sector is readily extended to the bottom one. Afterwards, the LECs could be either determined by fitting to the results obtained in the LQCD analyses of the D → π(K) decays carried out in Refs. [16,17,21] or to the LQCD and LCSR com-binedB → π andB s → K scalar form factors reported in Refs. [15,[18][19][20] and [12][13][14], respectively. In both scenarios LQCD and LCSR results are well described using the MO dispersive representations of the scalar form factors constructed in this work.
However, our best results are obtained by a simultaneous fit to all available results, both in the charm and bottom sectors. As mentioned above, all of the LECs involved in theB (s) φ interactions orB (s) →φ semileptonic transitions are related to those in the charm sector by making use of the heavy quark scaling rules [65], which introduce some constraints between the polynomials that appear in the different channels. Thus, assuming a reasonable effect of the HQFS breaking terms, a combined fit is performed to the D → π/K and theB → π and B s → K scalar form factors, finding a fair description of all the fitted data, and providing reliable predictions of the different scalar form factors in the whole semileptonic decay phase space, which turn out to be compatible with other theoretical determinations by, e.g., perturbative QCD [66,67]. The results of the fit allow also to predict the scalar form factors for the D → η, D s → K and D s → η transitions in the charm sector, and for the very first time for theB → η decay. In some of these transitions, the form factors are difficult for LQCD due to the existence of disconnected diagrams of quark loops.
Based on the results of the combined fit, and taking advantage of the fact that scalar and vector form factors are equal at q 2 = 0, we extract all the heavy-light CKM elements and test the second-row unitarity by using the |V cb | value given in the PDG [1]. We also predict the different form factors above the q 2 −regions accessible in the semileptonic decays, up to energies in the vicinity of the involved thresholds, which should be correctly described within the employed unitarized chiral approach.
This work is organized as follows. In Subsec. II A, we introduce the definitions of the form factors for H 3 decays. A general overview of the MO representation of the scalar form factors is given in Subsec. II B, while the inputs for the MO problem and the solutions are discussed in Subsec. II C. In Subsec. II D, we derive the scalar as well as vector form factors at NLO in heavymeson chiral perturbation theory and then perform the aforementioned matching between the MO and the chiral representations close to the corresponding thresholds. Section III comprises our numerical results and discussions, with details of the fits given in Subsecs. III A and III B. With the results of the combined charm-bottom fit, in Sec. III C, we extract the related CKM elements and make predictions for the values of f + (0) in various transitions. Predictions for flavour-changing b → u and c → d, s scalar form factors above the q 2 −regions accessible in the semileptonic decays are given and discussed in Subsec. III D. We summarize the results of this work in Sec. IV. Finally, the heavy-quark scaling rules of the LECs involved in the Hφ interactions are discussed in Appendix A, while some further results for b → u form factors, obtained with quadratic MO polynomials, are shown in Appendix B.

A. Form factors in H 3 decays
For a semileptonic decay of the type H(p) → φ(p ) (p )ν (p ν ), the Lorentz invariant Feynman amplitude is proportional to where G F is the Fermi constant, and V Qq is the CKM matrix element corresponding to the flavour changing Q → q transition. The terms in the first and second curly brackets stand for the weak and hadronic matrix elements, respectively. In the hadronic matrix element, the axialvector part vanishes due to the parity conservation, while the remaining vector part is parametrized in a conventional form as where f + (q 2 ) and f 0 (q 2 ) are the vector and scalar form factors, respectively, with q µ = p µ −p µ and Σ µ = p µ +p µ .
Note that both form factors should be equal at q 2 = 0. As discussed in Ref. [68], they specify the P -wave (J P = 1 − ) and S-wave (J P = 0 + ) of the crossed-channel matrix elements, Both the scalar and vector form factors contribute to the differential decay rate, see e.g., Ref. [22]. Nevertheless, when the lepton mass is neglected, the differential decay rate in the H−meson rest frame can be simply expressed in terms of the vector form factor via It is then possible to extract the CKM element |V Qq | even for a single value of the four-momentum transfer, provided one simultaneously knows the vector form factor and the experimental differential decay width. A possible choice is q 2 = 0, where the scalar and vector form factors coincide, f + (0) = f 0 (0). In the next two subsections, Subsecs. II B-II C, we give specific details on the MO representation of the form factors. For brevity, in some occasions we will focus on the formalism for the case of the charm sector (H = D). The extension to the bottom sector (H =B) is straightforward using HQFS, though some aspects are explicitly discussed in Subsec. II C 2 b.

B. Muskhelishvili-Omnès representation
We now discuss the dispersive representation of the scalar form factors within the MO formalism. Throughout this work, isospin breaking terms are not considered, and therefore it is convenient to work with the isospin basis. Before proceeding, we first discuss the relation of the form factors expressed in the particle and isospin bases. We start defining the phase convention for isospin states: while the other states are defined with a positive sign in front of the |I, I 3 state. The form factors involving the c → d transition are those appearing in the D 0 → π − , D + → π 0 , D + → η, and D + s → K 0 semileptonic decays. (Note that D 0 → π − and D + → π 0 transitions are related by an isospin rotation.) The details of these form factors close to the zero recoil point, where the outgoing Goldstone boson is at rest, are greatly influenced by the πD, Dη, and D sK scattering amplitudes in the (S, I) = (0, 1/2) sector. Note that to be consistent with the convention of Refs. [56,62], we use πD instead of Dπ to construct the isospin 1/2 state. We then define the vector-column F (0,1/2) as (1, 0) sector. We thus define 2 : for which we will also use the notation Here, and again to be consistent with the convention of Refs. [56,62], we use DK instead of KD to construct the isoscalar state, which is just the opposite convention to that used for the isospin πD state. With these definitions, the unitarity relation for any of the F(s) can be compactly written as: where T (s) stands for the coupled-channel S-wave scattering amplitude in the corresponding (S, I) sector, which will be discussed further in Sec. II C. The diagonal matrix Σ(s) contains the phase space factors. For (0, 1/2), one has Σ(s) = diag σ Dπ (s), σ Dη (s), σ DsK (s) , whereas in the (1, 0) sector Σ(s) = diag (σ DK (s), σ Dsη ). The function σ ab (s) is defined as: with λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2yz − 2xz, the Kählén function. Invariance under time-reversal together with the optical-theorem leads to In this convention, the T -and S-matrices are related by The unitarity relation of Eq. (8) can be used to obtain dispersive representations for the form factors. We start considering the D → π transition in the single channel case (elastic unitarity) where the form factor satisfies with t πD (s) the πD S-wave elastic scattering amplitude. Equation (12) admits an algebraic solution [69], where P (s) is an undetermined polynomial, and the Omnès function Ω(s) is given by with δ(s) the elastic 0 + πD phase shift, in accordance with the Watson final state interaction theorem [70]. This was the scheme adopted in the previous studies carried out in Refs. [30][31][32][33][34][35][36]38]. For coupled channels the solution for F(s) takes the form: being P(s) a vector of polynomials with real coefficients and Ω(s) the Omnès matrix that satisfies 3 Im Ω(s + i ) = T * (s + i )Σ(s)Ω(s + i ) (18) which leads to the following unsubtracted dispersion relation: with s th the lowest threshold, which is referred to as the MO integral equation [72]. Taking a polynomial P (s) of rank (n − 1) would require the knowledge of F(s) for n values of s. Unlike the single channel case, there is no analytical solution [in the sense of Eq. (14)] for the coupledchannel MO problem and it has to be solved numerically.
The MO equation can be written in an alternative form: Re Im Ω(s) = X(s) Re Ω(s) , with P denoting the principal value and which is expressed in terms of the T -matrix and encodes the information on the Dφ re-scattering. The linear MO system, Eq. (20), can be solved by following the appropriate numerical method described in Ref. [71]. 3 Taking into account that the Ω(s) matrix should have only a right-hand cut and it should be real below all thresholds, Eq. (18) is equivalent to with H(s) = (I + 2iT (s)Σ(s)). Furthermore since H(s)H * (s) = I, though H(s) is not the S-matrix in the coupled-channel case, it follows with exp 2iφ(s) = det [H(s)]. This is to say that the determinant of the matrix Ω(s) satisfies a single-channel Omnès-type relation [71], which is extensively used in this work to check the accuracy of the numerical calculations. Note that above all thresholds, det [H(s)] = det [S(s)] and therefore in the elastic case (η i → 1 ∀i), φ(s) = n i=1 , δ i (s), with n the number of channels.

C. Inputs and MO solutions
To solve the MO equation and obtain the Ω(s) matrix, the T -matrix is needed as an input. We will use here the amplitudes based on unitarized chiral effective theory as computed in Refs. [56,62]. Because of the normalizations used in Eqs. (8) and (10), the unitarized amplitude in Ref. [56], denoted here by T U (s), is related to the Tmatrix introduced in the above subsection by: This unitary T -matrix is written as: where the elements of the diagonal matrix G(s) are the two-meson loop functions [56], and the matrix V (s) contains the interaction potentials which, as pointed out above, are taken from the O(p 2 ) chiral Lagrangians deduced in Ref. [62] (see Ref. [52] for more details). 4 In Ref. [56] all the parameters involved in the T -matrix, which are a few LECs and one subtraction constant, have been pinned down by fitting to the lattice calculation of the S-wave scattering lengths in several (S, I) sectors. The obtained interaction potentials successfully describe [52] the charm (S, I) = (0, 1/2) finite volume energy levels, without adjusting any parameter, calculated by the Hadron Spectrum Collaboration [60]. Moreover, as mentioned in the Introduction, these well constrained amplitudes turned out to be fully consistent [42] with the LHCb data on the B − → D + π − π − reaction [63].
Since the amplitudes are based on chiral potentials, the obtained T -matrices are only valid in the energy region not far from the corresponding thresholds. We thus adopt such T -matrices only up to a certain value of s, denoted by s m . Above that energy, the T -matrix elements are computed as an interpolation between their values at s = s m and the asymptotic values at s = ∞. The interpolation still gives a unitary T -matrix since, as we will specify below, it is actually performed on the phase shifts and the inelasticities. Moreover, the approximation of quasi-two-body channels cannot hold for arbitrarily large energies and Eq. (8) is a reasonable approximation to the exact discontinuity only in a finite energy range. However, as we are interested in constructing the form factors in a finite energy region also, the detailed behaviour of the spectral function at much higher energies should be, in principle, unimportant. As we will see below, this is not entirely correct, in particular for theB → π semileptonic transition, because of the large q 2 −phase space accessible in this decay. Nevertheless, we will assume that Eq. (8) holds up to infinite energies, only requiring that the T -matrix behaves in a way that ensures an appropriate asymptotic behaviour of the form factors, and we will discuss the dependence of our results on the contributions from the high energy region. In general, except for theB decays 5 , the asymptotic conditions on the T -matrix are chosen such that: where n is the number of channels involved in the Tmatrix and δ i (s) are the phase shifts. These conditions ensure (in general) that the unsubtracted dispersion relation for the Omnès matrix in Eq. (19) has a unique solution, albeit a global normalization [72] (see also in particular Sec. 4.3 of Ref. [71]). The condition of Eq. (25) guaranties that as can be deduced from the discussion of Eqs. (16) and (17). Note that the normalization of the Omnès matrix is completely arbitrary, and the computed form factors do not depend on it 6 .
In what follows, we detail the T -matrices and the specific shape of the asymptotic conditions for the two coupled-channel cases [(S, I) = (1, 0) and (0, 1/2)] that will be analyzed below. In this sector, we will consider two coupled channels, DK (1) and D s η (2), and above all thresholds, the Tmatrix is parametrized in terms of two phase shifts and one inelasticity parameter, with the phase φ 12 = δ 1 + δ 2 + mod(π) and 0 ≤ η ≤ 1.
To solve the MO integral equation [cf. Eq. (20)], we use 5 We will specifically discuss the situation for these transitions below. a T -matrix of the form with s th = (m D + M K ) 2 the lowest threshold, T U defined in Eq. (23) and T H the asymptotic matrix that will be discussed below. Phase shifts, inelasticities and amplitude moduli from T Fig. 1 up to √ s = 2.8 GeV, slightly above s m = (2.7 GeV) 2 . Above this scale, the T -matrix elements are computed as an interpolation between their values at s = s m and the asymptotic values at s = ∞, given in Eqs. (24) and (25). Thus, T H is constructed from Eq. (27) using the following parameterizations for phase shifts and inelasticities: as suggested in Ref. [71]. As discussed above, the Omnès matrix is uniquely determined by choosing η(∞) = 1 and δ 1 (∞) + δ 2 (∞) = 2π. The only remaining freedom is the distribution of 2π over the two phase shifts. Note that, δ i (s m ) is defined modulo π and this ambiguity is fixed by continuity-criteria. Here for the DK − D s η coupled channels, we choose δ 1 (∞) = 2π, δ 2 (∞) = 0. Different choices of the asymptotic values or of the interpolating functions in Eq. (29) will modify the shape of the Omnès solution far from the chiral region. The numerical effect of such freedom on the derived scalar form factors should be safely compensated by the undetermined polynomial in front of the Omnès matrix.
In Fig. 2, we show the solution of the MO integral equation (20), with the input specified above, and the contour condition Ω(q 2 max ) = I, with q 2 max = (m D − M K ) 2 . We display results only up to s = s m that would be later used to evaluate the scalar form factors entering in the D →K and D s → η semileptonic transitions. Note that the imaginary parts are zero below the lowest threshold s th = (m D + M K ) 2 , and how the opening of the D s η threshold produces clearly visible effects in the Omnès matrix. At very high energies, not shown in the figure, both real and imaginary parts of all matrix elements go to zero, as expected from Eq. (26). parameters, i.e., three phase shifts and three inelasticities [74,75], Furthermore, the parameters in the off-diagonal elements are related to the diagonal ones δ i and η i by and α ij is determined as Note that the solutions for α ij can be either arcsin(X ij ) or π − arcsin(X ij ). The inelasticity parameters should satisfy the following boundary conditions: With the inputs specified above, the three-dimensional (S, I) = (0, 1/2) Omnès matrix can be numerically computed and its complex elements are shown in Fig. 4 up to √ s ≤ 2.6 GeV.
b. Bottom sector: In Figs. 5 and 6, we show phase shifts, inelasticities and the solution of the MO integral equation for the (S, I) = (0, 1/2) channel in the bottom sector. The chiral amplitudes 8 are used in Eq. (19) up to s m = (6.25 GeV) 2 , and from there on, the asymptotic forms of the amplitudes are employed. As we will show in the next section, in the case ofB decays, the accessible phase space is quite large, and q 2 varies from around m 2  at zero recoil [q 2 max = (m B − M φ ) 2 ] down to zero, when the energy of the outgoing light meson is about m B /2 far from the chiral domain. TheBπ scalar form factor decreases by a factor of five, and the LQCD results around q 2 max and the LCSR predictions in the vicinity of q 2 = 0 are not linearly connected. In the present approach, as we will discuss, we multiply the MO matrix Ω by a rankone polynomial, and thus the extra curvature provided by the MO matrix becomes essential. While Ω(s) around q 2 max is rather insensitive to the adopted asymptotic behaviour of the T -matrix, since it is dominated by the integration region close to threshold (s < s m ) where the chiral amplitudes are being used 9 , this is not the case for low values of q 2 close to 0, quite far from the twobody scattering thresholds. This unwanted dependence, due to the large extrapolation, could be compensated in the form factors by using higher rank polynomials, but that would introduce additional undetermined parameters. Conversely, this dependence of Ω(0), relative to the results at q 2 max , on the details of the amplitudes at high energies will be diminished by solving a MO integral equation involving several subtractions, instead of the unsubtracted one of Eq. (19). This, however, will also introduce some more free parameters [34]. The situation is better in the charm sector, where the needed q 2 −range is much reduced, and thus most of the contributions to the MO matrix come from integration region within the chiral regime. Indeed in the bottom sector we need to use δ 1 (∞) = 2π, δ 2 (∞) = 2π and δ 3 (∞) = 0, instead of the choice of Eq. (33) used for the charm decays, to find acceptable fits to the LQCD and LCSR predictions of theBπ andB s K scalar form factors. With this choice, we find theoretically sound fits where the LECs 10 that determine the rank-one Omnès polynomials describe the LQCD data close to q 2 max , within the range of expected validity of the chiral expansion, while the LCSR results are reproduced thanks to the non-linear behaviour encoded in the MO matrix Ω(s). This picture will be reinforced by the consistent results that will be obtained, assuming a reasonable effect of the HQFS breaking terms, from combined fits to the D → π/K and theB → π and 10 These are β P 1 and β P 2 , to be introduced in Subsec. II D, that appear in the chiral expansion of the form factors at NLO.   We do not really have an explanation of why the above choice of phase shifts at infinity works better in the bottom sector than the usual one in Eq. (33) and adopted in the charm meson decays. We would like, however, to mention the different behaviour of the unitarized chiral phase shifts in the charm and bottom sectors. In both cases, the chiral amplitudes give rise to two resonances [52]: the first one, the non-strangeness flavor partner of the D * s0 (2317), quite broad, and located around 100 MeV above the Dπ orBπ thresholds and the second one placed below the heaviest of the thresholds, D sK andB sK , respectively. In the charm sector, the second resonance does not produce clear signatures in the phase shifts of the two open channels Dπ and Dη, while it is clearly visible in the phase shifts of the bottomBπ andBη channels. Moreover, the second resonance is significantly narrower for the latter heavy-quark sector than for the former one (70 MeV versus 270 MeV).
Note also that now lim s→∞ Re{Ω} Im{Ω} expect the form factors to vanish in this limit [76]. To achieve such asymptotic behaviour one should employ order zero polynomials. However, since we are interested in the region 0 s s max , with s max in the vicinity of (m H − M φ ) 2 , we prefer to keep the linear behaviour of the polynomials, since this allows for a better matching of the coefficients α 0,1 with the LECs that appear in the NLO chiral calculation of the form factors.

Form factors in heavy meson chiral perturbation theory
The leading-order (LO) coupling of the charm (D and D s ) or bottom (B andB s ) mesons to the Nambuin the bottom sector and the Omnès matrix elements are expected to decrease slightly faster than 1/s. Goldstone bosons of the spontaneous breaking of the approximate chiral symmetry of QCD, through the charged- , is described by the following chiral effective Lagrangian [64,65,77] where P and P * are the pseudoscalar and vector heavylight mesons with content (Qū, Qd, Qs), respectively, which behave as SU (3) light flavor triplets. Herem denotes the degenerate mass of the P (s) and P * (s) mesons in the chiral and heavy-quark limits, and f P is the pseudoscalar heavy-light meson decay constant defined as  The chiral block is defined by with F 0 the pion decay constant in the chiral limit (we will take the physical value for the decay constant F 0 92 MeV). The relevant NLO chiral effective Lagrangian reads [65] We need the LO PP * φ interaction as well, which is given by [64,65,77] , where 6 a dimensionless and heavy quark mass independent constant. The topologies of relevant Feynman diagrams are shown in Fig. 7. The vector and scalar form factors, in the (strangeness, isospin) basis, at O(E φ ) (i.e., NLO) in the chiral expansion read [65] f where ∆ P φ = m 2 P − M 2 φ and Σ P φ = m 2 P + M 2 φ with P ∈ {Qs, Qd, Qū} and φ ∈ {π, K,K η}. The coefficients C (S,I) [P φ] are collected in Table I. Moreover, we have fixed m R to m D * (m B * ) and to m D * s (m B * s ) for the (S, I) = (0, 1/2) and (S, I) = (1, 0) charm (bottom) channels, respectively. In principle, at LO in the heavy quark expansion, m R should be set tom, however the use of the physical vector mass is quite relevant for the vector form factor, because of the propagator structure, though it has much less relevance for the scalar form factor that we study in this work. We should also note that all kinematical factors are always calculated using physical masses Re{Ω} Im{Ω} of the involved mesons. It is worthwhile to notice that s is of the order m 2 so that a small change of O(E φ ) in (∆ P φ + s) only leads to a higher order effect. Thus, (∆ P φ + s) should be regarded as basically a constant with s ∼ m 2 P , and the expression in Eq. (41) should be matched to a rank-1 MO polynomial as mentioned before (see below for details of the matching).
Finally, we should mention that the LECs β P 1 and β P 2 scale with the heavy quark mass as [65] β P 1 ∼ √ m P , neglecting logarithmic corrections. 7. Topologies of the relevant Feynman diagrams contributing to the hadronic matrix elements. The solid circle denotes the LO PP * φ interaction and the solid square represents the left-handed current.

Matching
At energies close to the thresholds, the scalar form factors in Eq. (15) should have the same structure as the ones obtained from chiral perturbation theory, given in Eq. (41). We match the two representations at a point s = s 0 located in the valid region of the chiral expansion. Namely, we take s 0 = q 2 max = (m P −M π ) 2 for the (0, 1/2) form factors, and s 0 = (m D − M K ) 2 for the charm (1, 0) case, since this is the point in which the momentum of the lightest meson is zero. Imposing that the dispersive form factors and their first derivative to be equal to the chiral ones at s = s 0 , the coefficients in the polynomials, Eq. (34), can be expressed as: where the stands for a derivative with respect to s. The vectors F χ contain the chiral form factors, with all the elements given in Eq. (41). Similar expressions are used for the (0, 1 2 ) channel in the bottom sector. In other words, the vectors F χ (s) contain the form factors defined in Eqs. (6) and (7), but computed according to the chiral expansion.
It is worth noting that the NLO LECs β P 1 and β P 2 determined from a fit to data using the MO scheme would have some residual dependence on the matching point. To minimize such dependence, we have chosen s 0 = q 2 max , where the momentum of the Goldstone bosons is close to zero and higher order chiral corrections are expected to be small. Different choices of the matching point, within the chiral regime, will amount to changes in the fitted (effective) β P 1 and β P 2 LECs driven by higher order effects.
In the charm (1, 0) sector, due to the presence of the D * s0 (2317) state as a bound state in the T -matrix, the solution of Eq. (15) gets modified. The contribution from D * s0 (2317) is easily incorporated as follows: where β 0 is an unknown parameter which characterizes the coupling of D * s0 (2317) to the left-hand current. Γ contains the couplings of the D * s0 (2317) to the DK-D s η sys- tem 12 , namely, Γ = (g DK , g Dsη ) T . This bound state is dynamically generated in the unitarized amplitudes given in Ref. [56], that we employ here. The couplings are computed from the residue of the amplitude at the pole, The D * s0 (2317) pole position s p , together with g DK and g Dsη , are collected in Table II.

III. NUMERICAL RESULTS AND DISCUSSION
So far, the theoretical MO representations of the scalar form factors have been constructed. In this section, we want to confront the so-obtained form factors to the LQCD and LCSR results. In what follows, we first fit to theB → π andB s →K scalar form factors, where we expect the 1/m P corrections to the chiral expansion in Eq. (41) to be substantially suppressed. Next, we will carry out a combined fit to all the data in both charm and bottom sectors, by adopting some (approximate) heavyquark flavor scaling rule [65] for the β P 1 and β P 2 LECs in Eq. (41). Using the results of our combined fit, we 12 Note that the first term in the bracket of Eq. are not exactly independent of the choice of the point sn where one normalizes the Omnès matrix, Ω(sn) = I. Nonetheless, we have checked that this choice, varying sn from zero to q 2 max , has no practical effect in the determination of β 0 , which indicates that our assumption is reasonable. We also remark that this discussion has no effect at all in the (0, 1/2) sector. will: i) determine the CKM elements, |V cd |, |V cs | and |V ub |, ii) predict form factors, not computed in LQCD yet, and that can be used to over-constrain the CKM matrix elements from analyses involving more semileptonic decays, and iii) predict the different form factors above the q 2 −regions accessible in the semileptonic decays, up to moderate energies amenable to be described using the unitarized coupled-channel chiral approach.
Masses and decay constants used in this work are compiled in Table III. In addition, the mass of the heavylight mesons in the chiral limit, see Eq. (35), is set to m = (m P + m Ps )/2, for simplicity the same average is used to definem P in the Appendix A and in the relations given in Eq. (53). The PP * φ axial coupling constantg in Eq. (39) can be fixed by calculating the decay width of D * + → D 0 π + [1], which leads 13 to g ∼ 0.58 and hencẽ g D * Dπ ∼ 1.113 GeV. In the bottom sector we use a different value for g, around 15% smaller, consistent with the lattice calculation of Ref. [78], where g ∼ 0.51 (or gB * B π ∼ 2.720 GeV) was found. Note that the difference is consistent with the expected size of heavy-quarkflavor symmetry violations. In addition, there exist sizable SU (3) corrections to the overall size of the P * pole contribution to both f + and f 0 form factors. Thus, such contribution is around ∼ 20% smaller forB * B s K than forB * B π [32,38]. According to [79] this suppression is mainly due to a factor F π /F K ∼ 0.83 [22]. We will implement this correction in the pole contribution to f 0 in Eq. (41) when the Goldstone boson is either a kaon or an eta meson (for simplicity, we also take F η ≈ F K ), and both in the bottom and charm sectors.

A. Fit to the LQCD+LCSR results in the bottom sector
We are first interested in theB → φ transitions induced by the b → u flavour-changing current, which in-cludeB → π,B → η andB s → K. The scalar form factors involved in those transitions can be related to the Omnès matrix through Within the present approach, and considering just rankone MO polynomials, there are only two undetermined parameters: the NLO LECs β P 1 and β P 2 that appear in the chiral expansion of the form factors in Eq. (41). We fit these parameters, in the bottom sector, to LQCD (UKQCD [20], HPQCD [15,18] and Fermilab Lattice & MILC (to be referred to as FL-MILC for brevity) [19]) 13 Errors on g determined from the decay D * + → D 0 π + are very small of the order of 1%. and LCSR [12,13] results for the scalar form factors in B 0 → π + andB 0 s → K + semileptonic decays. Lattice results are not available for the whole kinematic region accessible in the decays, and they are restricted to large values of q 2 ≥ 17 GeV 2 , where momentum-dependent discretization and statistical errors are under control. To constrain the behaviour of the scalar form factors at small values of q 2 , we take four LCSR points (equally-spaced) in the interval q 2 = 0 − 6 GeV 2 for each decay.
The UKQCD Collaboration [20] provides data for both B → π andB s → K form factors together with statistical and systematic correlation matrices for a set of three form factors computed at different q 2 (Tables VIII and IX of this reference). In the case of HPQCD [18]B s → K and FL-MILC [19] B → π form factors, we have read off four points from the final extrapolated results (bands) given in these references, since in both cases, originally only four momentum configurations (0 → 0, 0 → 1, 0 → √ 2 and 0 → √ 3) were simulated. Finally, we also include in the fit the five B → π points provided by the HPQCD Collaboration in the erratum of Ref. [15].
Thus, the χ 2 function reads where χ 2 is the usual uncorrelated Gaussian merit function, and χ 2 cov is defined as with the covariance matrix C constructed out the statistical and systematic correlation matrices and uncertainties given in Ref. [20]. Here f 0 (q 2 ) stands for the theoretical scalar form factor obtained from the MO representation. The chi-squared fit results are with χ 2 /dof = 4.2 for 25 degrees of freedom, and a correlation coefficient 0.999 between the two fitted parameters. The first set of errors in the parameters is obtained from the minimization procedure, assuming Gaussian statistics, while the second one accounts for the uncertainties of the LECs quoted in Ref. [56] that enter in the definition of the chiral amplitudes. Such a correlation coefficient so close to 1 indicates that the considered data can not properly disentangle both LECs, 14 and that different (β B 1 , β B 2 ) pairs belonging to the straight line in the vicinity of the best fit values (β B 1 , β B 2 ) quoted in Eq. (51) lead to similar descriptions of the data (see the dashed-blue line in the right panel of Fig. 10). The scalar form factors obtained are displayed in Fig. 8. We find a fair description of the LQCD and LCSR results for thē B 0 s → K + scalar form factor, while we face some problems for theB 0 → π + decay. The large value of χ 2 /dof reported in Eq. (51) is mainly due to the existing tension between the LQCD results from different collaborations in this latter decay. The disagreement between UKQCD and HPQCDB 0 → π + scalar form factors was already highlighted in the top-right panel of Fig. 23 of the UKQCD work [20], where it is noted that the HPQCD calculation used only a single lattice spacing.
In addition, as we discussed before, theBπ−scalar form factor decreases by a factor of five in the q 2 −range accessible in the decay, and the LQCD results around q 2 max and the LCSR predictions in the vicinity of q 2 = 0 are not linearly connected at all. In the current scheme, where only rank-one MO polynomials are being used, this extra needed curvature should be provided by the q 2 −dependence of the MO matrix, Ω, whose behaviour near q 2 = 0, far from q 2 max , is not determined by the behaviour of the amplitudes in the chiral regime. Indeed, it significantly depends on the high-energy input 15 . This is an unwanted feature, source of systematic uncertainties. To minimize this problem, in the next subsection we will perform a combined fit to transitions induced by the b → u and c → d, s flavour-changing currents. The latter ones describe D → π and D →K semileptonic decays for which there exist recent and accurate LQCD determinations of the scalar form factors. Moreover, in these latter transitions the q 2 −ranges accessible in the decays and the form factor variations are much limited, becoming thus more relevant the input provided in the chiral regime.
To finish this subsection, we would like to stress that given the large value found for χ 2 /dof , statistical errors should be taken with some care. Indeed, one can rather 15 The results displayed in Fig. 8 might suggest that the present approach hardly provides enough freedom to simultaneously accommodate the near q 2 = 0 (LCSR) and q 2 max (LQCD) determinations of theBπ scalar form factor. The situation greatly improves when only the HPQCD, among all LQCD calculations, Bπ results are considered in the q 2 max region, being then possible to find an excellent combined description of the LCSR and HPQCD results with χ 2 = 9.65 for a total of 18 degrees of freedom (see dashed-red curve in the right plot of Fig. 10), which leads to χ 2 /dof = 0.5. The parameters β B 1,2 come out still to be almost totally correlated as in Eq. (51), and moreover they lie, within great precision, in the straight line of Eq. (52), but in the β B 1 ∼ 0.7 GeV region. assume some systematic uncertainties affecting our results, that could be estimated by considering in the best fit alternatively only the HPQCD or the UKQCD and the FL-MILC sets of predictions. We will follow this strategy to obtain our final results for the CKM matrix elements and form factors at q 2 = 0 from the combined-fit to charm and bottom decays detailed in the next subsection.

B. Extension to the charm sector and combined fit
Besides the parameter β 0 introduced in Eq. (46) to account for the effects on the D * s0 (2317) state in the c → s decays, one should also take into account that the LECs β P 1 and β P 2 depend on the heavy quark mass. The scaling rules given in Eq. (42) can be used to relate the values taken for these LECs in the bottom (β B i ) and charm sectors (β D i ). We will assume some heavy quark flavor symmetry violations and we will use where, one should expect the new parameter, δ, to be of the order Λ QCD /m D . Note that we are correlating the heavy quark flavor symmetry violations in the LECs β 1 and β 2 . There is not a good reason for this other than avoiding to include new free parameters. On the other hand, at the charm scale, one might also expect sizable corrections to the LO prediction f 0 (s) ∼ C × f P /F 0 of Eq. (41), even more bearing in mind the large (40 − 50%) heavy-quark symmetry violations inferred from the ratio f B /f D quoted in Table III. (Note that at LO in the inverse of the heavy quark mass, this ratio should scale as (m D /m B ) 1/2 ). Thus, we have also introduced an additional parameter, δ , defined through the replacement when Eq. (41) is applied to the c → d, s decays. Thus, we have three new parameters β 0 , δ, δ , which in addition to β B 1,2 , will be fitted to the LQCD & LCSR results for the scalar form factors in theB → π,B s → K, D → π and D →K semileptonic decays.
First we need to incorporate the c → d, s input into the merit function χ 2 , which was defined in Eq. (49) using only bottom decay results. In the last ten years, LQCD computations of the relevant D → π and D →K semileptonic decay matrix elements have been carried out by the HPQCD [16,17] and very recently by the ETM [21] Collaborations. Compared with the former, the latter corrects for some hypercubic effects, coming from discretization of a quantum field theory on a lattice with hypercubic symmetry [80], and uses a large sample of kinematics, not restricted in particular to the parent D meson at rest, as in the case of the HPQCD simulation. Moreover, it is argued in Ref. [21] that the restricted kinematics employed in the simulations of Refs. [16,17] [20], HPQCD [15,18], FL-MILC [19] and LCSR [12,13]), and for comparison, predictions from the NLO perturbative QCD approach of Ref. [66] for theB 0 → π + decay are also shown. Statistical (stat) and statistical plus systematic (stat & sys) 68%-confident level (CL) bands are also shown. The systematic uncertainties are inherited from the errors on the LECs quoted in Ref. [56], that enter in the definition of the chiral amplitudes, and are added in quadratures to the statistical uncertainties to obtain the outer bands. To estimate the systematic uncertainties for each set of LECs we re-do the best fit.
obscure the presence of hypercubic effects in the lattice data, and these corrections can affect the extrapolation to the continuum limit in a way that depends on the specific lattice formulation. This might be one of the sources of the important discrepancies found between the D → π form factors reported by the HPQCD and ETM Collaborations in the region close to q 2 max = (m D − M π ) 2 , as can be seen in the left top panel of Fig. 9.
Here, we prefer to fit to the most recent data together with the covariance matrices provided by the ETM Collaboration. This analysis is based on gauge configurations produced with N f = 2+1+1 flavors of dynamical quarks at three different values of lattice spacing, and with pion masses as small as 210 MeV. Lorentz symmetry breaking due to hypercubic effects is clearly observed in the ETM data and included in the decomposition of the current matrix elements in terms of additional form factors. Those discretization errors have not been considered in the HPQCD analyses, and for this reason we have decided to exclude the results of these latter collaboration in our fits.
The scalar form factors involved in the D → π and D →K transitions are related to the Omnès matri-ces displayed in Figs. 2 and 4 through 16 Eq. (15) and Eqs. (6) and (7). Hence, the bottom-charm combined χ 2 now reads where we have added sixteen ETM points, eight for each of the two D → π and D →K decay modes. Each of the new eight-point sets is correlated and the corresponding covariance matrices 17 have been obtained from the authors of Ref. [21]. Thus, we are fitting five parameters to a total of 43 points. The best-fit results for the five unknown parameters and their Gaussian correlation matrix are collected in Table IV and the resulting scalar form factors are shown in Fig. 9.
The results for the bottom scalar form factors are almost the same as the ones shown in Fig. 8, while the ETM c → d, s transition form factors are remarkably well described within the present scheme. As in the former bestfit to only theB (s) results, the large value obtained for χ 2 /dof is mainly due to the existing tension between the LQCD results from different collaborations in theB → π decay.
Due to the hypercubic effects, there might be inconsistencies between the ETM and HPQCD analyses for the Dπ scalar form factor in the region close to q 2 max = (m D −M π ) 2 [21]. As one can see, our result disagrees with the HPQCD data in that region too. We have checked that if we fit to the HPQCD instead of the ETM data in the charm sector, the best fit still tends to coincide with the ETM data. This observation is important and it seems to indicate that the Lorentz symmetry breaking effects in a finite volume, due to the hypercubic artifacts, could be important in the LQCD determination of the form factors in semileptonic heavy-to-light decays, as pointed out in Ref. [21].
The HQFS breaking parameters δ and δ turn out to be quite correlated and their size is of the order Λ QCD /m c . As expected, δ presents also a high degree of correlation with β B 1 and β B 2 , and on the other hand, the combined fit does not reduce the large correlation between these two latter LECs, while the central values (errors) quoted for them in Table IV are compatible within errors with (significantly smaller than) those given in Eq. (51), and obtained from the fit only to b → u transitions. In addition, the values quoted for (β B 1 , β B 2 ) in Table IV perfectly lie in the straight line of Eq. (52), deduced from the fit to only bottom form factors carried out in the previous Subsect. III A. Indeed, the straight line that one can construct with the results of Table IV in the (β B 1 , β B 2 )−plane is practically indistinguishable from that of Eq. (52). All this can be seen in the left plot of Fig. 10, where both straight lines are depicted, together with the statistical 68% CL ellipses and the one-sigma-rectangle bands obtained by minimizing the merit function given in Eq. (49) or alternatively in Eq. (55), and considering only bottom or bottom and charm scalar form factors, respectively.
In the right plot of Fig. 10, we show the dependence of χ 2 on β B 1 for different situations. We display the combined charm-bottom and the bottom-only fits, and in both cases, we have considered results obtained when all B → π LQCD form factors or only the HPQCD or the UKQCD and the FL-MILC subsets of results are considsince a tiny error in the fitting function yields a huge χ 2 value (specific examples can be found in Ref. [81]). We have used the singular value decomposition (SVD) method to tackle this issue, which is widely used by a number of lattice groups [82][83][84]. ered in the fits. The circles stand for the different best-fit results, accounting for variations of χ 2 up to one unit from the minimum value 18 , while the dashed and solid curves have been obtained by relating β B 1 and β B 2 through Eq. (52) and minimizing χ 2 with respect to the other parameters, β 0 , δ and δ , respectively. Several conclusions can be extracted from the results shown in the figure: • The combined charm-bottom analyses (solid lines) provide large curvatures of χ 2 as a function of β B 1 , hence leading to better determinations of this latter LEC, always in the 0.2 GeV region, as we also found in Eq. (49) from the best fit to only the bottom results. A value for β B 1 close to this region, taking into account errors, is also found from a fit where only the bottom form factors are considered, but without including the HPQCDB → π results (dashedgreen curve). Only the dashed-red line (fit only to the bottom results, but without including in this case the UKQCD and FL-MILCB → π scalar form factors) turns out to be incompatible with the combined fit presented in Table IV. Thus, we find some arguments to support the range of values quoted in Table IV for the parameters (β B 1 , β B 2 ) that appear in the heavy meson chiral perturbation theory (HMChPT) expansion of the scalar form factors at the bottom scale.
• The existing tension between HPQCD, and UKQCD and FL-MILC sets ofB → π form factors leads to large values of χ 2 . Thus, as mentioned above, statistical errors should be taken with some care, and some systematic uncertainties would need to be considered in derived quantities, as for instance in the values of the form factors at q 2 = 0 or in the CKM mixing parameters. We note that this source of systematics also induces variations on the fitted parameters in Table IV which range between 50% (β 0 and β B 2 ) to 100% (β B 1 , δ and δ ) of the statistical errors quoted in the table.
Our predictions for the scalar form factors for the D s → η, D → η and D s → K transitions, for which there are no lattice results as yet, are also shown in Fig. 9. Note that transitions involving the η meson in the final state are more difficult to be evaluated in LQCD simulations. Interestingly, the D → η scalar form factor in the three-channel (0, 1/2)−case is largely suppressed, similar to the component regarding the K → η transition in the strangeness-changing scalar form factors as shown in Ref. [26].  Table IV for details). The three bottom panels are similar to those depicted in Fig. 8, but computed from the results of the combined best-fit. The four panels in the first two rows show form factors for c → d, s semileptonic transitions. Only ETM results, corrected for some hypercubic (discretization) effects [21], have been considered in the fit of Table IV. For comparison, predictions from the HPQCD [16,17] Collaboration are also displayed. Differences between ETM and HPQCD sets of D → π and D → K form factors are clearly visible in the vicinity of q 2 max , in particular for the D → π case. Statistical (stat) and statistical plus systematic (stat & sys) 68%-confident level bands are also given and are calculated as explained in Fig. 8. Finally, predictions for the D → η, Ds → K and Ds → η scalar form factors are also shown.

Further considerations
We have also obtained results using constant and quadratic MO polynomials. In the first case, the dispersive representations of the form factors should be matched to the LO chiral ones, where the terms driven by the β P 1 and β P 2 LECs are dropped out. The first consequence is that bottom and charm sectors are no longer connected since, in addition, we are not enforcing the heavy-quark scaling law for the decay constants. To better describe the data, one might perform separate fits to bottom and charm form factors with free α 0 parameters in Eq. (34). Fits obviously are poorer, and moreover, they do not necessarily provide reliable estimates of the form factors at q 2 max , since the fitted parameters are obtained after minimizing a merit function constructed out of data in the whole q 2 −range accessible in the decays. For charm decays, the description of the D → π form factor is acceptable, while that of the D →K is in comparison worse, mostly because the s−dependence induced by the D * s0 (2317) can not be now modulated by the MO polynomial. In the bottom sector, as one should expect, the simultaneous description of LQCD and LCSR form factors in the vicinity of q 2 max and q 2 = 0, respectively, becomes poorer. Indeed, since the LQCD input has a larger weight in the χ 2 than the LCSR one, the latter form factors are totally missed by the new predictions, which now lie below the lower error bands of the LCSR results.
The consideration of quadratic MO polynomials solves this problem, as shown in Fig. 12 of Appendix B. Indeed, it is now possible to improve the description of theB → π LCSR form factors, providing still similar results in the q 2 max −region, where the LQCD data are available. Thus, for instance, we get fB →π + (0) = 0.248(10) using the new fit to be compared with 0.211(10) obtained using the parameters of the fit of Eq. (51) (form factors displayed in Fig. 8). Nevertheless, as we will see in the next subsection, there exist some other systematic errors, which practically account for the latter difference, and thus this source of uncertainty will be considered in the determination of the CKM matrix element |V ub |. In addition, though the χ 2 /dof obtained with quadratic MO polynomials is better, it is still large (around 3.7) due to the tension between theB → π LQCD results from different collaborations. Moreover β B 1 and β B 2 are still fully correlated, and the quadratic terms of the MO polynomials that multiply the elements Ω ij , (i = 1, 2, 3, j = 2, 3) of the matrix displayed in Fig. 6 are almost undetermined (see the large errors in the parameters α 2,3 and especially α 2,2 given in Table VI). Finally, the centralvalue predictions, that will be shown in Subsec. III D, for the form factors above q 2 max and to moderate energies amenable to be described using the unitarized coupledchannel chiral approach, are not affected by the inclusion of quadratic terms in the MO polynomials, though errors are enhanced. For all of this, we consider our best estimates for the form factors those obtained using rank-one polynomials.
We do not discuss quadratic terms in the charm-sector because rank-one MO polynomials led already to excellent reproductions of the form factors (see Fig. 9), in part due to the smaller q 2 −range involved in these decays. Moreover a correct charm-bottom combined treatment will require the matching at next-to-next-to-leading order (NNLO) in the chiral expansion, which is beyond the scope of this work.

C. Extraction of CKM elements and predictions
Taking advantage that scalar and vector form factors are equal at q 2 = 0, the results of the combined charmbottom fit presented in the previous subsection can be used to extract the vector form factor, f + , at q 2 = 0 for various semileptonic decays studied in this work. Moreover, given some experimental input for the quantity |V Qq |f + (0), with Qq = bu, cd or cs, we can extract the corresponding CKM matrix element using the present MO scheme. Measurements of the differential distribution dΓ(H →φ ν )/dq 2 at q 2 = 0 will directly provide model independent determinations of |V Qq |f + (0), 19 while measurements of the total decay width could be used to estimate this latter quantity only after relying on some model for the q 2 −dependence of f + .
In the charm sector from the fit presented in Table IV and Fig. 9, we find f D→π f D→K where the first and second sets of errors are similar to those quoted in Table IV and account for statistical (propagated from the 1σ fluctuations of the fitted parameters) and chiral systematical (propagated from the errors of the LECs that enter in the computation of the MO matrix) uncertainties, respectively. The third set of errors (sys 2 ) takes into account the variations that are produced when in the best fit one considers alternatively only HPQCD or UKQCD and FL-MILCB → π form factors. The results of Eq. (56) and (57) are in good 19 Neglecting the lepton masses. agreement with our preliminary estimates reported in [85], where we fitted only to the charm ETM LQCD form factors.
In combination with the experimental values f D→π taken from the report by the Heavy Flavor Averaging Group (HFLAV) [86], we obtain for the corresponding CKM matrix elements. The dominant error is the theoretical one affecting the determination of the form factors at q 2 = 0, within the scheme presented in this work. As expected, these results nicely agree with those reported by the ETM Collaboration [21] since we describe rather well the LQCD scalar form factors calculated in this latter work. The values of Eq. (59) agree within around 1σ with the average ones 20 given in Ref. [1], |V cs | = 0.995 (16) .
On the other hand, the test of the second-row unitarity of the CKM matrix is satisfied within errors where |V cb | = 0.0405(15) from PDG [1] has been used. Likewise, in the bottom sector, we obtain from the combined bottom-charm fit fB →π where we see that the error budget is now dominated by the inconsistency between HPQCD and UKQCD and FL-MILC sets of results for fB →π 0 for high q 2 , above 17 GeV 2 . Dropping out the UKQCD and FL-MILC sets of results forB → π, the LCSR and HPQCD results for this transition can be significantly better described simultaneously, leading to values of the form factor at q 2 = 0 around 0.24 for the combined charm-bottom fit, in the highest edge of the interval quoted in Eq. (64), and compatible within errors with the result of 0.26 +0.04 −0.03 predicted in Ref. [12] using LCSR 21 . However, the description of theB s → K and D → π scalar form factors gets somewhat worse, being thus the situation unclear. 20 Determinations from leptonic and semileptonic decays, as well as from neutrino scattering data in the case of |V cd |, are used to obtain the PDG averages. 21 Fitting only to the b → u data and not considering UKQCD and FL-MILC sets ofB → π results, we find fB →π 0 (0) ∼ 0.27, even in better agreement with the LCSR determination.
In principle, based on the above values, the CKM element |V ub | could be determined as in the charm sector. However, the full kinematic region in the bottom case is very broad, and the experimental determination of f + (0)|V ub | might suffer from large systematical uncertainties. A customary way to extract |V ub | has been to perform a joint fit to the LQCD and LCSR theoretical results for f + (q 2 ) and to measurements of the differential decay width, with |V ub | being a free parameter, see, e.g., Refs. [33,87,88]. This is not feasible to us, since we only know the value of the vector form factor at zero momentum by using the relation f + (0) = f 0 (0). However the latest Belle [9] and BaBar [10] works reported accurate measurements of theB → π partial branching fractions in several bins of q 2 that are used to extract the f + form factor shape and the overall normalization determined by |V ub |. As a result, Belle and BaBar obtained values of (9.2 ± 0.3) × 10 −4 and (8.7 ± 0.3) × 10 −4 for f + (0)|V ub |, respectively. Though the latter values were extracted from direct fits to data, they might be subject to some systematic uncertainties, since they were obtained using some specific q 2 parameterizations (Becirevic and Kaidalov [89] and Boyd-Grinstein-Lebed [90] in the Belle and BaBar works, respectively). Nevertheless, we average both determinations and we take fB →π Using this latter value and our estimate for the form factor at q 2 = 0 given in Eq. (64), we get There exist tensions between the inclusive and exclusive determinations of |V ub | [1]: and combining both values, R. Kowalewski and T. Mannel quote an average value of in the PDG review [1], which is in good agreement with our central |V ub | result of Eq. (66). We should mention that it is higher than the typical values obtained from LQCD and LCSR determinations of theB → π f + (q 2 ) form factor, combined with measurements of the q 2 distribution of the differential width. Thus, the FLAG review [22] gives an average value (in 10 3 units) of 3.67 ± 0.09 ± 0.12. Nevertheless, this latter value is still compatible, taking into account the uncertainties, with our result. These extractions of the CKM elements rely strongly on the results either from LQCD in the high q 2 region or from LCSR in the vicinity of q 2 = 0 (the latter only in the bottom sector), which are used in the combined fit, and hence are not ab initio predictions. However, our extractions incorporate the influence of general S-matrix properties, in the sense that unitarity and analyticity are implemented in the MO representation of the scalar form factors. Moreover, one of the advantages of our approach is that we can make predictions for the channels related by chiral SU (3) symmetry of light quarks. In some of these channels, the form factors are difficult for LQCD due to the existence of disconnected diagrams of quark loops. The D → η, D s → K, D s → η andB → η scalar form factors were already shown in Fig. 9 for the whole kinematical regions accessible in the decays. On the other hand, their values at q 2 = 0 are particularly important, since they might serve as alternatives to determine the CKM elements when experimental measurements of the corresponding differential decay rates become available. Our predictions for the absolute values of the vector form factors at q 2 = 0 are (we remind once more here that vector and scalar form factors coincide at q 2 = 0) |f D→η |f Ds→K |f Ds→η |fB →η |fB s→K For the decayB s → K, we find, adding errors in quadrature, fB s→K + (0) = 0.30 ± 0.03 in perfect agreement with the results obtained from the LCSR (0.30 +0.04 −0.03 [13]) and HPQCD (0.32 ± 0.06 [18]) analyses, but about 1 sigma above the LQCD result of the UKQCD Collaboration [20]. The single-channel Omnès-improved constituent quark model study of Ref. [38] led to 0.297 ± 0.027, which is also in good agreement with our result. D. Scalar form factors above the q 2 max −region It is worth recalling here the relation between the results obtained for the form factors and the scattering amplitudes used as input of the MO representation. If we focus, for instance, on the charm form factors, the lightest open-charm scalar resonance, called D * 0 (2400) by the PDG [1], lies in the (S, I) = (0, 1/2) sector. In Refs. [43,45,47], two different states, instead of only one, were claimed to exist in the energy region around the nominal mass of the D * 0 (2400). These studies were based on chiral symmetry and unitarity. This complex structure should be reflected in the scattering regime of the form factors. Indeed, this can be seen in the first row of panels of Fig. 11, where form factors for different semileptonic transitions are shown above the q 2 max −region. As discussed in Sec. II C, here we use the O(p 2 ) HMChPT amplitudes obtained in Refs. [56,62], which also successfully describe the (0, 1/2) finite-volume energy levels reported in the recent LQCD simulation of Ref. [60] (see Ref. [52] for details) and are consistent with the precise LHCb data [63] for the angular moments of the B − → D + π − π − [42]. These chiral amplitudes predict the existence of two scalar broad resonances, instead of only one, with masses around 2.1 and 2.45 GeV, respectively [42,52], which produce some signatures in the D → π, D → η and D s → K form factors at around q 2 = 4.4 and 6 GeV 2 , as can be appreciated in Fig. 11. The effect of this two-state structure is particularly visible in the D s → K form factor. Note that this two-state structure should have also some influence in the region below q 2 max , where we have fitted the LQCD data. Below q 2 max , the sensitivity of the form factors to the details of the two resonances is however smaller than that of the energy levels calculated in the scattering region, since the former ones are given below the lowest threshold, while the latter ones are available at energies around and above it. Nonetheless, the success in describing the LQCD results for the D → π scalar form factor clearly supports the chiral input, and the predictions deduced from it, used in the current scheme. If better determined form factors were available in all of the channels, perhaps the two state structure for the D * 0 (2400) could be further and more accurately tested.
A similar pattern is found in the bottom sector [42,52], as expected from the approximate heavy-flavor symmetry of QCD. The two-state structure is clearly visible, more than that in the charm sector, in the corresponding form factors (three bottom plots of Fig. 11), and it has a certain impact in the form factors close to q 2 max , where LQCD results are available.
In the charm (S, I) = (1, 0) sector the effect of the narrow D * s0 (2317) resonance, which is the SU (3) flavor partner of the lighter one of the two D * 0 states, predicted by the unitarized NLO chiral amplitudes [42,52], is clearly visible in the scalar D →K and D s → η form factors, and it fully dominates these form factors in the vicinity of the pole, as can be seen in the second row of panels of Fig. 11. Indeed, this state also influences the D →K form factor below (near) q 2 max , where the LQCD results are available 22 , and the excellent description of the ETM results gives clear support to the coupled-channel MO representation of the D →K and D s → η scalar form factors derived in this work.

IV. SUMMARY AND OUTLOOK
We have studied the scalar form factors that appear in semileptonic heavy meson decays induced by the flavourchanging b → u and c → d, s transitions using the 22 Indeed, the existence of the D * s0 (2317) was suggested in [34] by fitting the single channel MO representation of the D →K scalar form factor, constructed out the unitarized LO chiral elastic DK amplitude, to LQCD results of the scalar form factor below q 2 max .
MO formalism. The coupled-channel effects, due to rescattering of the Hφ (H = D,B) system, with definite strangeness and isospin, are taken into account by solving coupled integral MO equations. We constrain the subtraction constants in the MO polynomials, which encodes the zeros of the form factors, thanks to light-quark chiral SU (3) and heavy-flavor symmetries. The Hφ interactions used as input of the MO equations are well determined in the chiral regime and are taken from previous work. In addition, some reasonably behaviors of the amplitudes at high energies are imposed, while appropriate heavy-flavor scaling rules are used to relate bottom and charm form factors. We fit our MO representation of the scalar form factors to the latest c → d, s and b → u LQCD and b → u LCSR results and determine all the involved parameters, in particular the two LECs (β P 1 and β P 2 ) that appear at NLO in the chiral expansion of the scalar and vector form factors near q 2 max , which are determined in this work for first time. We describe the LQCD and LCSR results rather well, and in combination with experimental results and using that f 0 (0) = f + (0), we have also extracted the |V ub |, |V cd | and |V cs | CKM elements, which turn out to be in good agreement with previous determinations from exclusive decays.
We would like to stress that we describe extremely well the recent ETM D → π scalar form factor, which largely deviates from the previous determination by the HPQCD Collaboration, providing an indication that the Lorentz symmetry breaking effects in a finite volume, due to the hypercubic artifacts, could be important in the LQCD determination of the form factors in semileptonic heavy-tolight decays, as claimed in Ref. [21]. As it is also pointed out in the previous reference, this is a very important issue, which requires further investigations, since it might become particularly relevant in the case of the determination of the form factors governing semileptonicB−meson decays into lighter mesons.
We have also predicted the scalar form factors, which are in the same strangeness-isospin multiplets as the fitted D → π, D →K,B → π andB s → K ones. Our prediction of the form factors in such channels (D → η, D s → K, D s → η, and B → η) are difficult for LQCD simulations due to the existence of disconnected diagrams. These form factors are related to the differential decay rates of different semileptonic heavy meson decays and hence provide alternatives to determine the CKM elements with the help of future experimental measurements.
Moreover, we also find that the D → η scalar form factor is largely suppressed compared to the other two components (D → π, D →K) in the three-channel (0, 1/2)−multiplet, which is similar to what occurs for the K → η strangeness-changing scalar form factor in Ref. [26].
Our determination of the form factors has the advantage that the constraints from unitarity and analyticity of the S-matrix have been taken into account, as well as the state-of-the-art Hφ chiral amplitudes. Thus, our pre-  [56,62], and the LECs, compiled in Table  IV, obtained from a fit to LQCD and LCSR results below q 2 max . Statistical (stat) and statistical plus systematic (stat & sys) 68%-confident level bands are also given and are calculated as explained in Fig. 8.
dictions for the flavour-changing b → u and c → d, s scalar form factors above the q 2 −region accessible in the semileptonic decays, depicted in Fig. 11, should be quite accurate 23 and constitute one of the most important findings of the current research. Indeed, we have shown how the form factors in this region reflect details of the chiral dynamics that govern the Hφ amplitudes, and that give rise to a new paradigm for heavy-light meson spectroscopy [42] which questions the traditional qq constituent quark model interpretation, at least in the scalar sector.
As an outlook, the scheme presented here will also be useful to explore the Hφ interactions by using the lattice data for the scalar form factors in semileptonic decays of B or D mesons. As pointed out in Ref. [73], more data are needed to fix the LECs in the NNLO potentials. Since the dispersive calculation of Dφ andBφ scalar form factors depend on the scattering amplitudes of these systems, the LQCD results for the form factors can be used to mitigate the lack of data and help in the determination of the new unknown LECs.
One might also try to extend the MO representation to a formalism in a finite volume with unphysical quark masses, such that comparisons to the discretized lattice data could be directly undertaken. On the other hand, the chiral matching of the form factors can be carried out at higher order to take into account the expected sizable corrections in SU (3) HMChPT. Moreover, this improved matching will in practice suppose to perform additional subtractions in the dispersive representations of the form factors, and it should reduce the importance of the high-energy input used for the Hφ amplitudes. The high energy input turns out to be essential to describe the scalarB → π form factor near q 2 = 0, and it represents one of the major limitations of the current approach.
Both improvements would lead to a more precise and model-independent determination of the CKM matrix elements related to the heavy-to-light transitions.  There is a total of 22 degrees of freedom and the best-fit gives χ 2 /dof = 3.7. The LECs β B 1 , β B 2 and α2,i are given in units of GeV, GeV −1 and GeV −4 , respectively.
where α 2 , together with β B 1,2 , are fitted to the merit function defined in Eq. (49). (Note that α 0,1 are still expressed in terms of β B 1,2 .) Best fit results are compiled in Table VI