Standard Model predictions for $B\to K\ell^+\ell^-$, $B\to K\ell_1^- \ell_2^+$ and $B\to K\nu\bar{\nu}$ using form factors from $N_f=2+1+1$ lattice QCD

We use HPQCD's recent lattice QCD determination of $B \to K$ scalar, vector and tensor form factors to determine Standard Model differential branching fractions for $B \to K \ell^+\ell^-$, $B\to K \ell_1^+\ell_2^-$ and $B \to K\nu \overline{\nu}$. These form factors are calculated across the full $q^2$ range of the decay and have smaller uncertainties than previous work, particularly at low $q^2$. For $B \to K \ell^+ \ell^-$ we find the Standard Model branching fraction in the $q^2$ region below the squared $J/\psi$ mass to exceed the LHCb results, with tensions as high as $4.2\sigma$ for $B^+\to K^+\mu^+\mu^-$. For the high $q^2$ region we see $2.7\sigma$ tensions. The tensions are much reduced by applying shifts to Wilson coefficients $C_9$ and $C_{10}$ in the effective weak Hamiltonian, moving them away from their Standard Model values consistent with those indicated by other $B$ phenomenology. We also update results for lepton-flavour ratios $R^{\mu}_e$ and $R^{\tau}_{\mu}$ and the `flat term', $F_H^{\ell}$ in the differential branching fraction for $\ell\in\{e,\mu,\tau\}$. Our results for the form-factor-dependent contributions needed for searches for lepton-flavour-violating decays $B\to K\ell^-_1\ell^+_2$ achieve uncertainties of 7%. We also compute the branching fraction $\mathcal{B}(B\to K\nu\bar{\nu})$ with an uncertainty below 10%, for comparison with future experimental results.


I. INTRODUCTION
The weak B → K transition involves a b → s flavour changing neutral current (FCNC) which proceeds at loop level in the Standard Model (SM), where contributions are suppressed by loop factors and associated elements of the Cabbibo-Kobayashi-Maskawa (CKM) matrix [1,2]. This makes FCNC processes promising places to look for effects from new physics.
Tests for new physics require the accurate calculation in the SM of quantities that can be compared to experimental results. Because the weak transition occurs between hadronic bound states of the b and s quarks, key inputs to the SM calculation are the hadronic matrix elements of the weak currents, parameterised by form factors. The form factors are functions of the 4momentum transfer, q, between B and K and depend on the strong interaction effects that bind the quarks inside the mesons. To calculate the form factors from first principles in quantum chromodynamics (QCD) requires the use of lattice QCD. Here we will update the phenomenology for B → K decays in the SM using improved form factors calculated using lattice QCD by the HPQCD collaboration [3] and compare to existing experimental results. We also provide tables of results in a variety of q 2 bins that will be useful for comparison to future experimental analyses.
The lattice QCD form factors that we use [3] improve on previous calculations in a number of ways. They are the first to be calculated on gluon field configurations * w.parrott.1@research.gla.ac.uk † chris.bouchard@glasgow.ac.uk ‡ christine.davies@glasgow.ac.uk § http://www.physics.gla.ac.uk/HPQCD (generated by the MILC collaboration [4]) that include u, d, s and c quarks in the sea (N f = 2 + 1+ 1). They use a fully relativistic discretisation of the QCD Dirac equation on the lattice [5] which allows accurate normalisation of the lattice weak currents to their continuum counterparts [6,7], rather than the O(α s ) perturbative normalisation used in earlier lattice QCD calculations [8,9]. The calculation must be done at a range of masses for the heavy quark, from that of the c quark upwards, and at a range of values of the lattice spacing. Very fine lattices must be used to get close to the b quark mass for the heavy quark. Fits to the results then enable form factors for B → K to be determined and smoothly connected (since only the heavy quark mass changes) to those for D → K. Our form factors for D → K agree well as a function of q 2 with those inferred from the experimental results for the tree-level weak decay D → K ν [10], also giving a better-than-1% determination of V cs . We do not expect to see new physics in this decay, so this is a strong test of the agreement of (lattice) QCD with experiment for a SM weak process. Our B → K form factors build on this. They are also the first to cover the full physical q 2 range of the B → K decay; this is possible because of the use of very fine lattices. This means in particular that our form factors have smaller uncertainties in the low q 2 region, allowing more accurate SM phenomenology here than has been possible before.
The most direct comparison between the SM and experiment is for the differential branching fraction as a function of q 2 . A complication is that of contamination by long-distance contributions from charmonium reso-nances in the intermediate q 2 region [24]. These are not easily included in the SM analysis (and will not be included here) and so the key comparisons to be made are at low q 2 , below the resonance region, or at high q 2 , above the resonance region. Tensions between the SM and experimental results have been seen in earlier work [8,9,[25][26][27][28]. Here we will provide an improved analysis with smaller theory uncertainties from the form factors, particularly in the low q 2 region. We will also update SM results for the lepton-universality violating ratios, R µ e and R τ µ (although form factor effects largely cancel in these ratios) and give results for the 'flat term', F H .
The lepton-number-violating decay B → K − 1 + 2 is of interest for new physics constraints and upper limits on this decay have been obtained by experiment [15,[29][30][31].
Here we provide the hadronic input needed for future analyses of this decay.
The B → Kνν decay promises to be a useful mode for study in future, because it is free of charmonium resonance contamination. There is a small amount of experimental data on this from Belle [32], Belle-II [33] and BaBar [34], with the promise of more to come [35]. Here we provide the SM predictions for the differential and total branching fractions for this decay mode with an improved level of accuracy compared to previous calculations.
The rest of the paper is organised as follows. Most of the paper is dedicated to the well-studied B → K − + decays. We begin with these decays in Section II, which outlines the calculation of SM observables in Section II A (including discussions of how we deal with resonances, and the effects from isospin, QED and lepton flavour), comparisons with experiment in Section II B and comparisons with previous theoretical results in Section II C. Section III looks at the lepton flavour number violating B → K − 1 + 2 decay, calculating several parameters relevant to this channel. In Section IV we calculate the branching fraction for the B → Kνν decay, comparing with experimental bounds and other theoretical work. We conclude in Section V. This version of our paper contains a correction of the scale used for α EW , compared to the published version. Instead of 1/α EW (M Z ) = 127.952 (9) [43], we take 1/α EW (4.2 GeV) = 132.32 (5), in line with the findings of [44]. This change affects our branching fractions for B → K + − and our results for B → K + 1 − 2 , by an overall multiplicative factor of α 2 EW (4.2 GeV)/α 2 EW (M Z ) ≈ 0.94. For most of our results, this leads to a change of approximately 1σ, which does not significantly impact our conclusions. All affected figures and tables have been updated accordingly, and will appear as an erratum. There is no effect on ratios of branching fractions such as R µ e and F H , nor on B → Kνν, as this uses α EW (M Z ) [45].
Additionally, we have removed a factor of η EW = 1.007 (2), which was incorrectly applied to G F in the original paper. This results in a modest reduction of 1.4% to all branching fractions.
The majority of theoretical and experimental work on B → K weak decays has focused on B → K + − . An example of a Feynman diagram for this process, sometimes called a penguin diagram, is shown in Figure 1. The analysis of such decays in the SM starts from the effective weak Hamiltonian constructed from 4-fermion operators multiplied by Wilson coefficients (C 1 -C 10 ), with the matrix elements of the 4-fermion operators expressed in terms of form factors [46]. This section discusses the calculation of SM observables related to these decays in Sec. II A, reviewing the inputs used in the evaluation of observables in Sec. II A 1 and focusing on the Wilson coefficients in Sec. II A 2. Section II A 3 then discusses how we handle resonance contributions in a way that allows comparison to experiment. Sections II A 4 and II A 5 then discuss the uncertainties to account for the fact that our form factors are determined in a pure QCD calculation on the lattice in the strong-isospin limit (m u = m d ). The discussion of uncertainties there is also relevant to our input form factors for Sections III and IV. Section II B gives our main results for SM observables and their comparison to experiment. Section II C closes the section with a comparison to earlier theoretical results.
A. Calculating SM observables The SM differential decay rate for B → K + − for lepton is constructed as follows, where we use the notation in [40,47].
where a and c are given by with F P,V,A are constructed from the scalar, vector and tensor form factors f 0 , f + and f T respectively, by The Wilson coefficient C eff,1 9 = C eff,0 9 + ∆C eff 9 + δC eff 9 includes nonfactorisable corrections in ∆C eff 9 , as well as O(α s ) and more heavily suppressed corrections in δC eff 9 . Similarly, C eff,1 7 = C eff,0 7 + δC eff 7 contains O(α s ) corrections in δC eff 7 . These corrections are discussed in detail in Appendix B.
C eff,0 where  [57] and running the scale to its own mass. mc and m b are the c and b quark pole masses obtained from the masses in the MS scheme at three loops (see Appendix A for details).

Inputs
Key inputs to our determination of the differential rate are our recently calculated scalar, vector and tensor form factors [3]. These are plotted as a function of q 2 in Fig. 2. The plot shows the ±1σ error bands, smaller than those for previous results, especially in the low q 2 region. See [3] for details on how these form factors were calculated. Note that we run f T (calculated in [3] at scale 4.8 GeV) to a scale µ b = 4.2 GeV, to match the scale used in the Wilson coefficients of Table I. This is achieved by following [7] and multiplying our form factor result f T (µ b = 4.8 GeV) by a factor of 1.00933 (16). The value for α EW (µ b = 4.2 GeV) is run from 1/α EW (4.8 GeV) = 132.14(4) [58].
Further external input parameters needed for Eq. (1) are provided in Table I. Except in the case of Equation (4), all expressions above (and in Appendix B) use the pole masses m c and m b given in Table I. These pole masses are derived from the masses given in the same table in the MS scheme, determined to high accuracy using lattice QCD calculations [48,57]. We convert to the pole mass using the standard perturbative matching factor at three-loop order [59,60]. For more details, see Appendix A. The perturbation series in this expression suffers from the presence of a renormalon in the pole mass. To account for renormalon-induced uncertainties [61] in the pole mass, we allow an error of 200 MeV which dwarfs any other uncertainty in the mass. We assess the impact of this uncertainty by shifting the central values of m c and m b down by 200 MeV. In the case of m b , we find almost no effect on the branching fractions in the important 1.1 ≤ q 2 /GeV 2 ≤ 6 and 15 ≤ q 2 /GeV 2 ≤ 22 regions.
In the case of m c , a maximal effect of 0.3 σ is seen.
To obtain the branching fraction B, we note that B = Γ τ B , where τ B is the lifetime of the B meson. The values for τ used are given in Table I.
The value of |V tb V * ts | that we use comes from the recent HPQCD lattice QCD calculation [55] of the mixing matrix elements for the B s meson, the first to be done with N f = 2 + 1 + 1 quarks in the sea. |V tb V * ts | is obtained by combining the lattice QCD result with the experimentally-determined B s oscillation rate [43]. The oscillation rate can be determined to high precision, so the uncertainty in |V tb V * ts | is dominated by that from the lattice QCD result, where the uncertainty is smaller than that for previous lattice QCD calculations. The value of |V ts | quoted in [55] used |V tb | = 0.999093 from [62]. The value of |V ts | that we use is about 1σ higher than values based on CKM unitarity that do not use the information from the B s oscillation rate, for example that in [63].

Wilson coefficients
The Wilson coefficients we use are from Table I of [56]. They are quoted at µ b = 4.2 GeV with uncertainties given for higher order corrections, as well as for a variation (doubling and halving) of the matching scale at the top quark mass, m t . These uncertainties do not allow for possible scale dependence in µ b , however. Whilst our differential decay rates should be scale independent, we test this over a reasonable variation of µ b and include an uncertainty to account for any variation. We look at the effect of changing µ b to 4.8 GeV on branching fractions (for both meson charges and ∈ {e, µ}) in the key ranges 1.1 GeV 2 ≤ q 2 ≤ 6 GeV 2 and 15 GeV 2 ≤ q 2 ≤ 22 GeV 2 . This involves the following changes: 1) Removing the running of f T (µ b ) described above, so we have f T (µ b = 4.8 GeV).
3) Shifting the central values of the Wilson coefficients to the 4.8 GeV values given Table I of [56]. 4) Setting 1/α EW (4.8 GeV) = 132.14(4) [58]. After making these changes, we observe a change of the same absolute size and sign in the branching fractions for both q 2 regions. This corresponds to a change of between 1.4% and 2.6% on the branching fraction, a maximum of 0.3σ. Though a modest effect, we account for this by including a 2.5% uncertainty on a and c across the full q 2 range to allow for reasonable variations in µ b .

