Dispersive Analysis of B → K ( ∗ ) and B s → ϕ Form Factors

,


Introduction
Rare B decays are decays of B mesons with small branching ratios of the order of 10 −6 or less.In the Standard Model (SM) these tiny rates are explained by the fact that they are loopmediated processes with CKM or GIM suppression.The fact that experimental measurements of these rates are of the same order magnitude of the SM predictions already poses strong indirect constraints on Beyond-the-SM (BSM) physics, most prominently in models where these processes arise at tree level.This ability to pose strong constraints on BSM physics is the same as discovery power [1].
A class of rare B decays that has been discussed extensively in over ten years of LHC operations are the Flavour-Changing Neutral Current (FCNC) exclusive processes of the type b → sℓ + ℓ − , such as B → Kℓ + ℓ − , B → K * ℓ + ℓ − , and Bs → ϕℓ + ℓ − .While the importance of these measurements was clear and had been discussed before [2][3][4][5][6][7], a generalised interest was driven by the "anomalies" observed starting from 2013 [8][9][10][11][12][13][14][15][16][17][18][19].While the deviations in the observables measuring lepton-flavour non-universality seem to be now more in agreement with the SM, the deviations in the muonic observables are still puzzling.Independently of this, the original motivations to study these rare decays have only been strengthened by the sheer amount and the precision of the experimental measurements.These measurements will clearly be one of the legacies of the LHC era.
On the theoretical side, efforts are focused on the control of hadronic uncertainties that are crucial to the unambiguous interpretation of the data.All such exclusive b → sℓ + ℓ − observables can be computed from a handful of decay amplitudes given by (see, e.g., [20]) Here, q is the four-momentum of the lepton pair, L/R indicates the lepton chirality, λ is the dilepton polarisation, C i are Wilson coefficients of the Effective Field Theory (EFT) below the weak scale, and F (T ),λ (q 2 ), H λ (q 2 ) are respectively local and non-local hadronic form factors (FFs).These hadronic FFs depend on the initial and final states and on QCD infrared scales, and cannot be calculated in perturbation theory.Both types of FFs have been the subject of intensive research [12,.Although the non-local type is considerably more complicated, they can be written in terms of local FFs at leading power in an operator product expansion (OPE) [36,37,39,43].Hence, understanding the local FFs precisely seems a more urgent matter.
The FFs F (T ),λ are functions of the momentum transfer q 2 and different determinations of the FFs apply at different values of q 2 .For instance, determinations based on Light-Cone Sum Rules (LCSRs) apply at low values of q 2 , while those based on Lattice QCD (LQCD) calculations have (until recently) only applied to the high end of the q 2 spectrum.Combining all the known information on the FFs in order to be able to predict observables in any given region of q 2 thus requires to understand the q 2 dependence.Fortunately, this q 2 dependence is simpler than the normalisation itself, and rigorous parametrizations of the q 2 dependence of the FFs can be derived from the fundamental principles of analyticity and unitarity [45].
These parametrizations are based on Taylor expansions of analytic functions, and while rigorous, in real applications they must be truncated in order to feature a finite number of parameters to be fitted.This truncation introduces a systematic uncertainty that is difficult to assess.Fortunately, one can derive absolute bounds on some combinations of FFs integrated over a kinematic region, by relating the exclusive to the inclusive rates and making use of a dispersion relation.These bounds are called "dispersive" or "unitarity" bounds, and effectively constrain the truncation error in specific parametrizations.The dispersive bounds have been used within the so-called BGL and BCL parametrizations [46][47][48], and have been proven extremely useful in B → D ( * ) [49] and Λ b → Λ ( * ) semileptonic transitions [50,51].
The purpose of this paper is to revisit the parametrizations and the dispersive bounds for B → P and B → V FFs.We shall demonstrate that a stronger version of the dispersion bounds of Ref. [47] exists, and that the parametrization of Ref. [44] is advantageous when below-threshold branch cuts appear.We will also apply these improvements to carry out a dispersive analysis of B → K ( * ) and Bs → ϕ FFs.
We begin in Sec. 2 with a detailed discussion of the theoretical framework, first defining the FFs and the parametrizations, and then discussing the dispersive bounds and our improved version.In Sec. 3 we present a complete numerical analysis where we fit the to the parameters of our FF parametrization including the improved dispersive bounds.Our conclusions are presented in Sec. 4. Supplementary material and a discussion of the polynomials used in the parametrization are provided in App.A and App.B.

Theoretical framework
In this section we provide the theoretical framework that we use to perform our analysis of Sec. 3. In Sec.2.1, we define the B → M form factors (FFs) starting from the corresponding matrix elements.Then, we briefly review the FF parametrizations written in terms of the conformal variable usually denoted by z, which have been commonly used in the literature.In Sec.2.2, we introduce the dispersive bounds, which constrain the parameters of certain z parametrizations.We derive new improved dispersive bounds that supersede the common ones of Ref. [47] by removing some of the spurious correlations among the FFs that arise in dispersive analyses and by giving stronger constraints for some FFs.In Sec.2.3, we give an accurate description of the FF parametrization that we adopt, which is based on the parametrization of Ref. [44] and allows for a straightforward implementation of the dispersive bounds.

Form factor definitions and z parametrizations
Hadronic FFs are scalar functions of the momentum transfer squared q 2 = (p − k) 2 that arise in the Lorentz decomposition of exclusive hadronic matrix elements.We define the B → P (seudoscalar) or B → V (ector) meson FFs following the same notation as used in Ref. [31], which relies on a commonly-used Lorentz decomposition.In this work we focus on the B → K, B → K * , and Bs → ϕ transitions -which we collectively denote as B → Malthough both the definitions given below and our results of Sec.2.2 and Sec.2.3 can be used for any B → P , V transition.Note that we do not make any attempt here to account for the instability of the vector meson [32].
The three local B → P FFs are defined by 2) The seven local B → V FFs are defined by where η is the polarisation four-vector of the vector meson, and we abbreviate (2.7) The local currents J µ Γ used above read (2.8) Throughout this work, we keep the q 2 dependence of the FFs implicit, unless necessary.In order to diagonalize the dispersive bounds (see Sec. 2.2), we introduce the helicity FFs where is the Källén function.The absence of dynamic singularities at q 2 = 0 and the antisymmetry of σ µν under the exchange µ ↔ ν imply the following relations Due to the symmetries of the helicity amplitudes at the kinematical endpoint q 2 = s − ≡ (M B − M V ) 2 , there are two additional constraints on the FFs [52]: (2.12) (2.13) The FFs are sensitive to infrared QCD scales and thus cannot be directly calculated in perturbation theory.Lattice QCD (LQCD) provides a means to perform such non-perturbative calculations.However, most of the presently available LQCD calculations for the FFs provide results only in the high-q 2 region of the semileptonic phase-space.Therefore, one needs to either extrapolate the LQCD results to the low-q 2 region or to combine them with results obtained at low-q 2 using other methods.In both cases a suitable parameterization of the FFs is needed.In the absence of LQCD data at all or in parts of the phase space, an alternative QCD-based method is the framework of light-cone sum rules (LCSRs).The LCSRs provide estimates of the FFs at low q 2 , which can be used to anchor their parametrizations.Their use incurs an irreducible systematic uncertainty due to the modelling of the corresponding distribution amplitudes and the use of semi-global quark-hadron duality [53].
The most frequently-used parametrizations in the literature describing the B → M FFs are the Boyd-Grinstein-Lebed (BGL) [46,47], the Bourrely-Caprini-Lellouch (BCL) [48], and the Illustration of the relevant regions in momentum space of the FFs both in the q 2 and the z plane.The semileptonic (SL) region extends from q 2 = 0 to , the branch cut starts at q 2 = s Γ , while the region for q 2 ≥ s + is the one relevant for the dispersive bounds (see Sec. 2.2).
Bharucha-Straub-Zwicky (BSZ) [30] parametrizations.Even though the BSZ parametrization was proposed more than thirty years later than the BGL one and is inspired by it, the former has a much simpler form than the latter since it is not devised to implement the dispersive bounds.For this reason, we first review the BSZ parametrization in this section in order to introduce the relevant notation, while we discuss the BGL parametrization in Sec.2.2 and Sec.2.3.To understand how these parametrizations have been introduced, one first needs to understand the analytical structure of the B → M FFs.Due to on-shell particle production, a FF arising from a specific current J µ Γ contains a number of simple poles from QCD-stable one-particle states, and develops a branch cut on the positive real axis for q 2 ≥ s Γ , where s Γ is the first multi-particle threshold.The dispersively bounded parametrizations BGL and BCL are only valid in the particular case where s Γ coincides with s + ≡ (M B + M M ) 2 .We note that s Γ = s + holds for processes like B → π or D → π, but it does not hold for the processes considered in this work: B → K, B → K * , and Bs → ϕ.In fact, the vector and tensor B → M FFs have s V = s T = (M Bs + M π 0 ) 2 and the axial and axial-tensor B → M FFs have s A = s AT = (M Bs + 2M π 0 ) 2 .In order to deal with the case s Γ < s + , we use the parametrization first proposed in Ref. [44], which we review in Sec.2.3.A recent alternative way to achieve the application of the dispersive bounds when s Γ < s + is discussed in Ref. [54] for B s → K transitions.The improvements to the dispersive bounds presented in this work can be straightforwardly applied to the approach of Ref. [54].
Following these considerations, we use the conformal mapping where s 0 is a free parameter that can be chosen in the interval (−∞, s Γ ).It is convenient to set , which minimises |z| in the semileptonic phase space.The transformation Eq. (2.14) maps the complex q 2 plane onto the complex z plane.In particular, the first Riemann sheet of the complex q 2 plane is mapped onto the open unit disk centred around the origin z = 0, and the second Riemann sheet is mapped onto the complement of the closed unit disk.The branch cut connecting both sheets is mapped onto the unit circle.
We provide an illustration in Fig. 1.The FFs are analytic functions of z on the open unit disk, except for a finite number of simple poles due to the aforementioned on-shell one-particle states, and which can be easily accounted for (see Sec. 2.3).Once the poles are removed, the resulting FFs can be expanded in a Taylor series that converges for |z| < 1.This is the essence of all z parameterizations, including BCL, BGL and BSZ.The latter parametrization reads where is a B → M FF and M F is the mass of the lowest-lying one-particle state with same spin, parity, and flavour quantum numbers as F B→M Γ,λ .These masses can be found in Tab. 3 of Ref. [30].
The simple form of Eq. (2.15) makes this parameterization convenient and its convergence in the semileptonic region is ensured by the analytic properties of the FFs.In any phenomenological application, the parametrization (2.15) would be truncated at some order n = N .However, it is not readily possible to estimate the truncation error, that is to estimate the uncertainty due to the truncated terms with coefficients c n for n > N .To estimate the truncation error, one needs to use additional constraints.These constraints can be obtained using unitarity and analyticity, and are usually called dispersive (or unitarity) bounds.In Sec.2.2 and Sec.2.3 we provide the derivation of these constraints, discussing both their original form of Refs.[46,47], where they have been applied to B decays for the first time.Then, we introduce our suggested modifications that strengthen these bounds and remove spurious correlations between FF parameters.

