Study of $B\to\pi\ell\nu_{\ell}$ and $B^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ decays and determination of $|V_{ub}|$

We reassess the $B\to\pi\ell\nu_{\ell}$ differential branching ratio distribution experimental data released by the BaBar and Belle Collaborations supplemented with all lattice calculations of the $B\to\pi$ form factor shape available up to date obtained by the HPQCD, FNAL/MILC and RBC/UKQCD Collaborations. Our study is based on the method of Pad\'{e} approximants, and includes a detailed scrutiny of each individual data set that allow us to obtain $|V_{ub}|=3.53(8)_{\rm{stat}}(6)_{\rm{syst}}\times10^{-3}$. The semileptonic $B^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ decays are also addressed and the $\eta$-$\eta^{\prime}$ mixing discussed.


Introduction
Quark flavour-changing transitions in the Standard Model are described by the Cabibbo-Kobayashi-Maskawa (CKM) matrix whose elements, V ij , weight the strength of the interaction. The CKM matrix satisfies unitarity imposing i V ij V * ik = δ jk and j V ij V * kj = δ ik . To verify these relations a precise determination of the magnitude of the CKM elements becomes of capital importance since an eventual deviation of unitarity of the CKM matrix may be a hint of new physics. The most common (unitarity triangle) combination to look at is which contains the best-known side quantity V cd V * cb but also involves V ub , one of the least-known elements. The inclusive, B → X u ν , and exclusive, B → π ν , semileptonic decays of a B meson represent an advantageous laboratory to determine the value of |V ub | yielding the most precise value up to date. Inclusive determinations are based, for example, on the Operator Product Expansion and perturbative QCD while exclusive determinations require knowledge of the shape of the participant meson Form Factor (FF) as a function of q 2 , yielding the hadronic transition. Numerically, the 2015 PDG reported values showed a 3.1σ deviation between the inclusive, |V ub | = (4.41 ± 0.15 +0. 15 −0.17 ) × 10 −3 , and the exclusive, |V ub | = (3.28 ± 0.29) × 10 −3 [1], determinations with a resulting average of |V ub | = (4.13 ± 0.49) × 10 −3 . The updated 2016 PDG version [2] reports, respectively, |V ub | = (4.49 ± 0.16 +0. 16 −0.18 ) × 10 −3 and |V ub | = (3.72 ± 0.19) × 10 −3 [3] for the inclusive and exclusive decays whose deviation, 2.6σ, has slightly been reduced due to the one-σ shift of the exclusive result. At present, the PDG reports [4], respectively, |V ub | = (4.49±0.15 +0. 16 exp−0.17 th ±0.17)×10 −3 and |V ub | = (3.70 ± 0.10 ± 0.12) × 10 −3 [5] for the inclusive and exclusive determinations. The origin of this long-standing discrepancy between the inclusive and exclusive determinations still remains unclear, demanding the resulting combined average, |V ub | = (3.94 ± 0.36) × 10 −3 [4], to be borrowed with caution. As pointed out already in Ref. [6], and recently adopted in Ref. [7], a new physics explanation of this tension is very unlikely and, therefore, it might be due to an underestimation of uncertainties in the experimental and/or theoretical analysis.
In this work, we reexamine the exclusive B → π ν decays and extract |V ub | following an alternative approach, regarding the parameterization of the participant vector form factor, slightly different than the most commonly used z-expansion and Vector Meson Dominance models, and profiting from the large set of experimental data and lattice simulations. A detailed scrutiny of each individual data set, explored bin by bin, allows us to identify agreements and tensions among them and propose a path towards further determinations.
The hadronic matrix current for the B → π ν decay can be written as where q = p B − p π = p + p ν is the transferred momentum to the dilepton pair while f + (q 2 ) and f 0 (q 2 ) are, respectively, the participant vector and scalar form factors encoding the dynamics of the strong interactions occurring in the heavy-to-light B → π hadronic transition. For light leptons (e and µ), one can safely take the m → 0 limit so that only the f + (q 2 ) is relevant and the corresponding partial decay width distribution is given by 1 where p π = (m 2 B + m 2 π − q 2 ) 2 − 4m 2 B m 2 π is a kinematical factor accounting for the momentum of the pion in the B meson rest frame. The main source of uncertainty in the extraction of |V ub | lies on the vector form factor which, in turn, requires a reliable parameterization in terms of q 2 .
From the experimental side, the CLEO Collaboration reported the first measurement of the B → π ν branching fraction in 1996 [9], later updated in 2003 [10], and released the partial branching ratio distribution measured in 4 q 2 bins in 2007 [11]. More recently, the q 2 decay spectra have been measured, respectively, in 6-and 12bins of q 2 by BaBar in 2011 [12] and 2012 [13], and by Belle in 13-bins in 2011 [14] and in 13-and 7-bins for the B 0 and B − mode, respectively, in 2013 [15].
On the lattice QCD side, results on the form factor shape at large q 2 were obtained by the HPQCD Coll. in 2007 [16] and by the FNAL/MILC Coll. in 2008 [17] in 5-and 12-bins of q 2 , respectively. In 2015 the RBC/UKQCD Coll. released new results in 3 bins of q 2 [18] and the FNAL/MILC Coll. presented an updated analysis [3].
In total, we have a set of five experimental measurements of the B → π ν decay spectra driving the form factor shape at small q 2 and a set of four lattice QCD simulations for the form factor dominating the large q 2 region. In order to determine |V ub | with good precision (beyond 10%), it is desirable to have a suitable parameterization of the intermediate energy region (15)(16)(17)(18)(19)(20) GeV 2 ) connecting both the small-and large-q 2 regions in a continuous and derivable way, under the constraints of unitary and analyticity.
From a theoretical perspective, parameterizations based on resonance-exchange ideas [19,20,21,22] have been widely used so far to describe the B → π FF shape. The parameterizations proposed by Bećirevic-Kaidalov [23] and Ball-Zwicky [24], incorporating some properties of the FF such the value of the kinematical constraint at q 2 = 0 and the position of the B * pole in the old-times spirit [19,20,21,22], became rather popular in the first decennial of this century. Both descriptions contain free parameters, such additional poles that pick up effects of multi-particle states, to be fixed from fits to experimental data. However, the election of these ansätze induces a source of theoretical (or systematic) uncertainty difficult to quantify. Moreover, as argued in Ref. [17], if the reconstruction of the FF obtained only from fits to experimental data is seen inconsistent with the shape derived by the lattice Collaborations, one would not unveil whether experiment and theory disagree or simple parameterizations are insufficient. To improve on that, the so-called z-parameterization was proposed [26,27,28]. This is based on a conformal transformation expansion which guarantees unitary constraints on its coefficients, even though in practice the constraints are rather weak.
Lets us return to the 2015 and 2016 PDG editions reported values for |V ub | from exclusive processes, (3.28 ± 0.29) × 10 −3 and (3.72 ± 0.19) × 10 −3 , respectively. They were obtained from simultaneous fits to the four most precise measurements of BaBar [12,13] and Belle [14,15] together with the 2008 and 2015 MILC Collaboration lattice simulations on the FF, respectively. While the 2015 PDG value corresponded to the determination provided by the HFAG as of summer 2014 [1], the updated 2016 PDG version reports the value obtained by the MILC Collaboration in 2015 [3]. These two values have been determined by using a z-parameterization as fit function and show a sizable deviation of 1.3 σ whose origin stems mainly from the following fact. While the lattice FF simulations from the MILC Collaboration in 2008 [17] were included into the fit in the HFAG analysis of 2014, and on top of that only 4 of the 12 points were used to avoid correlations between neighboring points, the result obtained by MILC in 2015 considers their updated FF simulation ones [3]. Moreover, while the 7 bins of the B − decay mode measurement reported by Belle in 2013 were not included into the HFAG 2014 fit, the MILC 2015 include them into their analysis.
Previous PDG reported values, e.g. |V ub | = (3.23 ± 0.30) × 10 −3 in PDG 2012, showed the corresponding HFAG fit results obtained from simultaneous fits to the existent experimental measurements at the time together with the MILC form factor predictions of 2008 using 6 of the 12 points instead of 4 of 12 as in the HFAG result of 2014. Both the choice of the MILC 2008 number-and bin-points to fit and the omission of the HPQCD form factor lattice simulations of 2007 is not rather clear to us. Accepting the FNAL/MILC lattice form factor calculation of 2015 presents several improvements with respect to their 2008 predictions, 2 still the theoretical error associated to the FF represents the largest uncertainty in |V ub |. In this respect, the lattice simulation provided by the RBC/UKQCD Collaboration [18] has been welcomed, obtaining |V ub | = (3.61 ± 0.32 stat+syst ) × 10 −3 from a combined fit to their results for the form factor together with BaBar and Belle experimental data.
In 2016, the FLAG working group reported |V ub | = (3.62 ± 0.14) × 10 −3 from a fit to lattice and BaBar and Belle experimental data [35] while the current PDG edition reports the HFLAV value of 2017 |V ub | = (3.70 ± 0.10 ± 0.12) × 10 −3 [4,5] obtained from an averaged q 2 spectrum of all BaBar and Belle data sets constraining the χ 2 minimization by averaged values for the coefficients of the form factor parameteriza-tion derived by the lattice groups and by the LCSR prediction at q 2 = 0.
Finally, a closer look to the plots of the corresponding fit results of the different analyses reveals a discrepancy between HFAG 2014 and Ref. [29], and lattice groups [3,18] on the position of the last experimental datum of both BaBar 2011 and BaBar 2012 measurements 3 .
Although the different |V ub | determinations are consistent with each other, we find the situation slightly unclear and without consensus among different groups regarding the use of experimental and theoretical data to fit.
The main purpose of this work is to reanalyze the B → π ν experimental data and discuss the impact of including into the fit each of the lattice-QCD simulations on the FF shape. We use the method of Padé approximants (PAs in what follows) to parameterize the B → π transition. These provide for a model-independent method, simple and user-friendly, with the important advantage of incorporating FF's unitary and analyticity constraints by construction, thus providing a systematic error.
We have discussed the Padé method in Refs. [36,37,38] and illustrated its usefulness as fitting functions in Refs. [39,40,41,42] applied to the description of the π 0 , η and η transition form factors. In these cases, the approximants showed an interesting ability to connect the low-and high-energy realms while improving the description of part of the intermediate-energy regime. The method allows us here to obtain a value for |V ub |, including both statistical and systematic uncertainties coming from the fit function, with a stamp of model independence. Constraints from unitary of the form factor will show up naturally and will provide for a roadmap towards next steps to follow both for theoretical as well experimental studies.
Although being the most precise, B → π ν only amounts to ∼ 7% of the B → X u ν decays. Measurements of the branching fractions distributions of B + → ω + ν and of B + → η + ν in 5 bins of q 2 were released in 2012, and the branching ratio of B + → η + ν reported, by the semileptonic charmless program of BaBar [13]. In the second part of this work, we will tackle the B + → η ( ) + ν decays, predicting the differential branching ratio distributions and extracting the η-η mixing angle taking advantage of the B → π form factor parameterizations obtained in the first part of this work.
As a final introductory remark, we shall mention that a method based on dispersion theory to extract |V ub | from the B 4 decay has been proposed in Ref. [43].
This article is structured then as follows: In section 2 we address the analytical structure of the participant B → π form factor and discuss the most common theoretical descriptions that have been considered in literature so far. In this section we also present our proposal, a parameterization based on the unitary and analyticity of the FF which allows us to use a sequence of PAs. In section 3.1, we show our fit results to the BaBar, Belle and CLEO differential branching ratio distribution experimental data which enables us to determine the product |V ub f + (0)| and extract, subsequently, |V ub | by using the LCSR prediction of f + (0) given in Ref. [30]. In section 3.2 we discuss the impact of including the different lattice QCD predictions on the FF shape into the analysis and determine |V ub | directly from a simultaneous fit.
In this section we present our central fit results, evaluate the role of introducing the value of f + (0) as an additional restriction in the χ 2 minimization, and perform fits to the lattice data alone. Unitary constraints on the Padé approximants are discussed in section 3.3. In section 4 we predict the B + → η ( ) + ν differential branching fractions distributions and determine the η-η mixing. Finally, our conclusions are devoted to section 5.
Preliminary results of this study have been presented in Refs. [31,33].