Resonances and vetoed regions
From Fig. 1 it is easy to see how a resonance that couples to the photon can affect the differential rate for B → K + − , causing a spike at q 2 = M 2 for a resonance of mass M . These 'long-distance' resonance effects are not included in the 'short-distance' SM calculation of Eqs. (1)- (6). The resonances of concern are uu and dd resonances (the ρ and ω), which appear at very low q 2 , below 1 GeV 2 , the ss φ close to 1 GeV 2 and cc resonances, which affect the intermediate q 2 region. The cc resonances that dominate are the J/ψ and ψ(2S) and experimentalists typically discard their data in the vicinity of these resonances for a B → K + − analysis. Whilst the regions blocked out vary somewhat from experiment to experiment, we will use the choices 8.68 ≤ q 2 /GeV 2 ≤ 10.11 and 12.86 ≤ q 2 /GeV 2 ≤ 14.18, which are adopted by several experiments, and mark these as 'vetoed regions' (greyed out in our plots). In these two regions and between them, the resonance contributions are large [22] and there are strong violations of quark hadron duality [24]. The decay rate is dominated by B → KΨ → K + − [64] and is approximately two orders of magnitude larger than the short-distance rate of Eqs. (1)- (6).
To avoid these resonances in the comparison between the SM short-distance calculations and experiment it is common practice to restrict the q 2 range to two 'wellbehaved' regions, giving separate results for the integrals over 1.1 ≤ q 2 /GeV 2 ≤ 6 and 15 ≤ q 2 /GeV 2 ≤ 22 (or similar). The lower q 2 region here is above the light-quark resonances but below the cc resonances. The higher q 2 region is above the dominant cc resonances, but small contributions are seen in the experimental data from higher resonances in this region [22]. Following [9], we allow for an uncertainty of 2% on the differential rate in this high q 2 region to allow for these additional resonance contributions. We do this by applying a factor of 1.00(2) to a and c (Eq. (1)) above q 2 = 14.18 GeV 2 . We will adopt these two q 2 regions here and our key results will be a comparison with experiment for these two regions. This allows us to test for the presence of new physics in as clean a way as possible.
It is also informative, however, to quote the total branching fraction for the decay, integrating over the full q 2 range. This means that some approach must be adopted for the q 2 regions where the resonances sit so that a comparison can be made between the shortdistance SM branching fraction and experiment. In the very low q 2 region this is not an issue; the light-quark resonances contribute only at the 1% level [22]. The narrow φ resonance can be cut out of the experimental data, where the data is sufficiently accurate to warrant it [19]. The cc resonances in the intermediate q 2 region have a bigger impact over a larger q 2 range. The integration over this region must then be handled carefully in a way that is consistent with what is done in experimental analyses. To cover the cc resonance region, we linearly interpolate a and c across the whole of the two vetoed regions and the gap between them, i.e. for 8.68 ≤ q 2 /GeV 2 ≤ 14.18. This is only necessary for the case ∈ {e, µ}, as = τ sits mostly above this region. This approach is used in [9] and similar interpolations are performed in experimental papers (see [19], the interpolation of which will be compared with ours below). This means that the total branching fraction is then also comparable between theory and experiment. Of course, comparisons with other theoretical calculations that are based on Eqs. (1)- (6) can be done at any q 2 value.
Finally, it should be noted that the m c value we use ( Table I) is such that 4m 2 c = 11.34 GeV 2 . This sits between the two vetoed regions and therefore inside the region in which we perform the linear interpolation described. This means that the cusp arising at this point from the function h of Eq. (6) is not visible in our plots.