Improved dispersive bounds
The application of dispersive bounds is commonplace for charged-current heavy-to-heavy transitions, such as B → D ( * ) [47,[54][55][56].Their application to the local FFs in B → K( * ) transitions has been discussed for the first time in Ref. [57].In this work, we improve upon this analysis in two ways.First, by analysing simultaneously the B → K, B → K * , and Bs → ϕ transitions, and second, by strengthening the BGL-style bounds [46,47] that were used in Ref. [57].To derive dispersive bounds for any of the B → M FFs, one commonly uses the two-point correlation functions where the currents J µ Γ have been defined in Eq. (2.8).It is convenient to decompose the tensor-valued correlator Π µν Γ (q) in terms of scalar-valued functions Π One possibility is to use structures with definite total angular momentum J: This decomposition has been frequently used in the literature [47,54,57].Its connection to the total angular momentum J provides insights into the FFs' properties.However, it is not the most general decomposition possible.
In this article, we suggest to decompose the rank-two tensor in terms of a full set of polarisation vectors ϵ(λ) ≡ ϵ(λ; q) for λ = t, 0, ⊥, ∥.These polarisation vectors arise in the description of a virtual vector boson with four-momentum q and form a basis of Lorentz vectors.Our proposal then reads One readily finds that the quantities Π vanish for λ ̸ = λ ′ .This is a consequence of conservation of angular momentum of the (hypothetical) virtual vector boson described by the vectors ϵ.Thus, we introduce the following short-hand notation for the diagonal terms: As we show below, our proposed decomposition brings along a tangible benefit to our analysis.
In fact, it splits the contributions from B → P and B → V FFs into individual helicity-specific bounds and hence it removes spurious correlations between FF parameters for a single exclusive process.At the same time, it retains the correlations between helicity FFs across different exclusive processes, e.g., between A B→K * 1 and A Bs→ϕ 1 .To derive these improved dispersive bounds, we work in a frame where the momentum q is aligned with the z axis.We define the polarisation vectors as Using that σ µα is an antisymmetric tensor, one finds that Π (t) Γ = 0 for Γ = T, AT .The discussion so far is only helpful if a connection can be made between the scalar functions Π (λ) Γ as defined in Eq. (2.19) and the hadronic FFs.Key to this connection is that the former fulfil a subtracted dispersion relation: (2.23) Bs (5.367, 0.2307( 13)) Γ considered in this work with the corresponding OPE results (taken from Ref. [57]), affected FFs, and sub-threshold poles.The reference value for the b-quark mass is m b = 4.2 GeV.The parameters of the poles are taken from Refs.[58][59][60][61][62].
Here Q 2 is the subtraction point and n is the number of subtractions, which is chosen a posteriori such that the integral on the r.h.s. of Eq. (2.23) converges.On the one hand, the discontinuities of Π (λ) Γ are known from a computation within a local OPE, and the results are valid for Q 2 ≪ (m b + m s ) 2 .These results, denoted here as Π (λ) Γ OPE , manifestly fulfil Eq. (2.22).We use the results provided in Ref. [57], which include perturbative corrections up to next-to-leading order in α s for the dimension-0 term and power corrections due to operators up to dimension 5.For convenience, we quote the numerical results of Ref. [57] in Tab. 1, alongside further information.
On the other hand, unitarity of the S-matrix implies that the imaginary part of the functions Π λ Γ can be expressed as an infinite sum of exclusive hadronic matrix elements of the currents J Γ (see, e.g., Refs.[47,54,57]): Here we sum over all hadronic states H with flavour quantum numbers S = −B = 1 and we integrate over the corresponding phase space.For every state H, we can isolate a positive semidefinite exclusive contribution to χ Γ , which we denote as χ (λ) Γ H .These terms are not expected to fulfil Eq. (2.22), since they can -in general -resolve the polarisation of the incoming/outgoing virtual vector boson, but they still obey the relation (2.25) Nevertheless, the one-particle contributions fulfil Eq. (2.22) since they cannot resolve the polarisation of the incoming/outgoing virtual vector boson: (2.26) For one-particle bound states H = B( * ) s , the contributions read (2.27) The decay constants are defined as , where η is the polarisation four-vector of the B * s meson.Some comments are in order: we do not include all known one-particle bs states.Rather, we consider only bound states, i.e., only states with masses below the current-specific threshold s Γ .This is necessary, since bs states above their respective threshold(s) are resonances within the two-particle spectrum.As such, they emerge on the FFs' second (or higher) Riemann sheets rather than the first sheets.As a consequence, the resonances' contributions are accounted for by contributions of the various FFs, both two-particle FFs and those for higher multiplicities.In this way, we avoid double counting of the resonances' contributions.
If H is a two-particle state, the H-to-vacuum matrix elements can be related through crossing symmetry to hadronic FFs.For this work, we restrict our analysis to the states H = BM ≡ BK, BK * , Bs ϕ.We obtain the exclusive contributions χ (λ) Γ H following the derivation of Ref. [47], and for λ = t we recover the results in that reference.For H = BK, our results read (2.30) For H = BK * , they read (2.31) (2.32) (2.34) (2.35) where λ kin was defined below (2.10).The contributions arising from H = Bs ϕ can be obtained from those arising from H = BK * with obvious replacements.We assume that isospin symmetry holds for the transition FFs discussed here.The isospin factor η B→M is equal to 2 for B → K and B → K * transitions, while is equal to 1 for Bs → ϕ.
The dispersive bound is now obtained by equating the OPE result and the unitarity relation in Eq. (2.24): Analogously, equating the OPE result and the unitarity relation for χ Γ , one obtains the common BGL dispersive bound: Here, the ellipsis denotes the contribution of further states with the right quantum numbers, such as Λ b Λ or BKππ.Each individual hadronic contribution to χ Γ is a manifestly positive quantity since it can be expressed as the modulus squared of a hadronic matrix element.Hence, Eq. (2.38) can be converted from an equality to an inequality simply by dropping unknown but positive exclusive contributions.Here, we drop the terms represented by the ellipsis.Using our knowledge of χ (λ) Γ OPE , Eq. (2.38) now provides an upper bound on some positive definite integral involving at least one of B → M FF.This inequality is commonly called the dispersive (or unitarity) bound.The dispersive bound as expressed in Eq. (2.38) has two clear advantages with respect to the common one as expressed in Eq. (2.39).First, each FF of any of the transitions considered in this work obeys an individual dispersive bound.This becomes apparent from the results shown in Eqs.(2.28)-(2.37),and it is one of the major results of this work.Obtaining independent bounds for each hadronic FF is beneficial also from a statistical perspective, since it removes spurious correlations between the parameters of any two FFs that happen to contribute to a common bound of the form of Eq. (2.39).Second, although the OPE results cannot resolve between the polarizations λ = 0, ∥, ⊥, we find that the individual exclusive hadronic representations do.As a consequence, the dispersive bounds are stronger when they are formulated for the polarizations λ = 0, ∥, ⊥. 2 To illustrate this statement, consider two FFs (e.g., A B→K * 1 and A B→K *

12
) that contribute jointly to a common correlator in the common decomposition of Eq. (2.18): 3χ Here, the contribution with λ = ⊥ emerges only once baryon-to-baryon FFs are taken into account [50,51].In the common setup of the bounds, the λ = ∥ term can overfluctuate, i.e. exceed our proposed improved bound, as long as it either • does not exceed the size allocated to the λ = ⊥ term by the OPE; or • is compensated by an underfluctuation of the λ = 0 term.
The latter can readily emerge, e.g., from anticorrelations amongst the FF parameters.In our proposed setup, neither of these methods of evasion are successful, since each of the FFs is individually bounded.
We conclude this section by emphasising that the improved dispersive bounds are not the spurious effect of additional assumptions.Instead, they arise naturally when considering the origin of the dispersive bounds from physical observables.The correlation functions Π µν Γ introduced for the purpose of the bounds are related, by a dispersion relation and the optical theorem, to certain kinematic moments of the inclusive cross section σ(e + e − → hadrons), where the hadronic final state must feature S = −B = 1.The improved bounds proposed here arise when considering kinematic moments of the cross section for polarized e + e − beams.Clearly, our improved bounds can be applied in the same fashion to other transitions like Bq → D ( * ) q .

Dispersively Bounded Parametrization
The dispersive bound of Eq. (2.38) is not written in a form that allows us to use it easily in phenomenological applications.To express it in a more convenient form, the dispersively bounded z-series parametrization has been developed [44,47,54,57,63].In the following, we revisit the parametrization of Ref [44] and propose minor modifications that improve its usefulness and accuracy.To this end, we first briefly review the analytic structure of the B → M FFs.As anticipated in Sec.2.1, the FFs arising from a current J µ Γ develop a branch cut at the first multi-particle threshold, that is q 2 ≥ s Γ .For vector and tensor FFs, this threshold corresponds to the Bs π 0 pair production: s V = s T = (M Bs + M π 0 ) 2 .For axial and axial-tensor FFs, simultaneous conservation of angular momentum and parity requires at least two pions in the final state, leading to s A = s AT = (M Bs + 2M π 0 ) 2 .
By applying the mapping Eq. (2.14) to the dispersive representation of the two-particle contributions in Eqs.(2.28)-(2.37),we can rewrite Eq. (2.38) as where the sum runs over the transitions considered in this work and we abbreviate where P F is the FF's Blaschke factor, and ϕ F is the FF's outer function.The Blaschke factors P F are needed to cancel any simple poles in within the open unit disk in the z plane.These poles are due to the bound states R F discussed in Sec.2.2, listed in Tab. 1 alongside the corresponding FFs.We define the Blaschke factors as The expressions (2.28)-(2.37),once normalised by their respective χ (λ) Γ OPE , fix the norms of the outer functions ϕ F on the unit circle Here the parameters {N F , p, n, m} are FF-specific and listed in Tab. 2. Given that the modulus squared of an outer function is fixed only for |z| = 1, there is some leeway in how to continue the outer functions into the open unit disk.One requirement is that they be free of (kinematical) singularities in the open unit disk.Another one is to choose their overall phase to be realvalued on the real z axis, reflecting the fact that their corresponding FFs only develop a branch cut at s Γ and their global phase can be chosen to be zero.Setting Q 2 = 0, one obtains (see e.g.Ref. [47,57,63]) Table 2: Parameters of the outer functions of the B → P and B → V FFs, with The BGL parametrization only works in the special case where where the equality follows from the fact that the z monomials are orthogonal on the unit circle.Since here we deal with the case where s Γ < s + , we follow Ref. [44] and instead of using Eq.(2.47) expand FB→M Γ,λ (z) in a series of polynomials p n (z) that are orthonormal on a symmetric arc on the unit circle {e and hence Eq. (2.43) can be written as Here it suffices to say that p F n is a polynomial of degree n.Some mathematical properties of these objects as well as a discussion on the convergence of the series a F n are relegated to App.B.
The advantage of using the expansion Eq. (2.49) is that the dispersive bound Eq. (2.41)including the one-particle contribution -can now be written in the simple form where the first sum runs over all the FFs F B→M Γ,λ .To summarise, using unitarity and analyticity one obtains the dispersive bound Eq. (2.38).The derivation of this bound is similar to the original derivation of Ref. [47], with the crucial difference that in our case FFs with different indices λ obey different dispersive bounds.Then, applying a conformal mapping, we recast the bound Eq. (2.38) in the more convenient form of Eq. (2.52).In doing this, we took into account the analytic properties of the FFs F B→M Γ,λ and in particular the branch cuts that appear for q 2 < (M B + M M ) 2 due to on-shell Bs π 0 (π 0 ) states.These branch cuts have always been neglected for B-meson decays, which introduces a hard-to-quantify systematic uncertainty.Here, we take them into account following the procedure proposed in Ref. [44] for the non-local FFs in B → M ℓ + ℓ − decays.The dispersive bound Eq. (2.52) is an extremely powerful and model-independent constraint on the coefficients of the expansion Eq. (2.49).In practice, it allows us to control the systematic uncertainties due to the truncation of same expansion.

Analysis setup
We now study the available theoretical data on the FFs within a Bayesian analysis.Its central element, the posterior PDF P (⃗ a | D), is a function of the data D and our parameters ⃗ a ≡ (a F 0 , . . ., a F N , a F ′ 0 , . . ., a F ′ N , . . .).We use the same value of N as the order of truncation of the series Eq. (2.51) for each FF.The posterior definition involves our choice of the prior PDF P 0 (⃗ a), the likelihood P (D | ⃗ a), and the data-dependent normalisation (also called the evidence) Z(D).
We use different choices of fit model and prior.We label all our fits models by N ∈ {2, 3, 4}, the truncation order applying to all FFs.The choice of the prior PDF is motivated by the weakest-possible form of the dispersive bound, i.e., |a k | ≤ 1 for all k.This prior factorises to where U −1,+1 is the uniform PDF on the interval [−1, +1].Since we have 17 independent FFs and 9 end-point relations, the number of independent parameters is K = 17(N + 1) − 9.
We face two types of likelihoods in this analysis: linear constraints on the FF parameters provided by LQCD and LCSR determinations; and quadratic constraints on the FF parameters provided by the dispersive bounds Eq. (2.52).The linear multivariate Gaussian likelihoods arising from a single hadronic transition are labelled appropriately, i.e., either B → K, B → K * , or B s → ϕ.The data underlying each likelihood are summarised as follows: B → K We use three LQCD analyses of the three B → K FFs.The underlying LQCD analyses have been performed by the Fermilab Lattice and MILC collaborations (FNAL + MILC) [64] and by the HPQCD collaboration [33,65].Although an average of the FNAL+MILC and the 2013 HPQCD results is available from the Flavour Lattice Averaging Group (FLAG) [66], we choose to use the individual results rather than their average here.This enables us to diagnose which of the individual analyses are in mutual agreement.As indicated in the FLAG report, the FNAL+MILC and the 2013 HPQCD analyses rely -at least partially -on the same underlying gauge ensembles.However, their overall uncertainties are dominated by systematic sources, rather than the purely statistical uncertainty.Hence, we combine the results of both analyses under the assumption that they are uncorrelated, just as has been done to obtain the FLAG average [66].The third analysis, HPQCD 2022 [33], uses N = 2 + 1 + 1 gauge ensembles that are statistically independent of the previous analyses.The results of the three analyses are provided as multivariate Gaussian distributions in the parameters of a BCL expansion [48].We use the parametrized results of the various LQCD analyses to produce synthetic data points across the three FFs, making sure that the number of data points matches the number of degrees of freedom in the parametrizations.For the FNAL+MILC and HPQCD 2013 analyses, we generate points close to the q 2 values of the lattice ensembles, q 2 = {17, 20, 23} GeV 2 .The HPQCD 2022 results provides for the first time access to an exclusive B → K FF at q 2 = 0. Consequently, we spread the points over the entire physical range, q 2 = {0, 12, 22.9} GeV 2 .
The B → K FFs are also available from LCSR analyses with kaon distribution amplitudes.The most recent results obtained in this framework are available from Ref. [67].
Compared to the LQCD results, these estimates feature large parametric uncertainties and hardly-quantifiable systematic uncertainties.As a consequence, we do not use these estimates to obtain any of our numerical results and only use them for illustrative purposes and cross checks.
B → K * The full set of FFs is available from a LQCD analysis and a subsequent addendum [68,69].The authors of that analysis have provided us with twelve synthetic data points per FF, corresponding to the extrapolation of their LQCD ensembles.Correlations between the FFs arising from (axial)vector and (axial)/tensor currents are not presently available and can therefore not be accounted for in our analysis.
Beside from this LQCD analysis, numerical results for the full set of FFs are also available from LCSR with B-meson distribution amplitudes [31].The authors of this calculation provide five synthetic data points at q 2 = {−15, −10, −5, 0, 5} GeV 2 for each FF, and the full correlation matrix across FFs and data points are available.Following the argument in our previous paper [20], we do not use LCSR with light-meson distribution amplitude for FFs with vector final states.
B s → ϕ The inputs are similar to the B → K * ones.Refs.[68,69] provide the mean to create two synthetic data points per FF, reflecting the smaller number of LQCD ensembles for this transition.Ref. [44] provides LCSR estimates of these FFs using a previous calculation published in Ref. [31].Five synthetic data points including their correlations are available at q 2 = {−15, −10, −5, 0, 5} GeV 2 for each FF.
The dispersive bounds discussed in Sec.2.2 are labelled by their respective spin structure Γ and polarisation state λ.The penalty for violating the dispersive bounds is implemented following Ref.[55]: Here σ reflects the relative uncertainty on the computation of χ Γ OPE , which is estimated to be 10% [57], and r Γ,λ is defined as The quantity r Γ,λ is the statistical estimator for the saturation of the bound with polarisation λ for the current J Γ defined in Eq. (2.52).
Imposing the dispersive bounds as expressed in Eq. (2.39) makes a combined analysis of the FFs worthwhile but also challenging.It is worthwhile since in the presence of the dispersive bounds some of the parameterization uncertainties are shared among the fitted transitions.It is challenging due the large number of fit parameters K, e.g., K = 59 in our nominal analysis with N = 3.Our proposed decomposition of the correlator Eq. (2.19) renders the dispersive bounds stronger (see Sec. 2.2).It also has an additional advantage at the level at which we are working: in the absence of baryon-to-baryon FFs, our overall likelihood factorises, thereby decoupling pseudoscalar-to-pseudoscalar from pseudoscalar-to-vector transitions.This allows to perform a dispersively bounded fit for B → K transitions separately from a joint fit to B → K * and Bs → ϕ transitions.Since a joint fit to the latter two involves only correlations due to the dispersive bound, we can further simplify the analysis.Separating B → K * and Bs → ϕ transitions into individual fits, we can achieve the results of the joint fit by re-weighting their individual parameter samples with their dispersive bounds' likelihood.Taking B → K * as an example, the current-specific sample weight is computed as where r 1pt Γ,λ and r B→K * Γ,λ are the saturations of the corresponding bound due to the one particle and the B → K * process respectively, and p Bs→ϕ Γ,λ is the PDF of the saturation due to the Bs → ϕ process.The global event weight is obtained by multiplying these statistically independent weights, ) .
This procedure can be extended to account for other transitions, e.g., baryonic transitions.
It permits us to split the production of Monte Carlo samples for our joint analysis into three individual samplings, one per transition, which facilitates the overall analysis significantly.
We carry out the fitting and sampling parts of our analysis using the open-source EOS software [70] version v1.0.7 [71].For each value of the truncation order N , we determine the best-fit point and overall goodness-of-fit diagnostics.Our a-priori p-value threshold for an acceptable fit is 3%.We draw all posterior samples using the nested sampling algorithm [72] as implemented in the dynesty software [73,74].The results are investigated in form of posterior-predictive distributions for a variety of pseudo observables, including all of the FFs and the relative saturation of the various dispersive bounds.

Numerical results
Our analysis is repeated for three values of the truncation order N ∈ {2, 3, 4}, as anticipated in the previous section.In agreement with the literature [30,31,33,57,64,65,68,69], we already find a good fit to all FFs at N = 2.The substantial uncertainties attached to both LQCD and LCSR estimates of the FFs lead to local p-values close to 1 for both B → K * and Bs → ϕ transitions; we find all local p-values to be in excess of 77%. 3 We summarise the goodnessof-fit diagnostics for N = 2 and N = 3 in Tab. 3 and we consider the N = 3 analysis to yield our nominal results.We find that higher truncation orders do not have an adverse impact on the overall goodness of fit for nearly all individual likelihoods.However, in the case of Bs → ϕ the χ 2 value is so low that increasing the truncation order beyond N = 3 can only lower the local p-value.
We illustrate the results of our analyses at the hand of a representative subset of FFs in Fig. 2. Plots for the full set of FFs are contained in the supplementary material [75], alongside the EOS analysis file that has been used to obtain all numerical results shown in this section.Using Eq. (2.52), we further compute the saturations of the full set of dispersive bounds and present them in Tab. 4 and Fig. 3.We use this information later on to determine if the parametric uncertainties account for the truncation error.
We summarise our findings arising from the individual fits as follows: Fit with N = 2 We confirm the observation of Ref. [57] that for this value of N most bounds are not yet saturated and only play a marginal role in the fit.We can quantify this finding by determining the saturation at 68% cumulative probability based on the samples' posterior-predictive distributions.We find that only the bounds with (Γ, λ) = (V, ⊥),  (A, ∥), and (T, ⊥) exhibit a saturation of more than 50% at 68% cumulative probability, while all other bounds are saturated at less than 40%.The bounds for (Γ, λ) = (V, 0) and (V, t), which receive two-particle contributions only from B → K FFs, exhibit a very low saturation below 10%.

Goodness of fit
Fit with N = 3 As shown in Tab. 3, increasing the truncation order leads to a decrease of both χ 2 and χ 2 /d.o.f.across all the transitions.Increasing the truncation order mostly impacts the FF values and uncertainties at low and at negative q 2 , which is expected due to the dominance of the LQCD information at high q 2 values.Since the FF results in the region of low and negative q 2 are crucial for the estimation of non-local effects [20], this behaviour raises concerns toward the uncertainty estimation for the non-local effects.
We observe that the saturations of the posterior predictive samples now peak at larger values.We quantify this change by comparing the saturations at 68% cumulative probability for both N = 2 and N = 3 in Tab. 4.Although the bounds for (Γ, λ) = (V, 0) and (V, t) increase slightly in their saturation, we still find them to be saturated well below 20% at 68% cumulative probability.All other bounds now reach substantial saturations of 48% or larger.
Fit with N ≥ 4 For completeness, we also perform the fit with N = 4.We find that for all transitions the minimal χ 2 reaches a plateau, i.e., it does not significantly reduce as we increase the truncation order.As indicated above, this inevitably leads to a lowering of its local p-value as N increases.All bounds now exhibit a saturation of at least 50% at 68% cumulative probability.
We find that the average saturation of all bounds increases with increasing truncation order N .This is expected, since the parametric uncertainties on the posterior predictions of the FFs start to accurately include the extrapolation uncertainties.We find that this is achieved for N ≥ 3. Our interpretation is supported by the "comparison"-type plots contained in our supplementary material [75]; see App. A. Nevertheless, two of the bounds exhibit low saturations below 20% even for N = 3, which requires further discussion.The bounds at hand feature Γ = V and λ = 0, t and receive only contributions from B → K FFs in our analysis.
In a general dispersively-bounded analysis, there is an intrinsic ambiguity in how to interpret the fact that these bounds only receive comparatively small increases to their respective saturation as we increase the truncation order.On the one hand, the uncertainties attached to the theory estimates of the FFs might prevent our fit model from accurately capturing details on the shape of the FFs, which can lead to small saturations. 4In this case, the problem can be addressed by increasing the precision of the FFs' theory estimates.On the other hand, our fit model might accurately capture the saturation.Hence the observed low saturation values might be a physical feature of the responsible FFs.In this case, one would expect that including further relevant two-particle transitions and transitions with even higher multiplicities would drive the saturations to higher values.
In our analysis specifically, the data leads us to conclude that the second interpretation is the correct one.We hence argue that including baryonic transition FFs in the analysis could yield a significant impact of B → K FFs.As discussed above the bounds for (Γ, λ) = (V, 0) and (V, t) are only weakly saturated even for N = 3.First, we have explicitly checked that this is not leading to underestimated uncertainties in the FFs f B→K + and f B→K 0 , which are the only FFs that contribute to these two bounds in our setup.Our check is performed by evaluating these FFs for N = 4 at two phase space points: q 2 = −5 GeV 2 and q 2 = +25 GeV 2 , since in this scenario, these two bounds start to be saturated.The largest increase in the uncertainties with respect to the scenario N = 3 is obtained for f B→K + (−5 GeV 2 ) and is of the order of 0.2%.This very small increase in the uncertainties can also be seen in the "comparison"-type plots contained in our supplementary material [75].It is in line with our conclusion that N ≥ 4 is generally not needed to achieve an accurate estimate of the truncation error for all other FFs, and therefore further supports our choice of N = 3 as the nominal fit.Next, we compare our results with those of dispersively bounded analyses for baryonic transition FFs as presented in Refs.[50,51].We find that in both of these analyses the saturation of the Γ = V bound complements our results.Our interpretation of these findings is that the Γ = V , λ = 0, t bounds receive substantial contributions to their relative saturation from the baryonic transitions.This points to a tangible benefit of analysing mesonic and baryonic transition FFs simultaneously and might lead to even more precise predictions for the FF f B→K +,0 .

Conclusions
We propose novel dispersive bounds for local form factors (FFs) in B decays that improve upon the commonly-used BGL bounds of Ref. [47] Our improved bounds exploit the fact that the vacuum-to-vacuum correlator is not sensitive to the polarisation of the incoming/outgoing virtual vector boson, while the exclusive two-particle contribution to the correlator's hadronic representation is.Following this consideration, we split the BGL bound for FFs with total angular momentum J = 1 into three different parts, one for each physical polarisation.As a consequence, each FF of a single exclusive mesonic process now obeys a distinct dispersive bound.This has two advantages.First, it makes the dispersive bounds more constraining.This is possible since overfluctuations of one FF cannot be hidden by underfluctuations of another FF.Second, it avoids the introduction of redundant correlations.These arise in the common BGL-bound by merging contributions arising from different transitions or from different polarisations.
We apply our dispersive bounds by fitting a recently proposed parametrization to existing theoretical data on B → K( * ) and Bs → ϕ FFs.This data is available from lattice QCD analyses and light-cone sum rule calculations.The latter data is used to anchor the respective FFs at low q 2 whenever lattice QCD results for this region are unavailable.The proposed parametrization is needed in our analysis due to the appearance of below-threshold branch cuts in the FFs, generated by Bs π 0 and Bs π 0 π 0 states.These branch cuts have not been considered in previous analyses, giving rise to unwanted and hard-to-quantify systematic uncertainties.We obtain numerical results for the fitted parameters from a comprehensive Bayesian analysis of the thee transitions B → K( * ) and Bs → ϕ.We provide the N = 2 results in a machine readable form as supplementary material [75].Our results include the central values, standard deviations, and correlations of both the coefficients obtained in our framework and the coefficients of the BSZ parametrization.The parameters within both parametrizations follow multivariate Gaussian distributions to a good approximation, with relative perplexities in excess of 95%.At the current level of precision on the estimation of the FFs describing B → K( * ) and Bs → ϕ transitions, we find that at present a simplified series expansion à la BSZ provides an accurate estimation of the FF uncertainties at positive q 2 , although this may change once more and more precise lattice QCD results become available.However, at negative q 2 , the coupled effect of a larger truncation order and the impact of the dispersive bounds are already essential to correctly estimate the extrapolation uncertainties.This fact is particularly important for the determination of the non-local FFs for these transitions since at low and negative values of q 2 they are proportional to the local FFs at leading power in the light-cone OPE.Therefore, we conclude that our proposed approach is crucial for this determination.
We emphasise that present and future lattice QCD analyses must fulfil the dispersive constraints.Hence we encourage the implementation of our framework to perform a posteriori checks.Nevertheless, we caution that dispersive bounds should not be applied directly or exclusively in theoretical calculations; doing so would prohibit their use in simultaneous FF fits.
Our present analysis is both a benchmark of what is possible within the framework of dispersively bounded FF parametrizations and a step stone towards more precise predictions of the non-local FFs at negative q 2 .A key advantage of this framework is the possibility of including further transitions in the analysis, in particular baryonic ones.We expect this to lead to a reduction of the uncertainties of all the transitions considered within a global analysis.Thus, we regard the work presented here as an important step towards a truly global analysis of rare B-decay observables that simultaneously uses dispersive bounds to control all systematic uncertainties, thereby advancing the precision frontier in indirect BSM searches.
• At q 2 < 0, we start to see that BSZ parametrization does not faithfully capture the uncertainty envelope of the dispersively bounded analysis.This is not unexpected, given the fact that the truncation error becomes relevant for this type of extrapolation of the available LQCD data.It gives rise to concern for the theory prediction of non-local FFs, which are crucially reliant on this extrapolation.The biggest concern here is the FF f B→K T , which shows deviations of the order of 10%.

B Orthonormal Polynomials
Our parametrization of the FFs involves orthonormal polynomials on the unit circle.For the purpose of our analysis, orthonormality is defined with respect to the scalar product Eq.(2.50), i.e., orthonormality with respect to integration over the arc A α = {e iθ , −α < θ < α}.We could not find a closed analytic formula for these polynomials in the literature.Instead, we evaluate them recursively using the Szegő recurrence relation [76]: where φ n is the nth orthogonal polynomial.We adopt here the common notation used in the mathematical literature φ * ,n n (z) ≡ z n φn (1/z) [76].The coefficients ρ n , known as the Verblunsky coefficients, uniquely define the set of polynomials and only depend on the integration measure.In our approach the relevant integration measure is the constant Lebesgue measure on the arc A α .
Since this arc is symmetric under complex conjugation, the Verblunsky coefficients are real numbers [76].The orthonormal polynomials p n appearing in Eq. (2.49) are obtained using [76]

Calculation of Verblunsky coefficients
The Verblunsky coefficients can be calculated by applying the Gram-Schmidt orthonormalisation procedure to the basis of monomials z k .A faster evaluation can however be obtained recursively and does not require any analytic calculation.To illustrate our approach to this evaluation, we first define Φ n ≡ (φ n (z), φ * ,n n (z)) T , which fulfils the recurrence relation This equation can be used to evaluate the orthogonal polynomials at any point z, provided that the set of Verblunsky coefficient is known.To compute the later, we introduce two sets of integrals, The Verblunsky coefficient ρ m can then be obtained by recursively filling the triangular matrices I n,k , J n,k for k < n ≤ m and using ρ m = I m,1 /J m,0 .This method is implemented in the open-source EOS software [70] used for numerical aspects of this paper.

Asymptotic behaviours
The asymptotic behaviours of the orthonormal polynomials have been studied in details, see e.g.[77] and references therein.Using Eq. (1.7) of Ref. [78] and applying the transformation θ ′ = π − θ, we find that at large n, The values of the polynomial thus form a finite alternating series at z = 1 and an alternating and exponentially divergent series at z = −1, p n (−1) = O (cot α 4 ) n .

Convergence of the form factor expansion
Following the approach of Ref. [79], we evaluate Eq. (2.43) at the branch point z(s Γ ) = −1, As shown in Eq. (B.8), the series p F n (−1) diverges exponentially with n.This, however, is not an issue for the parametrization.The power series element of the parametrization is analytic for |z| < 1.Moreover, this element is also finite at z = −1.By Abel's theorem, the value at z = −1 can be expressed as

FB→M
is analytic on the open unit disk, it can be expanded in a Maclaurin series,

Figure 2 :
Figure 2: Selection of local B → M FFs.f B→K +

Figure 3 :
Figure3: Saturations of the dispersive bounds due to one-and two-particle contributions.The lines represent Gaussian-smoothed distributions of the saturations in our samples that we combined following the discussion in Sec.3.1.The shaded area comprise the 68% probability interval for each scenario.The distributions are scaled with arbitrary factors for readability.

ρ n ∼ (−1) n cos α 2 . (B. 7 )
This result can be used when diagonalising the matrix in Eq. (B.3) to obtain asymptotic values for any z ̸ = 0.In particular one finds that at large n,

10 )
Together with Eq. (B.8), this implies that the coefficients a n must fall off at least exponentially to compensate the divergence of the series {p n (−1)}.

Table 3 :
Goodness-of-fit values for the fit models N = 2 and N = 3.

Table 4 :
Relative saturations of the dispersive bounds due to the one-particle and two-particle contributions.Two-particle contributions are shown for the different truncation orders and the table shows the saturation at 68% cumulative probability.