B → π form factor
A form factor is an analytic function everywhere in the complex plane except for isolated poles and branch cuts. Poles correspond to single particle intermediate states while branch cuts originate when the energy reaches a threshold for producing multiparticle intermediate states. For the B → π ν decay concerning us, the lightest production threshold is located at s th = (m B + m π ) 2 GeV 2 , lying slightly above the available kinematical energy range of the decay, 0 < q 2 < (m B − m π ) 2 GeV 2 . A first approximation to the form factor suggests a single pole description driven by the exchange of aūb intermediate state, the B * meson with mass m B * = 5.325 GeV (with very small width) and quantum numbers J P = 1 − . For illustrative purposes, let us consider a dispersive representation of the form factor in terms of q 2 , where q 2 is the invariant mass of the lepton pair, with Imf + (s) = πρ(s). The single pole description in Eq. (4) would correspond to using for the spectral function ρ(s) = f + (0)m 2 B * δ(s − m 2 B * ). This give raise to the Vector Meson Dominance model (VMD) with a B * pole appearing between the available phase space and the lowest production threshold, (m B − m π ) 2 < s p < (m B + m π ) 2 GeV 2 , where f + (0) is a normalization constant. However, this model obviates effects of heavier vector states. Bećirević and Kaidalov (BK) [23] proposed a modification of the VMD via including, above q 2 = (m B +m π ) 2 GeV 2 , a heavier narrow-width resonance, a B * , through adding Imf Implementing that the form factor behaves as 1/q 4 at large q 2 together with f + (0) = r 1 + r 2 , the standard expression for the BK form factor reads where α fixes the position of the second fitted effective pole. Later on, Ball and Zwicky (BZ) [24,25] proposed a similar expression in terms of three parameters {f + (0), r, α} by imposing the form factor to fall-off as ∼ 1/q 2 at large q 2 instead. The matching f + (0) = r 1 + r 2 and r = r 2 (α − 1) leads where r may be understood as a parameter which encodes the relative weight of the second effective resonance with respect to the first one. The above two parameterizations fix the position of the B * pole to its mass, m B * = 5.325 GeV, while the rest of free parameters, {f + (0), α} and {f + (0), r, α} in Eqs. (7) and (8), respectively, are inferred from fits to experimental data.
Exploiting the analyticity and positivity properties of the vacuum polarization functions, Okubo and collaborators proposed the method of unitary bounds [26] in the context of kaon decays, which later on was applied for semileptonic B decays [28,44]. This method, called z-parameterization and reviewed in Refs. [27,45], parameterizes f + (q 2 ) as a Taylor expansion in terms of a conformal complex variable z as follows: where with t + = (m B + m π ) 2 GeV 2 and φ(q 2 , q 2 0 ) a function given in Ref. [27]. The function P (q 2 ) = z(q 2 , m 2 B * ) is the Blaschke factor which accounts for the pole at q 2 = m 2 B * . The free parameter q 2 0 is chosen to optimize the fit. Assuming the spectral function driving the FF to be saturated by Bπ vector intermediate states, unitary and crossing symmetry guarantee the coefficients a n (q 2 0 ) to satisfy ∞ n=0 a 2 n (q 2 0 ) ≤ 1. In practice, Eq. (9) is truncated at a finite order (typically up to first or second order) which implies the FF to behave as ∼ 1/q 4 at large |q 2 | due to φ(q 2 , q 2 0 ), in contradiction with perturbative QCD scaling [21,22]. Beyond, as discussed in [28], the outer function has an unphysical singularity at the Bπ production threshold t + . This unphysical singularity may distort the behavior near the upper end of the physical region, where the FF is poorly known. These considerations triggered an alternative z-parameterization proposed in [28] by Bourrely-Caprini-Lellouch (BCL): where the pole included by hand ensures the correct analytic structure in the complex plane and the proper scaling, f + (q 2 ) ∼ 1/q 2 at large q 2 . Let us comment that the z-parameterization is not a zero-preserving transformation with respect of q 2 unless the particular choice q 2 0 = 0 is made, which implies z → 0 does not come from q 2 → 0, but rather from a large q 2 value. This poses a word of caution when using the z-parameterization to determine the behavior of the FF at low q 2 . We shall add here that the definition of z(q 2 , q 2 0 ) corresponds formally with a Quadratic approximant, a well-defined extension of a PA that includes square-root terms [31,32] 4 . As such, it is formally a PA of order given by the truncated series, either in Eq. (9) or (11), for which the PA convergence constraints must be applied.
To complete the overview, the AFHNV approach [46,47] is based on the Omnès representation which expresses the analytic function in terms of its phase along the boundary of the analyticity domain. If one takes into account the pole q 2 = m 2 B * , assumes that the FF have no zeros in the complex plane, and by Watson's theorem that the phase δ(t) is equal, below the first inelastic threshold, to the phase of the P -wave with I = 1/2 of the πB → πB elastic scattering, the representation reads (assuming n 1 which implies a multiply-subtracted dispersion relation-and neglects altogether the dispersive integral): where . Notice that after the multiply-subtracted dispersion relation the exponential behavior at large q 2 does not corresponds to the one from QCD and, due to the lack of the dispersive integral, the original branch cut q 2 ≤ t + .