Isospin and lepton flavour
The B → K decays we wish to study are for B 0 → K 0 + − and B + → K + + − . The charged and neutral B meson cases differ only in the flavour of the light spectator quark, which can either be an up or a down. The masses of the B 0/+ and K 0/+ differ slightly as a result of this (along with QED effects). We take these masses from the Particle Data Tables [43] and use these in the equations of Sec. II A for each case; it has negligible effect on the differential rates. The difference in branching fractions between the charged and neutral cases is much larger, however, because of the difference in B + and B 0 lifetimes, which we take from [54] and give in Table I. We must also include an uncertainty for the fact that our form factors [3] were calculated in pure lattice QCD, with no QED effects included and with degenerate u and d quark masses (as were previous lattice QCD calculations of these form factors [40,41]). In order to work out how much our form factors would change if the u and d quark masses were different we follow the procedure adopted in HPQCD's D → K analysis [10].
We can repeat our analysis to fit the B → K form factors from [3] but changing the physical value of m l so that it corresponds to either m u or m d (this means taking m s /m l = 27.18 × 3/2 or 27.18 × 3/4 respectively). Note that this will give an overestimate of the effect because this change refers to all the light quarks in the simulation, whereas in fact we just want to change the valence quark in the B and K. At the same time as making this change to m l we also change the meson masses that appear in the fit function to their appropriate values corresponding to a u or d spectator quark. For details on the fit function for the form factors see [3].
Looking at the effects of these tests on the f 0 , f + and f T form factors across the q 2 range, we find a maximal effect of ≈ 0.5% [3]. To allow for the full effect as an uncertainty on our form factors for the analysis here, we apply an uncorrelated multiplicative factor of 1.000(5) to each of our form factors. We include this uncertainty in our form factors for all of the calculations discussed here, both in this Section on B → K + − and Sections III and IV.
The lepton masses, m e , m µ and m τ as appropriate, also appear in the equations to determine the differential rate (see Eqs. (2) and (3)). One of the notable changes between leptons is in the minimum q 2 value, since q 2 min = 4m 2 . In practice, in the case of the differential decay rate, our final result is almost identical for = e and = µ because their masses are both so small. We sometimes treat these cases together, simply using the = µ case where a lepton is not explicitly stated. Larger differences are seen for the = τ case and we show results for that case in Sec. II C where we compare to other theory calculations.

QED effects
The largest QED effects in the decay rate for B → K + − are expected to be from final-state interactions of the charged leptons. The photon radiation produced can cause several issues within the experimental analysis, particularly in the hadronic environment of the LHC. These are handled using the PHOTOS software [65] to correct the differential rate to a fully photon-inclusive rate that can be compared with theory. It is also used in the LHC context to obtain a value for q 2 . Theoretical tests [66,67] indicate that this is a robust procedure.
We should allow an uncertainty for QED effects in our comparison with experiment, however, and we will base it on the size of QED final-state radiative corrections [66,67]. This will then comfortably include the uncertainty from other QED effects missing from our calculation, such as that from the fact that the quarks inside the B and K mesons have electric charge and that there are additional structure-dependent radiative effects from the K mesons. We will quote the uncertainty for QED effects as a separate uncertainty in Tables where we quote values so that its impact is explicit. We will not include it in plots of the differential rate, nor in comparison with other theory results that do not include a QED uncertainty.
From [66,67], we take an overall multiplicative 5% (2%) uncertainty to allow for missing QED effects in the case of final state electrons (muons). This is a conservative move, but we shall see that this uncertainty has minimal effect on our results. We take the QED uncertainty for final state τ to be negligible because the τ is so heavy.
For the case of the lepton-flavour-universality-violating ratio R µ e , there has been a lot of analysis of QED effects [67][68][69]. There is agreement that QED effects are well accounted for using PHOTOS [65]. An additional 1% uncertainty for QED effects in the comparison between theory and experiment is suggested in [68] and we adopt that here when discussing R µ e in Section II B 3.

B. Comparison to experiment
In this section we will give our results for the SM shortdistance differential branching fractions (and integrals thereof) for B → K + − using Eqs. (1)-(6) and the improved B → K form factors of [3]. We will compare to experimental results and show where tensions exist, and how significant they are.
There are numerous experimental results in the literature from BaBar [11][12][13], Belle [14,15], CDF [16] and LHCb [17][18][19][20][21][22][23], for decays B 0 → K 0 + − and B + → K + + − where ∈ {e, µ}. There are few for = τ , which we will treat separately. Experiments are labelled in figures and Table II provides a correspondence between references and labels. Naturally, different analyses from one particular experiment are likely to be correlated, and should be viewed as such. One could take the most recent data from each of the experiments. However, owing in part to the large number of different choices made with regard to q 2 limits, meson charge and lepton flavour, we find it most informative to display all available results in general, even in the case where older measurements have been superseded by more recent analyses. This allows the reader to observe the progress of experimental results, in comparison with our own numbers. In the discussion of branching fraction B (for ∈ {e, µ}) that follows, we will group experimental data by B/K charge, as our results are sensitive to this. We will not group results by e or µ lepton flavour because our results (in the SM) are not sensitive to this (see Appendix C for tables including our results for different ). In the case of lepton flavour, experiments sometimes give e and µ values separately and sometimes both together. In this latter scenario we will call the lepton and take an average of e and µ. All these will be treated together for a given B/K charge. Similarly, results are given for a fixed meson charge, as well as for the average of the charged and uncharged case. We will use the following labelling conventions, which are adopted in many experimental papers. Where a charge is specified, it will be labelled, with the unlabelled case indicating the average of the charged and neutral cases. Similarly, e and µ leptons will be labelled explicitly, with indicating an average of the two. Where possible in plots comparing values with experiment, one set of charges or lepton flavours will be given in red, the other in blue, and the average in purple.

Differential branching fraction and integrals over low
and high q 2 regions avoiding cc resonances Figures 3, 4 and 5 show our calculation (as a blue band) of the differential branching fractions for the cases of positively charged mesons, neutral mesons, and the average of the two, respectively, across the full physical q 2 range. Note that the different charges give shapes which vary not just because of the overall factor of the lifetime τ B 0/+ , but also because of the charge dependence of the correc- FIG. 4. Differential branching fraction for B 0 → K 0 + − , with our result in blue, compared with experimental results [16,17,19]. All experimental results take = µ. Horizontal error bars indicate bin widths. LHCb '14A LHCb '14A FIG. 6. Comparison of branching fractions with recent experimental results [15,19,23] in low and high regions of q 2 away from the charmonium resonance region. Here we show the ratio of the experimental branching fraction to our results, compared to the black vertical line at the value 1. The error bars are 5σ long, with markers at 1, 3 and 5σ. Note that the σ here are for the ratio, so not the same as those calculated for the difference in Table III. On the right, labels indicate the colours of the q 2 bins in units of GeV 2 . No uncertainty from QED is included in this plot.
tions to C eff,0 9 detailed in Appendix B. The black vertical bands indicate regions vetoed because of the J/ψ and ψ(2S) resonances, not included in our calculation. As discussed in Section II A 3, our result is interpolated across these regions and the gap between them.
Experimental results for decays to electrons, muons or both (averaged) are displayed in each case as coloured points, with the results shown for each experimental q 2 bin. The horizontal error bars on the experimental results reflect the width of the bin. Some of the experimental results are for = e and some for = µ; our results are insensitive to the difference. The experiments ignore data taken in the black vetoed regions, but there are results in between these regions. However, we cannot make a reliable comparison between our short-distance SM results and the experimental results between the vetoed regions.
Figures 3, 4 and 5 show that our results are somewhat higher than experiment in most cases, particularly in the region 4 ≤ q 2 /GeV 2 ≤ 8.68. This is most clearly visible in Figure 3, where the tension between our result and the most precise data from LHCb is obvious.
To examine this tension in more detail, we integrate over two well-behaved q 2 regions, one above and one below the cc resonances, as discussed in Section II A 3. For these regions we can make a reliable comparison with experiment. We show the results in Table III; these constitute our main numerical results. In Table III, we compare our branching fractions with the most recent experimental results available for B → Ke + e − and B → Kµ + µ − . Note that our relative uncertainties are comparable to those from the experiments for most of the values. We have larger uncertainties than those for LHCb '14A for

Channel
Result q 2 /GeV 2 range B × 10 7 Tension with HPQCD '22  [15,19,23] in well behaved regions of q 2 . Note that the B + → K + e + e − result quoted here from LHCb '21 is obtained using the B + → K + µ + µ − result from LHCb '14A, combined with the ratio determined in LHCb '21. In the fifth column, the tension is given by mean(Experiment -HPQCD)/σ(Experiment -HPQCD). An additional uncertainty for QED of 5% for e and 2% for µ is included in brackets on each branching fraction (see Section II A 5). Additionally, the final column gives the tensions with and without QED uncertainty -those including this uncertainty are in brackets.
B + → K + µ + µ − but smaller uncertainties than those from Belle '19. Our uncertainties are dominated by those from the form factors, followed by those from the CKM elements |V tb V * ts |. We find our partial branching fractions to be significantly higher than the experimental values, typically with tensions exceeding 2.7σ and with LHCb'14A [19] in the low q 2 region with a tension reaching 4.2σ for B + → K + µ + µ − . We do not see any tension with the Belle'19 results [15], which are themselves in tension with LHCb '14A for B + → K + µ + µ − . Belle '19 have considerably larger uncertainties than those from LHCb and those from our results. The uncertainties we allow for QED make little difference to the tensions that we see.
The same results, plotted in the form of a ratio of the experimental branching fraction to our result, are shown in Figure 6, with error bars indicating 1, 3 and 5σ differences from a value for the ratio of 1. This makes the tensions graphically clearer and also shows the difference in the size of uncertainties in the comparison with different experiments. Note that the tensions in the comparison of the ratio to the value 1 are not exactly the same as in the differences between our result and experiment shown in Table III.
It is worth noting that the corrections to C eff,0 7 and C eff,0 9 (discussed in Appendix B) shift our results by amounts of between 0.3 and 1σ, depending on the q 2 region. In the high q 2 region, this is mostly driven by corrections to C eff,0 7 . In the low q 2 region the size and direction of the shift depends on the meson charge. The impact of the corrections on the tension between our SM theory and experiment given in Table III then amounts at most to 0.6σ.
Our numerical results for integration of our differential branching fraction for a variety of q 2 ranges are given in Table VII in Appendix C.
In Section V we discuss the effect on our results, and the tensions with experimental values, of changing the values of the Wilson coefficients C 9 and C 10 that we use in Eq. (4). This tests the impact of possible new physics at a high scale that would change the effective weak Hamiltonian.