Our proposal: Padé approximants
The form factor f + (q 2 ) is a Stieltjes function, which is a function that can be represented by an integral form defined as [48] f + ( where φ(u) is any bounded and non-decreasing function. By defining R = (m B +m π ) 2 GeV 2 , identifying dφ(u) = 1 π Imf + (1/u) u du, and making the change of variables u = 1/s, Eq. (13) returns the dispersive representation of the form factor given in Eq. (4). To be strict, Eq. (4) is a meromorphic function of Stieltjes type.
Since the FF, and its imaginary part, is created by the vector current, Imf + (s) is a positive function, the requirement of φ(u) to be non-decreasing is fulfilled and the convergence of PAs to the FF is guaranteed.
Padé Theory not only provides a convergence theorem for a sequence of PAs to Stieltjes (or Stieltjes-type) functions, i.e., lim N,M →∞ P N M (q 2 ) − f + (q 2 ) = 0, but also its rate of convergence [48,49], given by the difference of two consecutive elements in the PA sequence. As we will see later, this error prescription will return very small theoretical uncertainties. To be more conservative, in Refs. [37,38,39,40] we designed a different method to extract such uncertainty which yields errors at the level of the statistical ones. 4 In particular, z(q 2 Padé approximants to a given function are ratios of two polynomials (with degree M and N , respectively) 5 , with coefficients determined after imposing a set of a accuracy-through-order conditions with the function one wants to approximate, which is f We would like to remark that Eqs. (5), (7), (8), and (9) can be seen as particular elements of the general sequence of PAs given in Eq. (14).
Besides ordinary sequences of PAs we will also consider Padé Type approximants T M N (q 2 ) and Partial Padé approximants P M Q,N −Q (q 2 ) in our study. T M N (q 2 ) have the denominator fixed in advanced (by imposing the location of the zeros of it), while P M Q,N −Q (q 2 ) only Q zeros of the denominator are fixed in advanced while the rest are left free. Strictly speaking then VMD, BK and BZ correspond, respectively, to the T 0 1 (q 2 ), P 0 1,1 (q 2 ) and P 1 1,1 (q 2 ) elements, while the z-parameterization corresponds to a polynomial in terms of a Q 1 1,1 as we argued before. The main advantage of fixing a pole is that the number of parameters to fit decreases by one and typically allows to reach higher elements of the sequence [37]. If the sequence is large enough and the position of the first singularity is accurately known, the convergence of the T M N (q 2 ) is faster than the convergence of ordinary PAs for Stieltjes functions [36,48,50].

Fits to the B → π ν BaBar and Belle data
Our first analysis consists of fitting the most recent B → π ν branching ratio distribution experimental data released by BaBar in 2011 [12] and 2012 [13] and by Belle in 2011 [14] and 2013 [15]. We will also briefly discuss the effect of including CLEO 2007 results [11] into the fit, which are usually neglected. In order to facilitate the reproduction of our results, we would like to write down from which Tables of the papers the experimental data we use come from. For CLEO 2007 we use the results reported in Table I Tables  III, IV and V of Ref. [14], respectively. Finally, for Belle 2013 we employ the data given, respectively, in Tables XVII, XVIII, XIX and XX of Ref. [15]. To this later data we have added, as suggested in Table XII of [15], a systematic uncertainty of 5.0% and 5.1% of the q 2 bin value for the B + and B 0 mode, respectively, and assumed 5 With any loss of generality, we take b 0 = 1 for definiteness. 6 CLEO 2007 result consists in measurements of partial branching fractions in only 4 unequal q 2 subregions and no bin-to-bin correlation matrix is reported. For our analysis, we have placed each experimental datum at the middle of each of the corresponding subregions and scaled the bin values accordingly. 7 For BaBar 2012, we use the combined analysis of both B 0 and B − modes assuming isospin symmetry. a systematic correlation of the 49% between the two modes as written in the paper below that Table. For our study, we assume, for convenience, isospin symmetry to translate the Belle 2013 data on the B − mode to the B 0 ones through where τ B 0 = (1.520 ± 0.004) × 10 −12 s and τ B − = (1.638 ± 0.004) × 10 −12 s, are, respectively, the mean life time of the neutral and charged B mesons [2]. In all, we will treat the five experimental data sets as independent measurements, i.e., no statistical either systematic correlations between the five different analyses is considered [3,18]. The χ 2 function minimized in our first fit is defined as where with Γ B the full width of the B meson, Cov ij denotes the corresponding covariance matrix and P M N (q 2 ) = P M N (q 2 )/P M N (0) the PA normalized to unity at the origin of energies whose coefficients will be determined by the fit.
We start fitting with ordinary PAs of the type P M 1 (q 2 ) and P M 2 (q 2 ) where the poles are left free to be fitted and we reach M = 2 and M = 0, respectively. Then, we proceed to fit with sequences of the type T M 1 (q 2 ) and P M 1,1 (q 2 ) by fixing the B * pole to m B * = 5.325 GeV reaching, respectively, M = 2 and M = 1. In Fig. 1, we provide a graphical account of the fit as obtained with P 2 1 (q 2 ) compared to data while our fit results are collected in the first row of Table 1. From the plot we observe that the uncertainty associated to the fit, given by the gray error band, is slightly larger in the low-q 2 energy region while from the Table we read that the values for |V ub | determined with approximants with two poles i.e. P 1 2 (q 2 ) and P 1 1,1 (q 2 ), give identical results than the single pole ones, P 2 1 (q 2 ) and T 2 1 (q 2 ). Then, we add CLEO 2007 experimental data into the χ 2 minimization of Eq. (16) and report the corresponding fit results in the second row of the Table 1, Figure 1: Simultaneous fit to BaBar [12,13] and Belle [14,15] B → π ν experimental data as obtained from the χ 2 data minimization of Eq. (16) with a P 2 1 (q 2 ) approximant (black solid line). CLEO data [11] is not included from the fit and rather shown for comparison.
comparison with the results shown in the first row we conclude that the effect of including this data into the fit is tiny.
In order to further improve on what can be learned from experimental data, we have also fitted experimental data of each Collaboration separately, an exercise that will be very illustrative in order to determine |V ub |.
The individual fits are displayed in Fig. 2 and the corresponding results shown, respectively, in the third (CLEO07), fourth (BaBar11), fifth (BaBar12), sixth (Belle11) and seventh (Belle13) rows of Table 1. From our set of fits collected in this Table, the diagonal and near diagonal P 2 2 (q 2 ), P 1 1,(1) (q 2 ) and P 2 1,(1) (q 2 ) approximants deserve special attention. For these approximants we find some extraneous poles placed either far away from the origin (marked with † in the Table) or pair up with a close-by zero in the numerator becoming what is called a defect or Froissart doublet (marked with † † in the Table), in accordance with the Nutall-Pommerenke's convergence theorem [36,48]. We would like to point out that the zeros of the numerator when individual fits to the BaBar 2012 data are performed tend to lie within the radius of convergence in the region of negative q 2 , a region we expected out of zeros. This feature may explain why the corresponding distribution is more rounded and with a sizable negative fall off at the origin in comparison with the other three individual fits. We also note that a Froissart doublet doublet appears at q 2 = −1 GeV for the   The product |V ub f + (0)| as obtained from fits to B → π ν data depending if the B * pole is let as a free parameter to fit or fixed at m B * = 5.325 GeV. The corresponding |V ub | value extracted using f + (0) = 0.261 +0.020 −0.023 [30], the pole(s) of the approximants and the χ 2 dof are also shown. Poles placed far away from the origin and Froissart doublet are denoted by † and † † , respectively. Errors are only statistical and symmetrized in the last column.
that a PA to a Stieltjes function is also a Stieltjes function as well [48]. As such, it must have an imaginary part positive defined. This feature does not correspond with what we obtain in our fits and disagrees with what we expect from a Stieltjes function. All zeros and poles of our approximants must lie along the unitary branch cut in order to fulfill the unitary requirements that the FF imposes. If a particular PA does not show this feature means the set of data fitted is not fulfilling the unitary requirements that must have. Thus, both defects and the appearance of poles and zeros outside the unitary branch cut are indications of a violation of unitary to a certain degree. We shall come back to this point later (see section 3.3).
We would like to notice that individual fits to CLEO data lead unrealistic results but for P 0 1,1 (q 2 ). Also notice that the fits to BaBar11 experimental data lead the worst χ 2 /dof, in agreement with Ref. [3], and the largest values for |V ub |, in line with Ref. [30] but in contradiction with Refs. [3,18]. These two features are somehow reflected in the left-top panel of Fig. 2 both in the error band and in the value of the branching ratio distribution at q 2 = 0 which are, respectively, wider and larger than in the other three panels of the figure. On the contrary, the fits to the BaBar12 data gives the best χ 2 /dof and tends to give smaller |V ub | values.
From each of the individual fits shown in Table 1 we can order the experimental Collaborations according to their bottom-up |V ub | values as: BaBar12, Belle11, Belle13 and BaBar11. This ordering is in line with the corresponding |V ub | values reported by the experimental groups from fits to their own experimental data.
The neat effect of fitting all experimental data sets together with respect to fitting data of each Collaboration separately can be seen in Fig. 3, where we represent the number of σ deviations of each experimental datum with respect the corresponding fits. In this figure, markers given by solid geometric figures accounts for the fit as given in Fig. 1 while empty geometric figures stand for the fits as shown in Fig. 2. This allows us to order the four experimental data sets according to their increasing degree of soundness with respect to the common fit i.e., BaBar11, Belle13, Belle11 and BaBar12. Clearly, BaBar11 data points suffer the largest deviation when including the other sets of data into the fit (see left-top panel in Fig. 3) while, on the contrary, BaBar12 and Belle11 data points seem to drive the χ 2 minimization dominating the fit (see right-top and left-down panels, respectively, in Fig. 3). Regarding Belle13 experimental data points, they show some oscillatory scatter lying in between BaBar11 and BaBar12/Belle11 cases.

Incorporating form factor lattice calculations
In the previous section we have not accessed the description of the form factor but rather its normalized version to unity at q 2 = 0. In order to achieve a parameterization of the form factor we include the form factor shape predictions at large q 2 obtained on the lattice as new data sets to be fitted. In particular, we consider  [18] and, finally, the updated analysis of the FNAL/MILC Collaboration of 2015 [3] 9 . The main advantage of performing a si- 9 While the FF predictions as obtained by HPQCD2007, MILC2008 and RBC/UKQCD in Refs. [16,17,18], respectively, are publicly available and the corresponding results given in the papers, the updated 2015 results of MILC are not. However, we have generated the FF from the fit as given in Table XIV of Ref. [3]. For the sake of comparison with their former 2008 predictions, we have generated 12 data points placed at the same q 2 -bins. For the extracted data points, the interested reader may contact the multaneous fit to all measured q 2 spectra experimental data supplemented by lattice QCD results on the FF shape is that not only |V ub | but also f + (0) can be determined directly from the fit since lattice data drive the height of the curve of the decay spectra. The χ 2 function to be minimized reads where χ 2 data corresponds to Eq. (16) and σ i is the uncertainty corresponding to the i-th bin.
The results derived from the minimization of Eq. (18) are collected in Table 2. In contrast to Table 1, Table 2 collects the results as obtained with each element of corresponding authors. We would like to thank Elvira Gámiz for correspondence along these lines.  the corresponding sequences going up to P 2 1 , P 2 2 , T 2 1 and P 2 1,1 , respectively. The final results, given in the last column, include both statistical, from the fit, and systematic uncertainties from the truncated PA sequence as the difference of central values of the element we have stopped the sequence and the preceding one. Notice that the systematic uncertainty increases when the B * pole is fixed.
The impact of including lattice data into the fit is evident and allow us to determine |V ub | with improved precision reducing the associated statistical uncertainty by ∼ 80% with respect to the case when only the decay spectra is fitted (cf. Table 1). In addition, the sequence P L 2 has been enlarged by one element. Again, we find some extraneous poles for the diagonal P 2 2 and P 2 1,1 elements. In the former, we find a complex-conjugate (c.c.) pole with a small imaginary part (see the dedicated discussion in section 3.3) while in the latter we find that one pole tends to pair up with a close-by zero in the numerator.
As a matter of example, we gather the coefficients of the Padé approximant P 2 1 (q 2 ) in Table 3. In this Table we also provide the series coefficients b (n) + corresponding to the BCL parameterization (cf. Eq. (11)) that are obtained by matching the Taylor series expansion of P 2 1 (q 2 ) with the power series expansion of the BCL parameterization at O(q 4 ). The coefficients thus obtained are not directly fitted to data but rather reconstructed from our rational function. These lie in the ballpark of the most recent RBC/UKQCD and FNAL/MILC lattice determinations [18,3] shown, respectively, in the fifth and sixth columns of the Table and are seen in nice agreement with the HFLAV fit values [5] given in the last column. A graphical account of the corresponding P 2 1 (q 2 ) combined fit result is depicted in Figs. 4 and 5 as compared to the decay spectra and FF lattice data, respectively. In the latter, our prediction for the BCL parameterization is also shown (purple dashed curve), accommodating pretty well all lattice data but the last datum and seen in nice agreement with the P 2 1 (q 2 ) (black solid curve) element it is reconstructed from. A closer look to the FF shape displayed in Fig. 5 reveals that the lattice simulations derived by the FNAL/MILC Collaboration in 2015 seem to dominate the large q 2 region (cf. Table 5).
Our preferred values for |V ub | and f + (0) after the simultaneous fit results shown in Table 2 is corresponding to P 2 1 (q 2 ) when the pole is let as a free parameter to fit. This choice is based on the fact that the second pole of the sequence of the type P M 2 is rare indicating that the single pole behavior for the form factor seems favored.
To compare on the same footing regarding the number of free parameters we choose |V ub | = 3.53(8) stat (5) syst × 10 −3 , f + (0) = 0.264(10) stat (5) syst , (20) corresponding to the partial Padé P 2 1,1 when the B * pole is fixed. Notice that the corresponding systematic uncertainties are large enough to cover the difference with P 2 2 and T 2 1 , respectively. The impact of including the value of the FF at q 2 = 0, f + (0) = 0.261 +0.020 −0.023 [30], as an external restriction in the χ 2 , Eq. (18), is probed through the fits displayed in Table 4 for our best fit sequences discussed previously when f + (0) was not included in the minimization (cf. Table 2), and repeated in this Table 4, second column, for ease of comparison. The corresponding fits are almost identical, which guarantees the independency of our results with respect of the model calculation of f + (0).
We have also performed fits including CLEO 2007 data [11] in the χ 2 , and found that their impact in the global fit is marginal and we hence refrain to show them, and explored the effect of fitting all experimental data together with certain groups of lattice FF simulations. We have considered three groups, HPQCD+RBC/UKQCD, HPQCD+RBC/UKQCD+MILC08 and HPQCD+RBC/UKQCD+MILC2015, and collected the results, respectively, in the second, third and fourth columns of Table 5 for the P 2 1 (q 2 ) (the other sequences yield almost identical results). While the second and third columns lead similar results, the fourth column, including the updated FNAL/MILC form factor simulation of 2015, clearly shifts upwards(downwards) the |V ub |(f + (0)) value by about 1.3σ yielding smaller statistical uncertainties and slightly enlarging the χ 2 dof .
Fits to BaBar and Belle data [12,13,14,15] Table 2: |V ub | and f + (0) values as obtained from a simultaneous fit to B → π ν decay data and lattice QCD form factor simulations. The pole(s) of the approximants and the χ 2 dof are also shown. Poles placed far away from the origin and Froissart doublet are denoted by † and † † , respectively, while c.c. stands for a complex-conjugate pole with a small imaginary part. The results in the last column include a systematic error coming from the difference of central values of the last two elements of the corresponding PA sequences. The errors are symmetrized.  Table 3: Coefficients of the Padé approximant P 2 1 (q 2 ), with the pole let as a free parameter, and of the reconstructed BCL parameterization, where the pole is fixed to the B * . The latter are compared with the fitted coefficients determined by the RBC/UKQCD and FNAL/MILC lattice groups [18,3] and with the HFLAV results [5].  Figure 4: Differential branching ratio distribution for B → π ν decays as obtained from a combined fit to experimental data and lattice predictions on the form factor shape with a P 2 1 (q 2 ) (black solid curve). CLEO data [11] is excluded from the fit but rather shown for comparison.
Upon comparison with last column, we conclude that FNAL/MILC simulations of 2015 drives the form factor while disagreeing for more than 1σ with respect to all the other lattice simulations (including their own determination but from 2008). This fact explains why the 2016 PDG reported value has been shifted with respect to the earlier edition by +1σ.
We close this section by performing fits to the lattice FF predictions alone and extract f + (0). This kind of exercise is new and, as a byproduct, allows us to determine |V ub | by equating the corresponding expression for the branching ratio (BR) to the measured ones, BR(B 0 → π − + ν ) = (1.45±0.05)×10 −4 [4]. We only obtain reliable results when, at least, the HPQCD 2007 and MILC 2015 predictions are included into the data sets to be fitted and for approximants with two poles. The corresponding fit results are gathered in Table 6.

Unitary constraints on PA's fits
All zeros and poles of our approximants must lie along the unitary branch cut in order to fulfill the unitary requirements that the FF imposes [48].  Figure 5: B → π form factor as obtained from a combined fit to experimental data and lattice predictions on the form factor shape with the approximant P 2 1 (q 2 ) (black solid curve). Our prediction for the BCL parameterization is also shown (purple dashed curve). a certain degree. Since we have performed a dedicated analysis Collaboration by Collaboration, bin by bin, and since we have found some cases which slightly violate these two statements, mostly happening when the BaBar 2012 data is involved, we are able to identify the source of the unitary deviation.
We find either complex-conjugate poles with an small imaginary part or a zero(s) within the radius of convergence for the P 0 2 (q 2 ) and P 1,2 1 (q 2 ) elements when fitting individually the BaBar 2012 data set, respectively. In particular, the complexconjugate pole of P 0 2 (q 2 ) is found to be at 5.72 ± i0.53 GeV while the zeros of the numerator are placed at −4.88 GeV and −4.40 GeV for the P 1 1 (q 2 ) and P 2 1 (q 2 ) elements, respectively (the second zero of P 2 1 (q 2 ) is placed at 12.97 GeV, far away from the origin). A complex-conjugate pole with an small imaginary part also shows up in the P 2 2 (q 2 ) element when performing the joint fit to data and lattice. In order to further explore on the origin of these extraneous poles and zeros we have also performed fits removing one experimental datum of each Collaboration e.g. those with more tension according to our Figures 2 and 3, and see what can we learn.
In particular, we remove the fifth datum of BaBar 2011, the tenth of BaBar 2012 and Belle 2011 and the bin located at 9 GeV 2 of Belle 2013. By doing this, we find that the zeros tend to move away from the radius of convergence while the complexconjugate poles become cancelled by a close-by zero in the numerator, a Froissart pole. Table 2 Constraining f + (0)
The impact of these four points is remarkable, inducing a positive shift on the |V ub | after removing them by about ∆|V ub | = 0.05 × 10 −3 , a 0.4σ deviation.
Breaking of unitary is then reducing the value of |V ub | an enlarging the discrepancy between inclusive and exclusive determinations. In view of this fact and the difficulty on deciding the best strategy to take this unitary violation into account when dealing with experimental data (other strategies beyond removing bins could be envisaged), we have decided to add in quadrature the difference ∆|V ub | as an extra source of error in our final determination of the CKM parameter and Eq. (19) should be sound |V ub | = 3.53(8) stat (6) syst × 10 −3 . This error could be removed as soon as the experimental Collaborations could take our observation into account and explore systematically the potential unitary violation within their data sets.

B + → η ( ) + ν decays and η-η mixing
In the previous section, the B → π form factor f + (q 2 ) has been parameterized using PAs to fit experimental data on the B → π ν differential branching ratio distribution w/o lattice FF simulations. In this section, we would like to take advantage of these parameterizations to describe the B + → η ( ) + ν decays as discussed in the following.
The expression for the differential B + → η ( ) + ν decay width is given by the same expression as for the B → π ν decay mode in Eq. (3) by replacing the final state pion by η ( ) where now f B + η ( ) + (q 2 ) represents the hadronic B + → η ( ) transition. What the B + → η ( ) transition is probing is the light-quark content of the η ( ) mesons since the ss component can only be accessed via a B s meson decay. This is so because from the quark-flavour perspective, η ( ) mesons are an admixture of uū, dd and ss components. Defining |η q = 1 √ 2 |uū + dd and |η s = |ss in this quark-flavour basis, one can relate the mathematical |η q,s states with the physical |η ( ) ones through the following matrix rotation where φ gives the degree of admixture.
Contrary to f Bπ + (q 2 ), there are no, to the best of our knowledge, FF simulations of f B + η ( ) + (q 2 ) on the lattice while only a few calculations at q 2 = 0 exist [8,24,54].
Therefore, we will relate the f B + η ( ) + (q 2 ) with the f Bπ + (q 2 ) using the quark-flavor basis. Assuming isospin symmetry between the u and d quarks, the form factor f B + η ( ) + (q 2 ) can be related to the f B + π 0 + (q 2 ) ones through [55] taking η uū π 0 as in Refs. [40,56]. From Eqs. (21) and (23), and taking f Bπ + (q 2 ) from any of the descriptions given in Table 2 together with the corresponding values for |V ub |, we can describe the differential branching ratio distribution of the B + → η ( ) + ν decays by setting the numerical value of the η-η mixing angle to, for example, φ = 38.3(1.6) [40]. Our prediction for the B + → η + ν differential branching fraction distribution is shown and compared with BaBar 2012 measurements in 5 bins of q 2 in Fig. 6 for P 2 1 (q 2 ) (black solid curve). Our description is seen quite in accordance with data although the second experimental datum seems to be slightly in tension. In Fig. 6 we also show our prediction for the B + → η + ν branching ratio distribution (blue solid curve), in this case without any experimental data to compare with.
The corresponding results are collected in Table 7. We observe that the central values show some scatter though they all agree within errors due to the large uncertainties in Eq. (25). In order to extract the mixing angle φ from B → η ( ) transitions with more precision, measurements of these decays with higher precision are required.
As a final exercise, we would also liked to fit either the individual BaBar 2012 B + → η + ν decay experimental data or perform a combined fit to B + → π + ν and B + → η + ν decays with the goal to provide an alternative semileptonic charmless B decay determination of |V ub |. However, due to the poor experimental situation in the case of B + → η + ν , we decide to postpone this analysis for the future.

Conclusions
In this paper we have reexamined the B → π ν decays to extract the CKM parameter |V ub | based on experimental data, lattice calculations and unitary constraints of the participant form factor. Contrary to the most commonly used z-expansion and Vector Meson Dominance models, we perform our analysis based on the method of Padé Approximants after realizing that most of the recent previous analyses belong to the Padé Theory, even though no one mention it. Thus, the rules and constrains imposed by the convergence theorems for Padé Approximants to the form factor, so far neglected, are fully exploited here, allowing to ascribe to our final result a new source of systematic or truncation error.
From our dedicated analysis we obtain |V ub | = 3.53(8) stat (6) syst × 10 −3 . This quantity includes both statistical, from the fitted data, and systematic, from the truncation of the Padé sequence, uncertainties, and has been obtained guaranteeing the independency with respect of the model calculation of f + (0) as external constrain.
On a first stage, after a detailed review of the state-of-the-art experimental data, determinations of |V ub | and theoretical representations of the analytical structure of the form factor, we have analyzed the measured q 2 differential branching ratio distribution experimental data released by the BaBar and Belle Collaborations. Our fitting strategy started by performing a combined analysis to all data sets using different types of Padé sequences. We thus have determined first the product |V ub f + (0)| directly from the fits and then extracted the CKM element |V ub | by invoking external theoretical information on f + (0). The resulting fit results are presented in Table 1 and a graphical account provided in Fig. 1. We then have carried out a detailed analysis Collaboration by Collaboration. The outcome of the individual fits is displayed in Fig. 2 and the neat effect on each experimental datum due to fitting all experimental data together with respect to fitting data of each Collaboration separately is shown in Fig. 3. This exercise allow us to classify the four differential experimental data sets according to their increasing degree of robustness: BaBar 2011, Belle 2013, Belle 2011 and BaBar 2012.
On a second stage, we have included into the analysis the four available lattice QCD predictions on the form factor shape. This data dominates the large-q 2 region and it is essential for a precise determination of |V ub |. The corresponding fit results are collected in Table 2 indicating that the statistical uncertainty associated to |V ub | is reduced by ∼ 80% after the inclusion of lattice data. We have also found that, out of the four lattice form factor simulations, the predictions released by the MILC Collaboration in 2015 tends to drive the form factor (see Table 5) but slightly enlarging the χ 2 dof . As a byproduct of our analysis, we have predicted the BCL form factor series coefficients that are obtained by matching the corresponding Taylor series expansion. The coefficients thus obtained are shown and compared with the determinations given by lattice groups in Table 3 while the q 2 shape of the reconstructed BCL parameterization is displayed in Fig. 5 proving the ability of the Padé Approximants in this transition.
On a third stage, motivated by the impact of the lattice data, we have also explored fits to the lattice predictions alone. The fit results are shown in Table 6 reflecting that only those approximants with two poles have the ability to extract first f + (0) and then determine |V ub | by equating the theoretical expression for the branching ratio to the corresponding experimental measurement.
Our central result, |V ub | = 3.53(8) stat (6) syst × 10 −3 , is presented and compared with other determinations using other methods and fitted data sets in Fig. 7. We would like to remark two features concerning this value that are related to the use of Padé Approximants. The first one, is that the central value tends to fall slightly downwards with respect to the values determined with the z-expansion parameterization in the studies carried out in the recent years. And the second one, is that the method allow us to ascribe a systematic uncertainty from the truncated Padé sequence. In fact, the z-paramaterizations do also allow to attribute a systematic error following the same reasoning. However, in practice, it has not so far usually been considered. For example, based on our criterion, the result as obtained by the FNAL/MILC Collaboration in 2015 would read |V ub | = 3.72(16) stat (9) syst , where the systematic uncertainty stems from the differing results for N = 3, 4 (cf. Eq. (11)). In our study, the ascribed systematic uncertainty includes, for the first time, an additional conservative source of error due to the unitarity constraints discussed in section 3.3. These constraints have to do with the appearance of extraneous poles and zeros outside the unitary branch cut and might indicate, to a certain degree, violations of unitarity.
As a final concluding remark for the B → π ν decays, we would like to point out that, contrary to the z-expansion and VMD models where the B * pole position is fixed to 5.325 GeV in advance, a very competitive value for |V ub | can be extracted without imposing any information regarding the position of it as we have shown along the lines of our detailed analysis.
In the second part of this work, we have addressed the B + → η ( ) + ν decays taking advantage of the B → π form factor parameterizations derived in the first part. In particular, we relate the participant Bη ( ) form factor to the Bπ ones by a single Euler angle rotation assuming that the light-quark component of the η ( ) is a qq pion to a large extent. Under this simple assumption, we obtain a reliable prediction for the differential branching ratio distribution of the B + → η + ν decay as shown in Fig. 6 compared to the BaBar measurement in 5 bin of q 2 released in 2012. As a byproduct of our study, we have also extracted the η-η mixing angle. This quantity, however, carries a large statistical error due to the large uncertainty on the measured B + → η ( ) + ν branching ratios. Regarding our prediction for the B + → η + ν decay distribution, there is no experimental data to compare with so far. In order to go beyond the simple quark-flavour basis decomposition and extract the η-η mixing angle with improved precision we would like to encourage experimental groups to measure these semileptonic B + → η ( ) transitions with improved precision.  [53] and this work ( ), from B → ω ν ( ) and B → ρ ν ( ) Bharucha 2015 [58] and from Λ b → pµν µ (•) LHCb [59], and from indirect fits (•) UTFit 2016 [60] and CKMfitter 2015 [61]. The solid and dashed error bar account, respectively, for the statistical and systematic uncertainties.