Total branching fractions
In this section we discuss the total branching fraction, i.e. dB/dq 2 integrated across the physical q 2 range from q 2 min = 4m 2 to q 2 max = (M B − M K ) 2 . By convention, experimentally measured branching fractions are extrapolated to the full physical q 2 range, ignoring resonances. This is achieved using correction factors based on an assumed q 2 distribution (see e.g. [71]) to the branching fraction integrated over the experimentally measured bins (see discussion in [19]). This achieves the same result as our interpolation across the cc resonance region (discussed in Section II A 3) and this should make our shortdistance theory results comparable to experiment here.
Dashed lines indicate the effect of adding QED uncertainty (see Section II A 5) to our result. Belle '19 BaBar '08 BaBar '08 BaBar '08 The total branching fraction for B → K + − . Our result (HPQCD '22) is given by the black star and grey band, as compared with experimental results [11,12,14,16]. Dashed lines indicate the effect of adding QED uncertainty (see Section II A 5) to our result.
Our result is in general larger than the experimental values, reflecting the pattern in Figures 3, 4 and 5, which do not contain exactly the same experiments. This is clearest in the case of the neutral meson decay, but both the charged and neutral LHCb '14A results with µ in the final state in Figures 7 and 8 show clear tension: 3.9(4.0)σ for B + and 4.2(4.2)σ for B 0 , with(without) QED uncertainty.
To test the effect of the LHCb '14A q 2 interpolation across the vetoed regions, we compare our total branching fraction (without QED uncertainty) with that when we exclude the regions vetoed by LHCb '14A (8 ≤ q 2 /GeV 2 ≤ 11 and 12.5 ≤ q 2 /GeV 2 ≤ 15), we obtain ratios 1.3814(84) and 1.3951 (84) for the B + and B 0 cases, respectively. This compares very well with the factor of 1.39 used by LHCb, and indicates that the interpolation used by them is consistent with our linear interpolation approach. This adds confidence that the comparison between theory and experiment here is a reliable one.
We conclude that, in agreement with what was seen in the low and high q 2 regions in Sec. II B 1, significant tension (exceeding 4σ) is also seen between the total shortdistance SM branching fraction and the most accurate For the decay with τ in the final state, experimental results are much less mature. The BaBar collaboration reported the first measurement of the total B + → K + τ + τ − branching fraction in [13]. Their result of (1.31 +0.66+0. 35 −0.61−0.25 ) × 10 −3 , (where errors are statistical and systematic respectively) is consistent with no signal. It is also consistent with our 7%-accurate SM result of B (+) τ = 1.68(12) × 10 −7 , where we do not include an uncertainty for QED effects since the τ is so heavy. Note that the q 2 range available for the case with the τ is much smaller than for the lighter leptons. We will discuss the τ case further in Sec. II C.
Our numerical results for the total branching fraction for ∈ {e, µ} are given in Table VII in Appendix C. Results for τ in the final state are given in Table X in the same Appendix.

R µ e
We define the ratio of branching fractions for decays to different leptons, We calculate this separately for the charged and neutral meson cases, R µ(0/+) e . We also give the charge-averaged result, R µ e , which is obtained from performing a charge average in each of the numerator and denominator of Eq. (7).
As discussed above, the effect of changing lepton flavour between electrons and muons on our branching Belle '19 BaBar '12 LHCb '21 HPQCD '22 FIG. 10. The ratio of branching fractions to muons and electrons, as defined in Equation (7), with various q 2 limits. Our result (HPQCD '22) has negligible uncertainty compared to experiment [12,14,15,23], and does not differ visibly when integrated over different q 2 ranges. Note that for the ratios of total branching fractions (top two panes), we have used a lower limit on the q 2 integral, q 2 min , that differs between e and µ cases in the ratio. This was done to match what was done in the experiment, but has no visible effect. We include a 1% uncertainty to allow for possible QED corrections [68] as discussed in the text. This completely dominates our theoretical error here. The vertical dotted line is at 1.
Our result (HPQCD '22) is given for bins between integer values of q 2 , with the uncertainty included but too small to be visible. Experimental results [18,20] are included for comparison.
fraction is very small. The factors of β in Eqs. (2) and (3) mean that R µ e < 1 at very small values of q 2 and the factors of m mean that R µ e > 1 at large values of q 2 . However, for almost the entire q 2 range R µ e is very close to 1. This also means that effects from the form factors and from the resonances cancel between numerator and denominator and the theory uncertainties in R µ e are small. Our results are tabulated in Table VIII in Appendix C for a variety of q 2 ranges; we find R µ e to differ from 1 by a few parts in 1000.
Radiative corrections from QED (see Sec. II A 5) af-fect the µ and e cases differently but are well understood [66,68]. Tests show that the corrections applied on the experimental side using PHOTOS [65] should be reliable up to ±1% [68]. We will take a 1% uncertainty on our results to allow for missing QED effects in the comparison of our results with experiment. However, where this uncertainty is already included in the experimental results this will amount to a double-counting of this source of error. This QED uncertainty dominates all other sources of uncertainty in our determination of R µ e . Figure 10 shows a comparison of our result for R µ e (including the 1% uncertainty) with those from experiment. The most notable experimental results are the recent values from LHCb '21 [23] with relatively small uncertainties for the charged meson decay, integrated over the low q 2 region. The LHCb '21 result differs from 1 (taken as the SM result when QED corrections are removed [68]) by 3.1σ. The relative uncertainty from theory is negligible so having improved form factors, as we do here, makes no difference to this tension. When we come to compare with other theoretical work in Section II C, potential deviations of R µ e from 1 will be studied in more detail.

F H
Two additional experimentally measured quantities with which we can compare are based on studying the angular distribution of decay products. The differential distribution in angle is given by [25] 1 Γ where θ is the angle between B and as measured in the dilepton rest frame. Both the flat term, F H /2, and the forward-backward asymmetry, A FB , are small, and sensitive to new physics. A FB = 0 in the SM, up to QED corrections, so we will not consider it here. We can construct F H from [25] Figure 11 shows F H for the B + → K + µ + µ − decay. The black stars show our result, calculated in bins between integer values of q 2 /GeV 2 , with the exception of the first bin, which starts at q 2 min . Experimental results from [18] are included for comparison, as are those from [20]. In the latter case, central values are not given for the data, and so in order to add the systematic error to the statistical confidence interval, this interval is assumed Gaussian, and the central value taken as the middle of the interval. The effect of this approximation is negligible, owing to the small systematic errors, and does not affect the comparison with our results, as the overall confidence interval is broad.
As with R µ e , the cancellation of correlated uncertainties in Eq. (9) means that our result here is very precise compared with experiment (our error bars are included in Figure 11 but not visible). We can see that a large reduction in experimental uncertainty will be required in order to reveal any meaningful inconsistency (if one exists) between SM theory and experiment.
Our results for F H for ∈ {e, µ} are tabulated in Table IX in Appendix C for a variety of q 2 ranges; we will discuss the flat term for the τ case in Section II C.

C. Comparison to theory
In this section we compare to previous theoretical results, both from previous lattice QCD calculations [8,9] and others [25-28, 36, 37, 70]. Again, Table II provides a correspondence between references and labels in figures. Different theory calculations differ both in the form factors they use and in their input parameters used to determine the differential rates. In [3] we compare the improved lattice QCD form factors that we use here with earlier lattice QCD form factors. Here we will show comparisons of differential and total branching fractions, flat terms etc. as for the comparison with experiment in Section II B. There are many more theory results for the case of τ leptons in the final state, so we can make more detailed comparison of that case.

Differential and total branching fractions
We compare differential branching fractions for ∈ {e, µ} in the same manner that we did for experiment, allowing us to factor in the different q 2 bin choices made by different authors. As before, lepton flavour makes an insignificant difference for ∈ {e, µ}, so we treat all cases together. However, we split the cases up by meson charge. The three charge cases (charged B, neutral B and the average of the two) are shown in Figures 12, 13 and 14. Again, note the effect of the charge dependence of the corrections to C eff,0 9 detailed in Appendix B. We see good agreement with previous work, although our uncertainties are generally a lot smaller. We do not agree with the results of AS'12 [36] at low q 2 shown in Figure 14. Our central values are slightly higher than most previous results. In the comparison with FNAL/MILC '15 [9], this ≈ 0.5σ effect can be traced to the larger value of |V ts | that we use, based on more recent B-mixing results (see Section II A 1); we agree well on the form factor central values (see [3]). GRDV '22 [70] fix their form factors from a combination of earlier lattice QCD results (at high q 2 ) and values derived from lightcone sum rules (at low q 2 ). Our form factors have lower central values than theirs but this is partly offset by the difference of values of |V ts | used. Note that GRDV '22 also include non-local form factors in their calculation that account for charmonium resonance contributions. The uncertainty in their differential branching fraction is dominated by that from their local form factors. We compare our results for the differential branching fraction for the case = τ to those of FNAL/MILC '15 [9] and HPQCD '13 [8] in Figure 15. The q 2 range is limited by the much larger value of q 2 min in this case (set by m τ ), and so no interpolation across the vetoed cc resonance region is required. We include both the charged and neutral meson cases on the same plot, in red and blue respectively. We see again that our results are somewhat, but not significantly, higher than the earlier results largely because of the change in |V ts |. Figures 16 and 17 show the total branching fraction for the ∈ {e, µ} and = τ cases respectively. The figures are broken down according to meson charge, with the integral extending from q 2 min to q 2 max , except as indicated otherwise in the plot. In HPQCD '13 [8] the branching fraction to electrons is also given, but (as is the case for us) is insignificantly different from the muon case shown, so we omit it. Again, our central values are somewhat higher than previous work, but not by a significant mar- gin. Our results provide an improvement in uncertainty over earlier values.

R µ e and R τ µ
We now look in more detail at the ratio of branching fractions to different flavours of leptons, defined in Eq. (7). Figure 18 shows our results (as black open stars) for the deviation of the ratio R µ e from 1.0, as a function of q 2 . We use q 2 bins of width 2 GeV 2 , following the pattern 2-4 GeV 2 , 4-6 GeV 2 etc, except for the first bin which covers 0.1-2 GeV 2 . We plot our results for the charged B case; values for the neutral B case are indistinguishable, except in the lowest q 2 bin (numerical values for each bin are given in Table VIII in Appendix C). As discussed earlier, the ratio is very well determined theoretically in a pure QCD calculation; uncertainties from possible QED FIG. 16. The total branching fraction for B → K + − ( ∈ {e, µ}), integrated from q 2 min to q 2 max compared with [8,9,26,37]. Different meson charges are treated separately. The lowest panel shows a comparison to BHDW '11 in which the integration over the differential branching fraction starts at q 2 min = 14.18 GeV 2 .
FIG. 17. The total branching fraction for B → Kτ + τ − , integrated from q 2 min to q 2 max compared with [8,9,26,37]. Different meson charges are treated separately. The lowest two panels show a comparison to BHDW '11 and HPQCD '13 in which the integration over the differential branching fraction starts at q 2 min = 14.18 GeV 2 .
corrections that would affect experiment are not included in Fig. 18. We see that R µ e differs from 1 at the level of 10 −3 , with a larger deviation at large q 2 . There is also a sizeable effect in the smallest q 2 bin where the kinematic cut-off from the µ mass starts to have an effect. Figure 18 also shows a comparison to earlier theory results. Good agreement is seen between the different calculations. Figure 19 shows our results for the ratio R τ µ (Eq. (7)) as open black stars, calculated in q 2 bins of width 2 GeV 2 . The lowest bin is from the top of the cc resonance region at 14.18 GeV 2 (slightly above the kinematic end-point of 4m 2 τ ) to 16 GeV 2 , the next 14-16 GeV 2 and so on. Results are shown for the charged meson case, but the neutral meson results are indistinguishable, see the nu-  FIG. 19. The ratio R τ µ (Equation (7)). Our results (HPQCD '22, open black stars) are evaluated in q 2 bins defined by 14.18, 16, 18, 20, 22 GeV 2 , for the charged B meson case. We compare with HPQCD '13 [8] and FNAL/MILC '15 [9]. No uncertainty is included for QED corrections that would affect experiment. merical values given in Table XI in Appendix C. R τ µ deviates from 1.0 much more strongly than R µ e ; given the lepton masses involved this is not surprising. The results show a similar pattern to R µ e with lepton mass effects pushing R values upwards as q 2 increases. R τ e results would not differ visibly and so are not shown. Again, the comparison to earlier results shows good agreement.
In Fig. 20 we give our result for R µ e using the full q 2 range from 4m 2 µ to q 2 max . This is shown as the deviation from 1 in units of 10 −2 . Note that in forming this ratio we integrate both the numerator and denominator of Eq. (7) through the vetoed region as discussed in Sec. II A 3. We also give our result for R τ µ using the range from 14.18 GeV 2 (the upper edge of the vetoed region) to q 2 max . Both  Open black stars show our results for F e H for the charged B case in bins of q 2 defined by 0.1, 2, 4, 6, . . . , 22 GeV 2 . We compare to HPQCD '13 [8] and FNAL/MILC '15 [9]. of our results are consistent with zero deviation from 1 at the 2σ level. We compare to HPQCD '13 [8] and see reasonable agreement; our result for R µ e has half the theory uncertainty of their earlier value.  effects that would affect experiment are included in the plots.
For the light, e and µ, lepton cases F H is small across most of the q 2 range but rises rapidly towards the low q 2 end. For small values of m and q 2 the integrand in the numerator of F H (Eq. (9)) is proportional to β(1 − β 2 ) (see Eq. (2)). For q 2 = O(4m 2 ), this has no suppression by m 2 and gives a large contribution to the integral. This means that great care must be taken in integrating the numerator in the region q 2 → 4m 2 (see Appendix C 1).
In Figures 21, 22 and 23 we also compare to earlier (mainly lattice QCD) results, plotting the values for the charged meson case where both are given. There is good agreement between the different calculations in general; our results have an improvement in uncertainty. We see some tension (2.8σ) for F e H and F µ H with results from FNAL/MILC '15 in the smallest q 2 bin for the charged B case, but agree well for the neutral B case. This can be seen by comparing our Table IX with Tables IX and X FIG. 25. F µ H , integrated over the full q 2 range, for the charged and neutral meson cases, and over a reduced range from 14.18 GeV 2 in the case of the average of the charged and neutral cases (lowest panel). Our results are shown as the filled stars and denoted 'HPQCD '22'. This work is compared with [9,26], shown as filled triangles and circles respectively. of [9]. For F τ H we see some tension with HPQCD '13 [8]. Figures 24, 25 and 26 show F H as defined above, integrated over the full q 2 range q 2 min ≤ q 2 ≤ q 2 max , except in cases where the range is stated explicitly. These plots are split up according to meson charge, although in most cases the difference is insignificant. A significant difference is seen for F e H between the charged and neutral cases when integrating from 4m 2 e up to q 2 max . This difference comes from the nonfactorisable corrections to C eff 9 discussed in Appendix B. These depend on meson charge and have a sizeable impact at very small value of q 2 , relevant to F e H (see below) but not for other lepton flavours or to other quantities. We FIG. 26. F τ H , integrated over the full q 2 range, for the charged and neutral meson cases, and over a reduced range from 14.18 GeV 2 in the case of the average of the charged and neutral cases (lowest panel). Our results are shown as the filled stars and denoted 'HPQCD '22'. This work is compared with [8,9,26], shown as filled hexagons, triangles and circles.
F τ H , with our results showing a considerable improvement in uncertainty in most cases. For F e H our results are in tension with those of FNAL [9] when we integrate from 4m 2 e up to q 2 max . Given the good agreement seen for most of the q 2 range in Fig. 21 it seems likely that the tension is a result of behaviour in the integrand of the numerator towards the lower limit of the q 2 integral, discussed above. We illustrate this in Fig. 24 by showing results also for the case where the lower cut-off on the q 2 integral is 4m 2 µ rather than 4m 2 e . Doing this gives considerably lower values for F e H because the integral in the numerator of F e H changes; the denominator is half the total branching fraction and very insensitive to small changes in q 2 min . We conclude that determining F e H at small q 2 values or integrated over the full kinematic range requires care over integration and is sensitive to small-q 2 corrections. However, away from small q 2 values there are no issues with determining F e H (see Figure 24) and this is the experimentally more accessible region.
We also investigate the rare lepton flavour number violating (LFV) B → K − 1 + 2 decay, with different final state leptons 1 = 2 (with masses m 1 and m 2 ). The differential branching fraction for this decay is given in Eq. (9) of [72] in terms of a number of Wilson coefficients and functions ϕ i (q 2 ) (Eq. (10) of [72]). The knowledge of these functions in the SM can be combined with experimental measurements to set limits on the Wilson Coefficients [73].
In the case that the LFV in B → K − 1 + 2 arises from purely vector or purely scalar operators, we can write simple forms for the branching fraction, in terms of P and lepton flavour dependent coefficients a 12 K , b 12 K , e 12 K and f 12 K . 10 9 B(B → K 1 2 ) = a 12 K |C 9 +C 9 | 2 +b 12 K |C 10 +C 10 | 2 (10) for the vector case and 10 9 B(B → K 1 2 ) = e 12 K |C S +C S | 2 +f 12 K |C P +C P | 2 (11) for the scalar. The factor of 10 9 makes the coefficients O(10).
where the ±s are matched with 9, 10 and S, P . We use our improved f + and f 0 form factors [3] to determine these functions, including uncertainties for strong-isospin breaking discussed in Section II A 4. For the MS s mass in Eq. (13) we take m MS s (µ b ≡ 4.2 GeV) = 0.07966(80) GeV, run from m MS s (µ = 3 GeV) = 0.08536(85) GeV [74]. m MS b (µ b ) is given in Table I, along with the other parameters used. None of the ϕ functions are sensitive to interchanging 1 ↔ 2, so neither are the resulting coefficients (a 12 K = a 21 K , etc.). We plot the functions φ i (q 2 ) ≡ |N K (q 2 )| 2 ϕ i (q 2 ) in Figures 27 and 28, for the choice 1 = µ, 2 = τ , and find good agreement with similar plots in [72]. Note that there is no issue with cc resonances in this case. The functions are insensitive to B meson charge, so the neutral case is shown here -the charged case differs mainly by the overall multiplicative factor of τ B 0/+ used in Eq. (14). Plots for lepton combination eτ are nearly identical to those for µτ , and are not shown.
The shape of φ S,P (q 2 ) varies little with lepton flavour choices, so we show only the one example in Fig. 28. The shape of φ 9,10 is more sensitive, however and we plot the case of eµ in Fig. 29. It shows the functions falling much more abruptly towards the minimal q 2 = (m 1 + m 2 ) 2 value than when one lepton is a τ .
Integrating φ 9,10,S,P across the full physical q 2 range for each 1 , 2 combination yields a 12 K , b 12 K , e 12 K and f 12 K respectively, and we report their values in Table IV. These are updates to the values in [72] (which took LCSR form factors from [38]). We agree well with the results found in [72], with a modest improvement in a 12 K and b 12 K uncertainties. Table IV gives results for both charged and  [15]. These limits start to provide constraints on beyond the Standard Model (BSM) physics, for example models with leptoquarks [75]. Current experimental limits for B + → h + τ are a few times 10 −5 [30] and so need to be reduced by several orders of magnitude to provide a test of this decay mode.

IV. B → Kνν
Lastly, we study the rare decays B → Kνν, which are of phenomenological interest [76]. Precise measurements in this channel could shed light on new physics models to explain lepton flavour universality violation [77,78].
There is limited experimental data at the present time, and currently no experimental results have conclusively demonstrated a non-zero branching fraction. However, more experimental results are anticipated in the near future, and the upper limit on the branching fraction may soon be driven down close to theoretical determinations [35].
We can calculate the short distance (SD) contribution to the differential branching fraction (summed over neutrino flavours) [9,45,79], which depends only on the vector form factor.
with X t = 1.469(17) [80], sin 2 θ W = 0.23121(4) [43] and 1/α EW (M Z ) = 127.952(9) [43]. Other parameters are given in Table I. We evaluate this expression using our improved vector form factor determined in [3], including uncertainties for strong-isospin breaking discussed in Section II A 4. If the reader wishes to reproduce this, f + (q 2 ) can be easily obtained using the code attached to [3]. Note that no QED uncertainty is required in this case.
Our results for the differential branching fractions for the short-distance contribution are shown in Figure 30 for both the charged and neutral B meson cases. Note that there is no issue with cc resonances or nonfactorisable terms in this case so the SM calculation is very straightforward and simply requires an accurate vector form factor. From the differential branching fraction we calculate the total branching fractions, adding in the long-distance contribution to the charged case. These values are given in Table V, and show a total uncertainty of around 7%. Of this uncertainty 5% comes from the vector form factor, 4.4% from the CKM elements |V tb V * ts | and 2.3% from X 2 t . Table V also compares our values to experimental limits and earlier theoretical results. Our results sit slightly higher than previous theoretical determinations and have smaller uncertainty. The comparison of our B + → K + νν branching fraction result with other theoretical values is shown in Figure 31. Results in various q 2 bins are presented in Table VI. Current experimental bounds on the B → Kνν branching fraction are roughly an order of magnitude larger than the SM expectation given by the theory results of Table V. The uncertainty that we have achieved here, however, complements the ≈ 10% error expected from Belle II with 50 ab −1 [35], allowing for a more stringent test of this quantity in future. This is demonstrated in Figure 32, which compares our result with the 90% confidence limits set by BaBar [34] and Belle II [33] in Table V, as well as the precision forecast by Belle II with 5 ab −1 and 50 ab −1 [35], assuming a result centred on our value.  [37] 3.98 (47) [45] 4.94 (52) [9] 4.53 (64) [86] 4.65 (62) [87] 5.58 (37) HPQCD '22     ) for B(B + → K + νν), summing short-and long-distance effects, compared with current 90% confidence limits from BaBar [34] and Belle II [33] (blue lines). Red filled squares with error bars show the expected precision to be achieved by Belle II with 5 ab −1 and 50 ab −1 [35], assuming our result is obtained for the central value.

V. CONCLUSIONS
We have used our improved scalar, vector and tensor form factors calculated in N f = 2+1+1 lattice QCD [3] to determine SM observables for B → K + − , B → K − 1 + 2 and B → Kνν decays. The form factors were calculated in a fully relativistic approach which enables the full q 2 range of the decay to be covered. This approach also gives better control of the current normalisation than in previous calculations. Our form factors then have smaller uncertainties, particularly at low q 2 , than earlier work so that we can improve QCD uncertainties on the SM observables. We have tabulated all of our results across a variety of q 2 bins in Appendix C, for future use.
For B → K + − we determine the differential rate for both charged and neutral meson cases, including an analysis of nonfactorisable corrections at low q 2 (see Appendix B). We then compare to the wealth of experimental information (see Figures 3, 4 and 5) focussing on the low and high q 2 ranges below and above the region where cc resonances make a large contribution to experimental results, but are not included in our calculation. We allow an uncertainty for QED effects missing from our calculation in this comparison.
Previous lattice QCD calculations by the Fermilab/MILC collaboration [9] saw a 2σ tension between the SM branching fraction and LHCb results [19] in the low and high q 2 regions, with the SM results being higher than experiment in both regions. Our improved form factors significantly sharpen this tension, which becomes particularly strong in the region below the charmonium resonances. For 1.1 < q 2 < 6GeV 2 we find a tension of 4.2σ for B + → K + µ + µ − and 3.1σ for B 0 → K 0 µ + µ − with LHCb '14A [19] and 2.7σ for B + → K + e + e − with LHCb '21 [23] (see Table III and Fig. 6). We see no tension with results from Belle '19 [15] in this low q 2 region, but they have much larger uncertainties. The tension between our results and LHCb in the high q 2 region (15 < q 2 < 22GeV 2 ) amounts to 2.7σ, again with our SM results being higher than experiment.
The tension between the SM and LHCb results in the 'safe' q 2 regions away from narrow cc resonances points to the possibility of BSM physics that would have the effect of changing the Wilson coefficients of the effective weak Hamiltonian. Recent analyses (see, for example, [70,88]) focus on BSM changes to C 9 and C 10 and include a number of different B decay processes in providing a best fit result for these BSM effects. Gubernari et al [70] give values ∆C BSM 9 = −1.0 and ∆C BSM 10 = +0.4 at µ = 4.2 GeV from comparison of their SM results to experiment in three separate analyses, one of which is a combined one using B → Kµ + µ − and B → µ + µ − and the other two are for B → K * µ + µ − and B s → φµ + µ − . If we repeat our analysis changing C 9 and C 10 by these BSM shifts we see significantly lowered values for the differential rate for B → K + − because C 9 and C 10 now have smaller magnitudes. We find the tension with LHCb '14A in the low q 2 region (1.1-6 GeV 2 ) to be reduced from 4.2σ to 1.1σ. All of the other tensions seen in Table III where the SM result exceeds the experimental value are now reduced below 1σ. At the same time a 2.2σ tension with the Belle '19 result for B + → K + µ + µ − develops in the opposite direction. We conclude that BSM shifts to C 9 and C 10 of -1.0 and +0.4 respectively give improved agreement between our results and the experimental values for the B → K + − branching fractions.
We update results for ratios of B → K + − branching fractions to leptons of different flavour, , giving results for q 2 bins across the kinematic range in Section II C 2. The 'flat term' in the lepton angular distribution is also examined in detail in Section II C 3 for different .
We have also provided improved results for the formfactor-dependent pieces that, multiplied by appropriate Wilson coefficients, give the branching fraction for the lepton flavour number violating decay B → K + 1 − 2 . Table IV gives our results for all 3 lepton flavour pairs and for charged and neutral B mesons, with uncertainties of 6-7%. These can be used to constrain BSM models in future searches for this decay.
Finally, Table V gives our improved SM values for B(B + → K + νν) and B(B 0 → K 0 νν). We have uncertainties below 10%, commensurate with that expected from Belle II with 50 ab −1 worth of data. This decay mode, expected to reach 3σ significance at Belle II with 5 ab −1 of data, is an exciting possibility for future tests for BSM physics, since it is theoretically 'cleaner' than B → K + − across the full q 2 range with no need for vetoed regions.
We have shown here the improvements to B → K semileptonic decay phenomenology that result from our improved lattice QCD form factors covering the full kinematic q 2 range. For B → K + − our relative uncertainties are comparable with those from experiment. Our uncertainties are dominated by those from the form factors so further improvements in the lattice QCD calculation of the form factors will reduce them further. This is possible in future with increased computational power enabling smaller statistical errors and a push towards smaller values of the lattice spacing reducing residual discretisation effects. A more detailed comparison will then be possible with future experimental results with smaller uncertainties for B → K + − and future differential rates for B → Kνν. In combination with analysis of other decay processes for hadrons containing b quarks, this is a promising way forward to uncover physics beyond the Standard Model.

VI. ACKNOWLEDGEMENTS
We thank D. van Dyk, P. Owen and E. Lunghi for clarification regarding data in their papers and Rafael Coutinho for help understanding the use of PHOTOS in LHCb results. We also thank L. Grillo and R. Zwicky for useful discussions, as well as Quim Matias and Sébastien Descotes-Genon. Computing was done on the Cambridge Service for Data Driven Discovery (CSD3) supercomputer, part of which is operated by the University of Cambridge Research Computing Service on behalf of the UK Science and Technology Facilities Council (STFC) DiRAC HPC Facility. The DiRAC component of CSD3 was funded by BEIS via STFC capital grants and is operated by STFC operations grants. We are grateful to the CSD3 support staff for assistance. Funding for this work came from STFC.
Note added in proof -Whilst this paper was being prepared for publication, a set of new results for R K and R * K appeared from LHCb [89]. These results agree with Standard Model expectations (see Sec. II B 3) and supersede their earlier values.

Appendix A: Calculating pole masses
The three loop relation between quark masses in the MS scheme and the pole mass scheme is [60] m(m) m  Table I. Plugging these results into Equation (A2) results in an updated value for m. The initial guess for m is adjusted to reduce the difference between it and the value obtained from Equation (A2). This process is repeated until the values of m converge.
Using this method, we obtain the pole masses m c = 1.684 (22) GeV and m b = 4.874(32) GeV. The perturbation series in this expression suffers from the presence of a renormalon in the pole mass [61], so we take a 200 MeV uncertainty on both numbers. The effect of this uncertainty is described in Section II A 1. We note that 4m 2 c = 11.34 GeV 2 falls within the vetoed charmonium resonance region.
where C eff,0 7 and C eff,0 9 are given in Table I. The corrections are discussed in Appendix B of [9], which compiles results from [26,27,64,[91][92][93][94][95][96][97][98][99][100][101]. We do not repeat formulae from that reference here, but direct the reader there for more detail; our aim here is to assess the significance of these corrections to the decay rates. Below we outline the form of the corrections, plot them versus q 2 and discuss their sizes relative to C eff,0 7 and C eff,0 9 . All numerical inputs not explicitly stated below can be found in Table I.
The leading contribution to the correction δC eff 7 is from O(α s ) effects, 1,c + C 8 F where the expression for F 1,c is lengthy and provided in the C++ header files of [92]. F is given in Appendix B of [9]. We use α s (4.2 GeV) = 0.2253 (28), which is run from α s (5.0 GeV) = 0.2128 (25) [90]. The next higher order contribution, which we neglect, is suppressed by a factor of λ  The leading order contributions to the correction δC eff 9 are given by We neglect the O(α s λ (s) u ) term, which is even smaller (see Equation (B11) of [9] for more details). The function h(q 2 , m) is defined in Equation (6) and F (9) 8 is given in Appendix B of [9]. Expressions for F (9) 1,c and F (9) 2,c are also provided in the C++ header files of [92]. The corrections δC eff 7 and δC eff 9 are applicable across the full q 2 range. We plot the real and imaginary parts of δC eff 7 and δC eff 9 in Figure  , which is approximately 4 across the full q 2 range (Table I). The O(α s ) contributions to δC eff 9 peak at q 2 ≈ 10 GeV 2 , due to the behaviour of the functions F 1,c , F 2,c , F 1,c and F (9) 2,c . This peak occurs within the experimentally vetoed J/Ψ resonance region and is largely contained within the region of q 2 between the J/Ψ and Ψ(2S) resonances, outside of which the contributions are modest. As a result, it has minimal impact on results in the well behaved regions of q 2 below the J/Ψ and above the Ψ(2S). The uncertainty in Re[δC eff 7 ] and Re[δC eff 9 ] grows rapidly towards q 2 max because of the behaviour of F (7/9) 8 [9]. This effect is suppressed in observables by the fact that the differential decay rate vanishes at q 2 max and so it is not noticeable in plots of observables versus q 2 or in uncertainties of observables in the largest q 2 bins.
Non-factorisable corrections are accounted for via ∆C eff 9 , which is dependent on the meson charge (see be- low). Following the notation of [9], where C F = 4/3, and T K,+ = 0. The functions Φ B,± (ω), Φ K , and T (0/nf) K,± are given in [91] and Equation (B5) is discussed in detail in Appendix B of [9]. We evaluate the expressions for Φ B,± (ω), Φ K , and T (0/nf) K,± using the following inputs: ω −1 0 = 3(1) GeV −1 [91], a K 1 = 0.0453(30), and a K 2 = 0.175 (50) [102] (we take Φ K to second order). The non-factorisable corrections are valid for small q 2 , and we turn them off at q 2 = 8.68 GeV 2 , the start of the vetoed J/Ψ resonance region. We do not calculate through this region. Instead, we linearly interpolate, beginning from the point where the corrections are turned off through the Ψ(2S) resonance, so the differential branching fraction is a smooth function of q 2 (see, e.g., Figure 13).
The contribution from ∆C eff 9 to observables is via the term f + ∆C eff, 1 9 in Equation (4). This has the effect of cancelling the dependence of f + ∆C eff 9 on the form factor f + . The real and imaginary parts of the three nonzero contributions to f + ∆C eff 9 in Equation (B5), corresponding to T In these plots we remove a decay channel-specific factor of the light quark charge, e q , which is 2/3 for B + → K + and −1/3 for B 0 → K 0 . Among these terms, the T The combined effect of these terms to f + ∆C eff 9 is shown in Figure 36, where the real and imaginary parts are plotted separately for both the B 0 and B + cases. Both the real and imaginary parts are smooth functions of q 2 in the region below the J/Ψ resonance where they are considered (4m 2 ≤ q 2 ≤ 8.68 GeV 2 ), and are small for q 2 > 1 GeV 2 . For q 2 < 1 GeV 2 , where the nonfactorisable corrections are largest, they also have little impact. This is because the differential decay rate falls rapidly as q 2 approaches m 2 l for kinematic reasons (see the factor of β = 1 − 4m 2 /q 2 in Eq. (3)). As a result, the corrections do not make a significant contribution to the differential branching fraction.
The overall modest contributions of the above corrections are shown for C eff 7 in Figure 37 and for C eff 9 in Figure 38. In each plot, both corrected and uncorrected values are plotted. In the case of C eff 9 , the difference between the B 0 → K 0 and B + → K + corrections, arising from the charge dependence of ∆C eff 9 , is responsible for the slightly different shapes of the differential branching fractions (see e.g. Figure 3 and 4). Cusps in C eff 7 and C eff 9 either occur within the vetoed J/Ψ resonance region near q 2 ≈ 10 GeV 2 or at q 2 = 4m 2 c ≈ 11.3 GeV 2 , between the J/Ψ and Ψ(2S) resonances. In this region    (Table I)  between the resonances, we linearly interpolate the differential decay rate and are therefore unaffected by the cusp at q 2 = 4m 2 c . The most significant effect of the corrections is an approximately 20% shift to Re[C eff, 1 7 ] arising from Re[δC eff 7 ]. In fact the total effect of all corrections on our results in the well-behaved regions is small (Section II B), so we do not include additional uncertainty on the corrections. The growth of corrections and in their uncertainties at high and low q 2 is suppressed by kinematic factors in the decay rate, resulting in modest impact in the well-behaved regions of q 2 where we give our main results (Table III), as discussed.

Appendix C: Numerical Results
The Tables in this appendix list our numerical results in a number of q 2 bins, aimed at matching those appearing in previous work and providing good coverage for future work to compare with. The charged and neutral meson case, as well as the average (no charge label) is given. Similarly the e, µ, average of e and µ (labelled as ) and τ cases are all given. In all cases, the average is taken (including correlations) after all other calculation and integration has taken place. In most cases the meson charge and/or light lepton flavour makes no significant difference. We include all cases for completeness, and ease of comparison with future results.

Integration
When performing integrals to generate the data in the following tables, as well as to generate the figures above, we use the trapezoid rule. That is to say, where f is evaluated at N + 1 evenly spaced q 2 values in total. In order to ensure accurate results, we start with N = 16, and double N repeatedly until the result for 2N differs from that for N by less than 0.02 σ (as uncertainties are given to two significant figures, this is the point at which on average the last digit of the mean has only changed by one). Using this method allows us to confirm that our integrals have converged, and most do so for N ≤ 1024. An exception is F e H . In this case, the integrand in the numerator (Eq. (9)) changes very rapidly as we approach the lower limit 4m 2 e ≈ 10 −6 GeV 2 and a naive approach will not allow the integrals to converge without very large values of N . This is not the case for F µ H , as the integrals do not span so many orders of magnitude for 4m 2 µ ≈ 0.05 GeV 2 . To solve this problem, we split the integral limits up into different orders of magnitude, integrating from q 2 low = 4m 2 e to 10q 2 low , then from 10q 2 low to 100q 2 low and so on, up to q 2 upp . In each case we increase N until the result is stable, as above, and then we sum the (correlated) sub-integrals to get the final value for both the numerator and denominator. This method concentrates the number of evaluations of f in regions where they are most needed, and so reduces the total number of evaluations required for the whole integral, making it tractable.
To assess the suitability of this method, and check that we are not introducing bias with our algorithm, we checked the integration against the vegas python package [103], and confirmed that the results agree.    7) integrated over some commonly used q 2 binning schemes. We give results for the charged meson case, the neutral meson case and the charge-averaged case (defined as the ratio of the two charge-averaged integrals). Numerical values for q 2 bins are in units of GeV 2 . The 1% uncertainty from QED effects [68] (see Sec. II B 3) is not included in these numbers.    (7)) for the τ case to the electron, muon and , (which is the average of the two) cases, integrated over some commonly used q 2 bins. Numerical values for q 2 bins are in units of GeV 2 . No uncertainty to allow for QED effects is included in these numbers.   (9)) integrated over some commonly used q 2 bins for charged and neutral B meson decays. Numerical values for q 2 bins are in units of GeV 2 . No uncertainty to allow for QED effects is included in these numbers.