B-physics from Lattice Gauge Theory

We discuss the main issues in dealing with heavy quarks on the lattice and shortly present the different approaches used. We discuss a selection of computations covering first the b-quark mass and the B(s) meson decay constants as the consolidated results (neglecting isospin breaking corrections). In the second part we consider recent calculations of form factors for tree-level semileptonic decays with emphasis on the tensions between the results produced by different collaborations. We propose benchmark quantities and tests suited to investigate the origin of such tensions. Finally, we review computations of the bag parameters parameterising neutral meson mixing and provide an overview on a few recent developments in the field.


Introduction
The large b-quark mass and the comparably long lifetime of B (s) mesons, allow for a plethora of experimental observables which can be measured precisely (see for example Refs.[1] and [2] for a discussion of the opportunities at LHC and Belle II, respectively).If precision Standard Model (SM) predictions are available, these can be used to test the SM as well as to search for and constrain New Physics (NP) beyond the Standard Model (BSM).In addition to direct comparisons between experiment and theory, the self-consistency of the SM can be tested by over-constraining the CKM matrix [3,4] and testing its unitarity.
Lattice QCD is the only known tool to provide non-perturbative, ab initio, systematically improvable precision predictions.In lattice QCD a finite Euclidean space-time volume (L) is discretised (with lattice spacing a) and a representative ensemble of the gauge field configurations is sampled via Monte Carlo methods.On these configurations, correlation functions are computed from which hadronic masses and matrix elements can be extracted.This is repeated for multiple choices of the simulation parameters (a, L and the quark masses am q ) and to make contact with experiment the observables are inter/extrapolated to the physical world (a → 0, L → ∞, am q → am phys q ).It is noteworthy, that since lattice QCD simulations take bare quark masses as inputs, predictions at unphysical quark masses can be made, which can be used to test effective field theories (EFTs) such as Chiral Perturbation Theory (χPT) or Heavy Quark Effective Theory (HQET).Modern simulations include contributions from sea effects of two degenerate light quarks and the strange quark (N f = 2+1) and the charm quark (N f = 2+1+1).
In recent years lattice QCD calculations of hadronic observables containing a b quark have made huge advances with precision predictions with all systematic uncertainties under control for several quantities. 1 Here, we comment on some particular challenges when simulating heavy quarks (Sec.2) before summarising the status of computations of the b-quark mass and leptonic decay constants (Sec.3), semileptonic decay form factors (Sec.4) and neutral meson mixing (Sec.5).In Sec. 6 we comment on recent, more exploratory developments before concluding in Sec 7.

Lattice challenges for heavy quarks
Including heavy quarks, such as the b quark, in lattice simulations of QCD poses a multi-scale problem and currently all approaches still rely on the use (even though possibly at different stages) of EFTs, typically HQET [6,7] or Non-Relativistic QCD (NRQCD) [8].The infrared scale is set by the lattice extent L and the dynamics of the light degrees of freedom, such as the pions, should not be distorted by the finite size of the lattice.The standard requirement is to have m π L > 4, with m π the pion mass.The ultraviolet scale is instead set by the lattice spacing a, the finest resolution in the system.In order to properly resolve the propagation of the heavy degrees of freedom, such as the b quark, the mass m b should be far from the cutoff 1/a or in other words am b should be smaller than one.Substituting the physical value for m b , one concludes that lattices with L/a significantly larger than 100 are needed to keep both finite size and discretisation (or cutoff) effects under control.While such simulations are becoming feasible with current machines and algorithms, in order to perform a controlled continuum extrapolation (a → 0) at the b-quark mass one will need simulations at even finer resolutions.This might require new algorithms in order to ensure an ergodic sampling of the configuration space (see Ref. [9]).
In this context EFTs are used essentially in two, non-exclusive, ways.Either they are implemented directly on the lattice by (formally) expanding physical quantities in inverse powers of the scale associated to the heavy degrees of freedom (e.g., the b-quark mass), or they are used to extrapolate results obtained for heavierthan-charm-quark masses to the b-quark mass.In a hard-cutoff regularisation, such as the lattice regularisation, the first approach often produces power-divergences (divergences in inverse powers of the lattice spacing), due to the mixing between operators of different dimensions.Those indeed naturally appear in an expansion in terms of couplings of negative dimension such as 1/m b .As pointed out in Refs.[10,11], these power divergences need to be removed non-perturbatively if one wants to perform a continuum limit extrapolation.
In practice one considers two EFTs, one being HQET (or NRQCD) and the other the Symanzik Effective Theory (SyEFT) [12,13], therefore implementing the O(a) improvement programme to systematically reduce cutoff effects.In SyEFT higher-dimensional operators are introduced depending on the symmetries of the lattice action.Dimensions are compensated by powers of the lattice spacing and the coefficients are tuned in order to remove the leading lattice artifacts.Roughly speaking, different approaches implement the two EFTs in different order.
In O(a)-improved HQET, one considers the static Lagrangian density for the heavy quarks.This locally preserves the heavy flavour number and is symmetric under continuous SU(2) rotations in Dirac space where ⃗ ϕ is the transformation parameter.Such symmetries are only realised in the infinite mass limit; they are not symmetries of QCD and are broken by 1/m b corrections.At the static order a term δm ψ(x)ψ(x) can be added to the Lagrangian with δm being a power divergent (as 1/a) parameter producing a shift in the energy levels of heavylight systems.This is the first occurrence of the power divergences mentioned above.The subtraction of this divergence can be performed through a non-perturbative matching between HQET and QCD, as discussed in Ref. [14] for the example of the computation of the b-quark mass.
Higher orders in 1/m b appear with operators of dimension 5 and larger.Including those in the Lagrangian would produce a non-renormalisable theory and therefore such corrections are treated order-by-order as space-time volume insertions in correlation functions computed in the static theory (see Ref. [15] for a general discussion, applied to the calculation of the b-quark mass including 1/m b corrections).Similarly, when considering matrix elements, local operators have a static expression and higher order (and higher dimensional) corrections that appear with appropriate powers of 1/m b .Those have to be treated together with the higher order terms in the Lagrangian and a framework for a complete and non-perturbative matching of the action and the vector and axial heavy-light currents between QCD and HQET has been put forward in Ref. [16].
The Symanzik improvement programme is implemented in a straightforward way for the action and operators.As only the static action is directly simulated (higher orders being treated as insertions) the symmetries of the static theory can be used in classifying the mixing between operators of different dimensions for the O(a)improvement.
Concerning the quantities reviewed here, results from non-perturbative HQET are available for the B s → Kℓν form factors [17,18], the b-quark mass [15,19], B (s) -meson decay constants [11,20,21] and mixing parameters [22,23].Since most recent results use different formulations for the b quark, we will not further discuss this approach here.
By interchanging the order in which the Symanzik expansion and the Heavy Quark expansion are performed, one is led to the relativistic heavy-quark formulation such as the one introduced in Ref. [24], which is now known as the Fermilab action.The main property is that the coefficients of the higher dimensional operators appearing through the Symanzik improvement programme are allowed to depend explicitly on the heavy-quark mass m h .In this way the relativistic heavy-quark actions, at fixed lattice spacing, interpolate between the massless limit and the infinitemass (static) limit.In particular, this implies that for am h ≫ 1 one expects to recover the power-like divergences of lattice HQET.The starting point is the anisotropic clover action density where the anisotropy parameter ζ, the clover coefficient c SW and the mass parameter m 0 are tuned to reproduce experimental values for B-mesons spectral quantities such as the dispersion relation and the vector-pseudoscalar splitting [24][25][26].The heavy-quark symmetries emerge naturally in the action in eq. ( 3) and therefore HQET can be used to model and estimate cutoff effects [27][28][29].A variant of the Fermilab action, where the parameters m 0 a, ζ and c P (a generalised version of c SW ) are tuned non-perturbatively [30], is known as the Relativistic Heavy Quark (RHQ) [26,31] action.Other effective approaches have been used to treat heavy quarks on the lattice.NRQCD is formally an expansion of QCD in powers of the heavy-quark velocity v and a lattice version has been introduced in Refs.[32,33].Leading operators as well as operators suppressed by O(v2 ) are included in typical applications together with operators needed for O(a 2 )-improvement.The expression for the lattice action is rather lengthy and can be found in Ref. [33].Operators up to dimension 7 appear at the order mentioned above.Dimensions are compensated by powers of the heavy-quark mass and the coefficients are determined through matching with QCD, typically performed at tree-level or at most at one-loop. 2  NRQCD is hence clearly non-renormalisable by power counting and the continuum limit at fixed heavy-quark mass does not exist.The accuracy that can be reached depends on the existence of a window in the lattice spacing where both cutoff effects (positive powers of a) and power-like divergences (negative powers of a) are under control.In actual computations the precision is at the percent level and the approach is becoming less and less used, since it would be computationally quite demanding to improve it (by including operators of even higher dimension and performing the matching at high loop orders).
Finally, in the most recent applications, heavy quarks are discretised using regularisations originally introduced for the light flavours.This is often referred to as fully relativistic formulations and is possible because of the very fine lattice spacings that can be reached in modern simulations (around 0.04 fm) and because of the use of highly improved actions, where mass-dependent cutoff effects start at O((am h ) 2 ).Since the bquark mass typically only satisfies am b ≲ 1 for the finest lattice spacing, HQET is still used for the simultaneous continuum and heavy-mass extrapolations.
The HPQCD Collaboration has introduced the use of highly improved staggered quarks (HISQ action) [35] at very fine lattice spacings and for b quarks in [36].The first study concerned the B smeson decay constant but by now the method has been applied to a variety of quantities that will be discussed in the following.
In a less direct approach, the ETM Collaboration has been using automatically O(a)-improved twisted mass fermions with masses in the heavierthan-charm region, together with results in the static approximation to interpolate to the bquark mass.The main implementation is through the ratio method [37], where suitable ratios of heavy-light quantities, for different heavy-quark mass, are constructed such that they possess a well defined static limit (typically equal to 1).B-physics quantities are then obtained as an interpolation (rather than an extrapolation) between the static value and the results of the simulations around the charm.
Results from the use of relativistic heavy-quark actions or fully relativistic actions represent the majority of recent results, including those that we are going to review here in the remaining part of this contribution.
3 The b-quark mass and leptonic decay constants

b-quark mass
The b-quark mass is a parameter of QCD and hence, from the theoretical point of view, knowing its value is of fundamental importance.Lattice determinations of the b-quark mass are three times more precise than continuum ones [38] and they drive the accuracy on the parameter.In the latest HPQCD21 F/M/T18 Gambino17 ETM16 HPQCD14B HPQCD14A FLAG21 Fig. 1 Summary of recent results for the b-quark mass in the MS scheme for calculations including the dynamical effects of 2 + 1 + 1 flavours.Individual results shown are HPQCD14A [40], HPQCD14B [41], ETM16 [42], Gam-bino17 [43], F/M/T18 [44], and HPQCD21 [45].The magenta band corresponds to the recommended value from the FLAG collaboration [39] based on these results.edition of the FLAG review [39] the values for the mass in the MS scheme (illustrated as the magenta band in Fig. 1), and for the Renormalisation Group Invariant (RGI) mass, are quoted using results with 2 + 1 + 1 dynamical flavours.As a remark, in the RGI case the error (which is around one percent) is dominated by the RG-evolution required in the definition and therefore by the uncertainty on the Λ-parameter of QCD.
The most recent computations of the b-quark mass entering the FLAG average for N f = 2 + 1 + 1 are shown in Fig. 1.These results use the ratio method described above [42,43] or tune the heavy-quark mass, possibly in some conveniently defined intermediate scheme, in order to reproduce (or extrapolate to) the B (s) -meson mass [44] or extract the b-quark mass from the analysis of moments of heavy current-current correlation functions [40,41,45].The latter method is rather novel and has already produced some of the most precise results and hence deserves a more detailed discussion.The idea is first introduced in Ref. [40] and consists of computing moments of the zero-momentum two-point function G(t) of the heavy pseudoscalar density j 5 = am 0h ψ h γ 5 ψ h , where am 0h is the bare heavy-quark mass.More precisely, in order to reduce discretisation errors, one computes the ratios Rn defined as for n ≥ 6 , (7) with G (0) n being the lowest perturbative order in the expansion of the correlation function.The key observation is that for low values of n the moments are short-distance dominated and they become more and more perturbative as the quark mass is increased.By matching the lattice results to the continuum perturbative prediction, one can simultaneously extract the heavy-quark mass in the MS scheme and the strong coupling α s .The coefficients in the expansion of Rn , for low n, are known to order α 3 s included [46][47][48][49][50], and the method can therefore provide rather accurate estimates of the strong coupling and the b-quark mass.
From the phenomenological point of view a precise knowledge of the b-quark mass is required for testing the properties of the SM-like Higgs particle.Its dominant decay mode is bb, with a branching ratio of about 60%.The theoretical estimate of that depends quadratically on the b-quark mass, which therefore contributes about 2% to the total uncertainty.The signal strengths (ratio of experimental to theoretical estimates) based on CMS and ATLAS data [51][52][53] are consistent with the SM with an error of about 10%.The PDG [38] averages the results to a signal strength of 0.99 (12), the error being currently dominated by the experimental uncertainties.Summary plots for the different production channels, taken from Refs.[51,52] are shown in Fig. 2.

Decay constants
Within the Weak Effective Hamiltonian framework, and neglecting electromagnetic interactions, the decay rate for tree-level process B + → ℓ + ν ℓ can be expressed as Similarly, the decay rate for the loop-mediated process B 0 q → ℓ + ℓ − (q = d, s) depends on the product |V * tb V tq | 2 and the decay constant f 2 Bq .These hadronic parameters are defined by the QCD matrix element (q = u, d, s, c) where the left-hand-side of the first equation is exactly what is computed on the lattice.By combining experimental measurements with results for the decay constants one can therefore obtain exclusive determinations of CKM matrix elements such as |V ub |.The measured channel is B − → τ − ν with uncertainties around 20% on the branching ratio from the two experiments Belle and BaBar [54,55] and these results show a tension a bit below the level of two combined standard deviations.The averages of lattice results for the decay constants have an uncertainty below one percent (for N f = 2 + 1 + 1) [39], hence the determinations of CKM matrix elements from leptonic decays are currently dominated by the experimental error.In that sense B-mesons decay constants have become a benchmark computation for new methods dealing with heavy quarks on the lattice.A much more precise exclusive estimate of |V ub | is extracted from the semileptonic channel B → πℓν, whose relevant form factors are discussed in the next Section.There theoretical and experimental errors are of comparable size.QED and radiative corrections to leptonic decays of B mesons are of interest as they may be enhanced.Refs.[56,57] discuss terms of O(m b /Λ) and logarithmic enhancements, for B (s) → µ + µ − .In general one expects large collinear logarithms of the form log (m b /m ℓ ) because of the different scales involved.It would obviously be very important to confirm such effects in a fully nonperturbative lattice computation.The situation is not as advanced as for pions and kaons where the strategy put forward in Ref. [58] could be used.This method relies on the use of a pointlike approximation in part of the computation (the real photon emission), which is equivalent to neglecting structure dependent contributions.Such contributions are large for heavy mesons as shown in Refs.[59][60][61], where radiative decays of D mesons have been studied on the lattice.The effect should be even more pronounced for B mesons in which case one expects a faster deterioration of the noise to signal ratio for the relevant correlation functions [61] in Euclidean time.

Exclusive semileptonic decays
Differential decay rates for semileptonic decays of a heavy-hadron H into a lighter hadron h are commonly parameterised as a sum of products of non-perturbative form factors f X and known kinematic factors where q µ = p µ H −p µ h is the momentum transfer onto the lepton pair.Depending on whether the process is tree-level or loop-induced, whether the final state is a pseudoscalar (PS) or a vector (V) state, and whether it is a baryonic or a mesonic decay the number of form factors X that are summed over differs.For the simplest case of a mesonic tree-level PS → PS transition (e.g.B → πℓν) there are two form factors (f + , f 0 ); for a loop-level induced PS → PS transition (e.g.B → Kℓ + ℓ − ) there are three form factors (f + , f 0 , f T ).For treelevel induced PS → V decays (e.g.B → D * ℓν) there are four (V, A 0 , A 1 , A 2 )3 and for loop-level induced PS → V (e.g.B → K * ℓν) there are seven (V, A 0 , A 1 , A 2 , T 1 , T 2 , T 3 ).Depending on the decay, there are also kinematic constraints, relating form factors at particular values of q 2 to each other.

Challenges for form factors
The data analysis relating simulated correlation functions to physical form factors is lengthy and challenging.Whilst it is in principle understood how to perform the simulations required for the prediction of semileptonic form factors, it is costly to generate data sets which enable data-driven control (guided by theory) over the required extrapolations listed in the following: Excited states: On the lattice, form factors can be related to Euclidean matrix elements ⟨h| J |H⟩ of some current J between the initial H and final state h and are extracted as the ground state matrix-element of a three-point function . Accurately isolating the correct matrix element is a trade-off between sufficiently large separations ∆T ≫ t ≫ 0 allowing for the excited state contributions to be exponentially suppressed, and small enough values of ∆T to maintain good control over statistical uncertainties.Based on chiral perturbation theory (χPT), Ref. [62] recently advocated that the expected first excited state energy to this type of three-point functions is expected to be of the size m B + m π and might therefore potentially be hard to disentangle from the desired signal.As an example, JLQCD in their computation of form factors for B → D * semileptonic transitions [63], which we discuss in detail later, shows results for the three-point functions for different source-sink separations, concluding that results stabilise for a separation of about 1.5 fm.(Heavy-quark)-chiral-continuum extrapolation: Typical calculations of these form factors either take place at the physical b-quark mass using an effective action inspired discretisation, or at lighter-than-physical heavy-quark masses.All of the results currently available in the literature utilise ensembles with heavier-than-physical light-quark masses, so that a chiral extrapolation is required in addition to the continuum limit.These extrapolations are typically guided by heavy meson χPT (HMχPT) and therefore rely on inputs stemming from the world at physical quark masses and on low energy constants (LECs), which have to be determined.This makes controlled predictions very challenging when simulations take place at heavier-than-physical light-quark masses and lighter-than-physical heavy-quark masses.Kinematic coverage: Accessible Fourier momenta are of the form a⃗ p = ⃗ n2πa/L where ⃗ n is a vector of integers and |a⃗ p| is required to remain small in order to control discretisation effects.Due to the heavy b-quark mass, it is not typically possible to cover the kinematically allowed range q 2 ∈ [0, (m H − m h ) 2 ≡ q 2 max ] whilst maintaining control over discretisation effects.As a consequence, current calculations only cover a portion [q 2 min,dat , q 2 max ] of the kinematic range at physical kinematics. 4This necessitates extrapolations of the lattice form factors over the full kinematical range.

z-expansions and unitarity constraints
The extrapolations over the full kinematic range are typically carried out as a model-independent z-expansion, a conformal mapping based on Ref. [64] (or variants thereof [65,66]).After removing terms corresponding to sub-threshold poles (P X ) and an appropriate normalisation (ϕ X (q 2 )), the form factor is expressed as a polynomial in the conformal variable z, Since the conformal variable satisfies |z| < 1, this is a convergent series that, for a given truncation, can then be fitted to available lattice data.Based on dispersion relations, a unitarity bound |a n | 2 ≤ 1 can be derived. 5Whether the unitarity bound is satisfied can either be checked a posteriori or the constraint can be imposed directly in the fit, either via the dispersive matrix method [69] or via a recently proposed Bayesian Inference framework [68].The latter allows to trade truncation effects for statistical noise encoding our ignorance of the neglected higher order terms, which are however regulated by the unitarity constraint.In general, multiple decay channels with the same quark-level transition can be incorporated into the same unitarity constraint, strengthening the constraints on a n coefficients.
As an example, since the current in B → πℓν and in B s → Kℓν transitions is the same, the unitarity constraint involves a combination of the corresponding form factors, which should therefore be fitted simultaneously.The unitarity bound will be closer to saturation (compared to imposing the constraint on the individual decay channels), resulting in less freedom and hence tighter constraints on the higher order z-expansion coefficients.This can be systematically improved by including further decays of the same current (e.g.Λ b → pℓν, B s → K * ℓν, ...).

Analysis strategies
Most collaborations perform their analyses in a multi-step procedure, consisting of 1.The extraction of the desired masses and hadronic matrix elements from the correlation function data.2. The extrapolation to zero lattice spacing and physical quark masses and a continuous description of the momentum transfer in the kinematic range covered by the data.The fit is then evaluated at a number of reference q 2values and a full error budget including the correlations between the form factors at the reference values is assembled.3.An extrapolation over the full kinematic range, typically using a model-independent zextrapolation.
Some collaborations have advocated combining the second and the third step into the modified z-expansion [70] in which the coefficients of the z-expansion are changed into functions of the lattice spacing and the quark masses.However, since this is no longer based on an underlying effective field theory, this approach might introduce model dependences, which are hard to quantify with current data.
In the following, we restrict the discussion to form factor calculations with results of multiple collaborations available and cases where tensions have been observed.In particular we will address B (s) → PS and B (s) → V tree-level decays that can be used for the determination of CKM matrix elements and which tend to be most mature.
All other results utilise effective action approaches for the b quark.Figure 3 summarises the status for B → π, in the range where lattice data typically exist.One can clearly see a tension near q 2 max (the vertical black dashed line) for the f 0 form factor.This tension is more severe than is apparent from a plot, since the reference values underlying the individual datasets are highly correlated.When attempting to jointly fit synthetic data points for these results, the FLAG collaboration reports χ 2 /dof = 43.6/12[39], necessitating an error inflation by the PDG scale factor χ 2 /dof.The situation is similar for B s → Kℓν as is demonstrated in Ref. [68] 7 , where the world data cannot be jointly fitted with an acceptable pvalue and χ 2 /dof = 3.89 is the best value that can be achieved.A possible explanation for some of the discrepancies between results, stemming from the assignment of external parameters to the pole structure imposed in the HMχPT has been 6 The result by HPQCD [75] uses gauge field configurations which are also part of FNAL/MILC15 calculation but span a smaller range in the key parameters such as the lattice spacing and pion masses. 7The FLAG report does not yet include the recent RBC/UKQCD23 [76] result.
put forward in Ref. [68], but further independent computations are needed to confirm this.
For the b → c transitions B (s) → D (s) there are fewer independent results.B → D form factors have been published by HPQCD [79] and FNAL/MILC [80].The calculations based on two and four lattice spacings respectively are compatible.The FNAL/MILC result dominates the combined fit, but due to overlapping gauge field configurations (the N f = 2 + 1 asqtad ensembles) between the two computations these results are expected to be statistically correlated.
Two results by the HPQCD collaboration [81,82] exist for B s → D s form factors.The former based on the N f = 2 + 1 asqtad ensembles with NRQCD b-quark and HISQ c-quark, the latter based on the N f = 2 + 1 + 1 HISQ ensembles with HISQ for all quark flavours.Whilst the result covers the entire kinematic range at unphysical kinematics, this is not the case at the physical values of the b-quark mass.To predict form factors at physical masses, an extrapolation of the form factors in the heavy-quark mass is required.This is complicated since the momentum transfer q 2 is also a function of the heavy-quark mass.Ref. [82] uses a modified z-expansion to simultaneously perform the extrapolations to physical light and heavy-quark masses, zero lattice spacing and a continuous description of the momentum transfer q 2 in the physical range.
Additional independent results for all of the above would be very desirable in order to address tensions (where present) and to test the choices that have been made in the literature.Computations are ongoing by several collaborations, in particular by the RBC/UKQCD [83,84], the JLQCD collaboration [85] and the FNAL/MILC collaboration [86].
In recent years the process B → D * ℓν (as well as B → Dℓν covered in the previous section) has received much attention, due to long-standing tensions between experiment and theory concerning these decays [87].In 2021 FNAL/MILC [88] produced the first comprehensive calculation of the four form factors involved away from the zero-recoil point.In 2023 two additional results appeared (HPQCD [89] and JLQCD [63]), which, at the time of writing this report, are only available as pre-prints.Due to the interest in these very recent results, we compare the three calculations in detail in the following Lattice set-ups: The FNAL/MILC [88] work uses 15 N f = 2 + 1 ensembles with the asqtad fermion action with five lattice spacings a ∈ [0.045, 0.15] fm and a range of (root-mean-square) pion masses 8 down to m RM S π ∼ 250 MeV.Both b and c quarks are simulated with the Fermilab action.The HPQCD [89] computation uses five N f = 2+1+1 ensembles with the HISQ action for all sea and valence quarks including three lattice spacings a ∈ [0.044, 0.090] and two ensembles with physical Goldstone-pion masses The JLQCD [63] calculation uses nine N f = 2 + 1 ensembles with the domain wall action for all quarks and three lattice spacings a ∈ [0.04, 0.08] fm.Simulations take place in the range 226 MeV ≲ m π ≲ 500 MeV.HPQCD and JLQCD simulate a range of heavyquark masses below the physical b-quark mass with the constraint am q ≤ 0.7 and 0.8, respectively and therefore require an extrapolation in the heavy-quark mass.
Heavy-to-heavy transitions are often described as a function of the kinematic variable w = v B • v D * .Since all the mentioned works take place in the B-meson rest-frame this reduces to w = E D * /M D * .The range of this parameter approximately corresponding to the semileptonic range 8 For staggered quarks several light states have the quantum numbers of the pion.Their masses differ by taste-breaking O(a 2 ) effects.The lightest state is called Goldstone-pion.
is w ∈ [1, 1.5].The FNAL/MILC, and JLQCD results cover the range from w = 1 to w ∼ 1.175 and 1.1, respectively on all their ensembles.The HPQCD result covers the range from w = 1 to 1.05, 1.20, and 1.39 on their a ≈ 0.09, 0.06, and 0.044 fm ensembles, respectively.Lattice to continuum: Each of the three works perform a fit to simultaneously extrapolate to the continuum and to physical quark masses and to obtain a continuous description of the form factors as a function of w.
All three collaborations use χPT (or staggered variants of this) to extrapolate to the physical pion mass and describe the kinematic behaviour as an expansion in (w − 1) k up to quadratic (JLQCD, FNAL/MILC) or cubic (HPQCD) order.
FNAL/MILC includes discretisation effects of order α s aΛ, (aΛ) 2 and (aΛ) 3 .JLQCD accounts for terms (aΛ) 2 and (am q ) 2 since due to the use of domain wall fermions discretisation effects from odd powers of the lattice spacing are absent (up to terms proportional to the residual mass which are negligible).JLQCD's heavy quark extrapolation is performed by first dividing out the matching factor between HQET and QCD and then parameterising this result in powers of 1/m h up to order 1.HPQCD parameterises discretisation effects and the heavy-quark mass dependence as products of (am c ) 2i (am h ) 2j [(Λ/M Hs ) k −(Λ/M Bs ) k ] for i, j, k = 0, 1, 2, 3.One such term is present for each order of (w − 1) n for n = 0, 1, 2, 3.
Per form factor, the resulting fits typically have 8 parameters for JLQCD, 10 parameters for FNAL/MILC and more than 250 parameters for HPQCD.JLQCD performs the fits as frequentist χ 2 -minimisations, whilst FNAL/MILC and HPQCD use a Bayesian framework with Gaussian priors for their fit parameters.
One caveat concerning the HPQCD computation [89] is that for each form factor the B → D * and B s → D * s data are fitted jointly, imposing that all the fit parameters are the same.The only freedom for the fit to distinguish between B → D * and B s → D * s is by some of these parameters multiplying an SU(3) breaking term.Since Ref. [89] does not present fits to the individual channels it is not possible to quantify the impact of this choice.Ref. [89]   of only B s → D * s form factors [90] in the helicity basis, which will be described below.z-expansions: Fermilab/MILC and JLQCD perform a two-step analysis, evaluating their chiralcontinuum-heavy-quark extrapolation at three values of w in the range of data where the simulations took place (w They each assemble a full error budget for all four form factors at these data points and quantify all correlations.These reference values serve as inputs for a subsequent z-expansion, which extrapolates the form factors over the full kinematic range.
HPQCD provides results for a z expansion, but does not clarify what input data was used for this z-expansion.It is unclear whether the fit was performed to reference w-values, and if so, how they were chosen, how the full error budget for these reference values is computed and what the correlations between these values are.
Since all three collaborations use the same conventions for the z-expansion, the results can be directly compared and the coefficients from the three collaborations are listed in Table 1 and displayed in Figure 4. 9 There are clearly several tensions between the results for the lower order coefficients which will need to be understood.B s → D * s form factor results by HPQCD: We now turn our attention to the comparison of the results for the B s → D * s form factors from the two HPQCD calculations from 2023 [89] and 2021 [90].
The only difference in the ensembles underlying this dataset is one additional ensemble in [89] at the intermediate lattice spacing with physical pion mass providing additional data in the range w ∈ [1.0, 1.2].Since there are no valence light quarks in the decay B s → D * s this is not expected to significantly alter the result.Furthermore, stemming from largely the same ensembles, the two results are highly statistically correlated, so that any discrepancies between the results are expected to stem from systematic uncertainties.
The form factors f and g largely remain within one standard deviation, but with a reduction of the quoted uncertainties by a factor of approximately 3.5.A similar reduction of uncertainties is quoted for the form factor F 1 .Whilst this form factor agrees with the previous result near w ∼ 1, the slope with w is very different and in the intermediate w region, there is an approximately 2σ tension between the two results.For the final form factor F 2 the error reduction is less significant, but the slope near w ≈ 1 is different between the two calculations.
The authors of Ref. [89] explain the observed differences with the change in analysis strategy from an expansion in powers of the conformal variable z to an expansion in the variable (w − 1) and by the joint fitting of B → D * and B s → D * s data.Since the joint fits in Ref. [89] are dominated by the statistically more precise B s → D * s data, the latter seems unlikely to explain the differences.According to Ref. [89], the main contribution to the total uncertainty is statistical, with systematic uncertainties accounting for at most 15%, 10%, 25%, and 20% of the variance of f , g, F 1 , and F 2 , respectively.Further investigations are required to understand the differences between Refs.[89] and [90].It would be interesting if correlated differences between the two results could be produced.

Desirable benchmark quantities, comparisons and checks
The extrapolations required to predict form factors that can be related to the physical world from the underlying lattice data points are more complex than those for quantities that do not depend on a kinematic variable.This is mainly due to the fact that the kinematically allowed range (q 2 max or equivalently w max ) changes as a function of the mass of the initial and the final states.This is further aggravated by the increase of statistical noise as the heavy-quark mass becomes heavier and as larger momenta are induced.In particular when simulating multiple heavy-quark masses, this often leads to a comparably small number of precise data points (at small m h , small ⃗ p) and a lot of data points with sometimes orders of magnitude larger uncertainties (large m h , large ⃗ p).These in turn need to be described by (often complicated) multidimensional fits.As a result a small portion of data points drive the fit, typically those furthest away from the desired physical parameters.The fit functions that are employed often have many parameters.In order to assess the weight of the different data in the fit and to determine the relevant fit parameters, it could help to i) fit only the less precise data, to see what their effect is, and/or ii) start by only including the most precise data (and therefore using simple fit ansätze) and then adding less precise data until results stabilise.This would provide an assessment of which portions of the covered parameter space constrain the fit.On a related note we remark that whilst datapoints with very large uncertainties (relative to the majority of the data) do not significantly contribute to the χ 2 , they alter the interpretation of how good a particular fit result is, by reducing the χ 2 /dof without adding much information.In the extreme case where the relative uncertainty of data points differs by orders of magnitude, it might be worth to develop a notion of effective degrees of freedom.
Given the tensions in several form factor computations, it would be desirable to devise some easier-to-compute benchmark quantities that all collaborations could provide, somewhat analogous to the window quantities in the hadronic vacuum polarisation contribution to the anomalous magnetic moment of the muon [91].These should be designed to allow for comparisons between different computations with reduced sources of systematic uncertainties.Particularly well suited are quantities which disentangle the effects of the different dimensions of the fit, such as the chiral, continuum, kinematic and heavy quark extrapolations.We conclude this section with the incomplete list of suggested benchmark quantities below: • A full error budget for the form factors f 0 (q 2 max ) for PS → PS transitions and f (w = 1) for PS → V transitions solely based on the zero momentum data points.Since all hadrons are at rest, these can be obtained from the most precise data points and no interpolation in the final state energy is required.
• Generalising the previous point: the form factor at kinematic reference points.Whilst several collaborations tend to provide these as input for subsequent kinematic extrapolations, they are typically obtained from a joint fit of all simulated momentum values.Performing such checks only from data in the vicinity of the kinematic reference point would allow more direct comparisons which do not rely on the broader choice of kinematic description of the data.For example, differences between choosing a modified z expansion, an expansion in (w − 1) or an expansion in the final state energy could be assessed.It would also shed light on how much of the information content stems from simulated data near the kinematic reference point as opposed to potentially more precise data which is kinematically far away.• For simulations which take place at lighterthan-physical heavy-quark masses, it would be valuable to perform the chiral-continuumkinematic extrapolation at fixed (unphysical) heavy-quark mass.If data was generated with this in mind, one could simulate data at, for example, half the b-quark mass and make predictions of the form factors at this choice of kinematics.This would remove the heavy-quark extrapolation dimension from the fit and hence disentangle the two most complicated extrapolations: the heavy-quark mass extrapolation and the continuum-limit.Furthermore, since the simulated heavy-quark mass in lattice units is smaller, better control over the continuum limit would be expected.These results could be compared between different collaborations.• Assessment of the contamination of ground state matrix elements and masses by excited states is a recurring feature of lattice QCD calculations.The excited state energies obtained from correlation function fits are physical quantities and should (up to discretisation effects) agree for different fermion formulations at equivalent quark masses.If this data was included in publications, it would help to build confidence in the absence of excited state contamination or to understand their nature.• In connection with the discussion on the varying size of statistical uncertainties, it would be valuable if the statistical correlations between all data points stemming from the same ensemble were given.Since this is a vital ingredient of the fit extrapolating to physical parameters, this information is required in order to reproduce results.Furthermore, these correlations should be universal (up to discretisation effects), a fact that could be checked if these results were provided by all collaborations.

Neutral meson mixing
Neutral mesons such as the B 0 q (q = d, s) mix with their antiparticles B0 q through box diagrams such as the one displayed in Figure 5. Due to this mixing, the mass and flavour eigenstates of the B q − Bq system do not coincide, resulting in experimentally measurable mass and width differences.Since the box diagrams are top-quark and therefore short distance dominated, the relevant non-perturbative matrix elements are calculable on the lattice.An operator product expansion of the box diagrams above yields 5 independent parity-even dimension-6 operators, whose matrix elements can be computed in LQCD.
For historical reasons, the matrix elements under consideration are typically cast into the form of bag parameters B (i) Bq where i = 1, ..., 5 and q = d, s, quantifying the departure of the matrix element from the vacuum saturation approximation (VSA).These are defined as where the η q ensure that the expression is unity in the VSA.For phenomenological applications, typically the quantities f Bq B (i) Bq are of relevance, but depending on the application, the bag parameters or the matrix elements themselves are also of interest.
In the SM, the experimentally measured mass difference ∆m q is related to Bq by known multiplicative factors and the product of CKM matrix elements V tq V * tb (q = d, s), and hence precise knowledge of the non-perturbative inputs enables the determination of CKM matrix element containing the top quark and thereby contribute to tests of CKM unitarity.Since the nonperturbative matrix elements are independent of the UV-properties of the theory under consideration they can also be used to probe BSM models.Finally, several theoretical uncertainties cancel in the SU(3)-breaking ratios B (1)

Bs /B
(1) so that non-perturbative determinations of this quantity provide more stringent bounds on |V td /V ts |.
Full results for all 5 operators (and in both systems) are available from ETM [92] (N f = 2), FNAL/MILC [93] (N f = 2 + 1) and HPQCD [94] Table 2 Some key properties of the results discussed in the text.The abbreviations for the various fermion actions stand for Wilson twisted mass (Wtm), a 2 -tadpole improved staggered quarks (asqtad), highly improved staggered quarks (HISQ), domain wall fermions (DWF) and Osterwalder-Seiler (OS).† The root-mean-square pion mass is listed.The corresponding goldstone pion masses ranges are [177,555] Bs /B (1) B d and ξ from N f = 2 + 1 QCD [95].One difference between these computation is which quantities are directly computed and which ones are inferred using available results for f Bq from the literature.For maximally exploitable phenomenological applications of these results, it would be desirable to have access to Bq and f Bq B (i) Bq from the same computation, in order to be able to quantify all correlations between these observables.
Table 2 summarises some of the key properties of these results, which are highly complementary: In addition to differences in the choice of sea, valence light, and valence heavy quark actions, the range of lattice spacings and simulated pion masses differs significantly.When chiral symmetry is maintained, the five bag parameters mix in a continuum-like fashion under renormalisation, so the renormalisation pattern is block diagonal.In the formulations used by the ETM and the RBC/UKQCD collaborations, this property is preserved in the lattice simulation and hence simplifies the renormalisation procedure.The treatment of the heavy-quark also differs, with FNAL/MILC and HPQCD employing an effective action approach and ETM and RBC/UKQCD using a fully relativistic set-up.
Observables are often quoted in the Renormalisation Group Invariant (RGI) scheme.Results for the phenomenologically relevant quantities B(1) Bq , f Bq B(1) Bq for q = d, s and their SU(3)-breaking 0.9 1.0 1.1  ratios are compared in Fig. 6.In addition to the works already described, this also includes HPQCD09 [96] and RBC/UKQCD14 [97].Results for the remaining 4 bag parameters in the MS scheme at µ = m b are shown in Figure 7.The red left facing triangles representing the ETM N f = 2 result are shown with open symbols, since the calculation does not quantify any uncertainty related to the missing sea strange quark, which in many observables has been shown to be significant.Since missing sea-charm effects are assumed to be negligible at this level of precision, N f = 2 + 1 and N f = 2 + 1 + 1 results can be directly compared and the magenta bands correspond to weighted averages of these.In cases where χ 2 /dof ≥ 1.25, the dashed magenta lines indicate the uncertainties obtained after the application of the PDG scale factor χ 2 /dof.The results from FNAL/MILC and HPQCD dominate the average, but show some tension.The χ 2 /dof for the weighted average of just these two results is 1.2 and 2.4 for B(1) , respectively and larger than 3 for both quantities in the B s case.It is noteworthy, that these tensions disappear when considering SU(3)-breaking ratios, since the effect is correlated between the B d and the B s system.The weighted averages of the observables shown in the right hand column of Fig. 6 have a precision of 2.1%, 2.0% and 0.8%, respectively.We note however, that if the average was only taken between the FNAL/MILC and the HPQCD result, the uncertainties for the first two quantities would be 3.5% and 3.3% due to the required rescaling of the uncertainties.The uncertainty on the CKM matrix elements |V td | and |V ts | and their ratio , which can be extracted by combining these observables with the experimental measurements of ∆m d and ∆m s , is clearly limited by the theoretical uncertainty, due to the per mille-level precision of the experimental inputs.This necessitates further theory improvements.
Turning our attention to the remaining bag parameters (see Figure 7) there is agreement between most results, but for B Bq we notice a tension between the FNAL/MILC and the HPQCD results.The former is very similar to what is observed in neutral kaon mixing [98] and might be related to the choice of RI/MOM instead of RI/SMOM.Clearly further independent calculations are desirable to address and resolve these tensions.
A joined effort between the RBC/UKQCD and the JLQCD collaboration [99,100]  Ref. [95] to predictions of the full set of operators in a fully relativistic set-up.This is achieved by complementing the existing data set described in the last column of Table 2 by JLQCD domain wall fermion ensembles, albeit with slightly different action parameters, with three lattice spacings in the range a ∈ [0.04, 0.08] fm.This allows simulations up to close to the physical b-quark mass, controlling the remaining extrapolation to its physical value.The chiral extrapolation is controlled by the inclusion of the physical pion mass ensembles of the RBC/UKQCD collaboration.The use of chirally symmetric domain wall fermions throughout, guarantees a fully non-perturbative renormalisation procedure with continuum-like mixing pattern, addressing one of the main systematic uncertainties present in current calculations.

Recent developments
Many (actually all) of the PS to V transitions described in the Section on semileptonic form factors are decays to multiple (two) hadrons in the final state.The D * for example decays strongly into Dπ.Such effects are typically included in the form factor chiral extrapolations through χPT or HMχPT [101].A formalism to treat such processes on the lattice however has been introduced several years ago by Lellouch and Lüscher in Ref. [102] and extended for the cases discussed here in Ref. [103].First numerical applications [104] have appeared only recently for the B → ρ(ππ)ℓν process, which provides an alternative exclusive channel for the extraction of V ub .Further results, also for other channels, can reasonably be expected in the near future.The formalism relies on the relation between two-particle energy levels in a finite volume and the infinite-volume scattering lengths [105].Since energy levels are directly continued to Euclidean time this allows to compute properties of scattering processes on the lattice.In a second step a perturbation is introduced inducing the transition between a single particle and two particle states (K → ππ is the process considered in Ref. [102]).By computing (to leading order) the effect of such a perturbation on the energy levels and the scattering, one obtains a proportionality relation between the finite volume matrix element and the (say, K → ππ) decay amplitude in infinite volume.The proportionality factor, known as "Lellouch-Lüscher" factor is a function of the momenta of the particles and the derivative of the scattering phases with respect to them.In Ref. [103] a further step is taken by extending the approach to the case of currents (the perturbations above) inserting energy, momentum and angular momentum for systems with an arbitrary number of mixed two particle states.In this case the application in mind is the B → K * (Kπ)ℓ + ℓ − transition.
Another topic in which considerable progress has been made, is the study of inclusive decays on the lattice, with the aim to shed some light on the tensions between inclusive and exclusive determinations of CKM matrix elements such as V cb and V ub .In a series of papers [106][107][108][109] an approach has been devised, taking the process B (s) → X c ℓν as prototype.The central quantity is the hadronic tensor where J µ is the weak current inducing the b → c transition, q is the lepton-pair momentum and the sum over the charmed final state X includes a spatial-momentum integral (in ⃗ p X ).At fixed p B , for example by choosing the rest-frame of the B meson, the tensor above is a function of the spatial components ⃗ q and ω = M B − q 0 .What can be computed on the lattice is essentially the Laplace Transform of such tensor at fixed ⃗ q C µν (⃗ q, t) = ∞ 0 dωW µν (⃗ q, ω)e −ωt . (15) Since C µν in finite volume is given by a sum of exponentially falling (in time) functions, trying to invert the relation above for arbitrary values of ω is an ill-posed problem.For ω smaller than the energy of the lowest state X the relation however can be inverted.In Ref. [106] the tensor W µν (⃗ q, ω) is computed on the lattice for that particular unphysical kinematic choice (ω < E min X ), where the final hadronic state can not go on shell and is then related to derivatives of the tensor in the physical region through a dispersion integral.The approach is very similar to the derivation of the moment sum rules in the case of Deep Inelastic Scattering.
Improving on this first approach a completely new method has been introduced in Ref. [107], based on the observation that in order to compute the inclusive decay rate what is needed is not the hadronic tensor but rather a smeared version of it with functions resulting from the leptonic tensor.The building blocks are the quantities X(l) (|⃗ q| 2 ) = ∞ 0 dωW µν (⃗ q, ω)K µν,(l) (⃗ q, ω) , (16) where the functions K µν,(l) (⃗ q, ω) are known and defined in Ref. [107].If one can find good polynomial approximations of the K-functions in powers of e −aω , then the X(l) can be obtained by taking suitable combinations of the correlator C µν at different times.Obviously the larger the time one can compute the correlation function with controlled errors, the higher the degree one can consider the polynomial approximation and the better the accuracy in X(l) (|⃗ q| 2 ), and eventually in the inclusive decay rate.In Ref. [107] Chebyshev polynomials have been studied while in Refs.[108,109], in an attempt to optimise the approximation balancing truncation errors and statistical noise from the correlator C µν , also the Backus-Gilbert method has been considered.All such studies are still exploratory, however the results are very encouraging.

Conclusions and outlook
The combination of Lattice QCD and B-physics constitutes an exciting and active field.In this review, we have commented on the challenges faced when making predictions for B-physics observables from lattice QCD calculation.We have reviewed the recent literature for the determination of the b-quark mass, leptonic decays constants, several semileptonic decay channels and neutral meson mixing parameters.To address some of the observed tensions, we have put forward a number of suggestions for benchmark quantities, which will help to understand systematic effects associated with complicated multidimensional fits.
Some quantities (b-quark mass, decay constants) are in a mature state, with many results using different methodologies agreeing, producing comparable uncertainties, and the corresponding tests of the SM are limited by the experimental precision.The more complicated calculation of semileptonic form factors has made tremendous progress, but due to to the wealth of decay channels and the increased computational and dataanalysis complexity, fewer independent results exist for the individual observables.Several of the results in the literature have not yet reached a community consensus and currently display some level of tension.However, there are major ongoing efforts by several groups to consolidate these calculations.We hope that the benchmark quantities we suggest in Sec.4.6 will be adopted by the community and will aid in shedding light on current tensions between different lattice computations.The situation for neutral meson mixing parameters is similar with only a small number of results currently available.
Finally, we comment on recent progress in calculations of semileptonic decays into two final state hadrons and inclusive semileptonic decays.Whilst these calculations are currently exploratory they are very promising and we expect they will mature to further ab-initio tests of the SM.
It is worth emphasising the increase in complexity in the quantities presented here.Masses and decay constants can be extracted from twopoint functions, form factors require three-point functions while for two-hadron decays and inclusive transition rates four-point functions are needed.The statistical noise in n-point functions unavoidably grows with n and in addition the combinatorics in the number of Wick contractions becomes more complex.It is a remarkable achievement of the advances in numerical methods and computer power combined with smart field theoretical ideas that has led us to tackle new challenges in B-physics on the lattice, which seemed un-treatable only a few years ago.

Fig. 4
Fig. 4 Comparison of the z-expansion coefficients for the B → D * form the factors g, f , F 1 and F 2 obtained by FNAL/MILC (blue circles), HPQCD (red triangles), and JLQCD (green squares).

Fig. 6
Fig. 6 Summary of the first bag parameter and f √ B in the RGI scheme.The magenta bands correspond to the values of a weighted average of N f > 2 results with its value given in the top right corner of each panel.

10 4 ) 5 )
The results from ETM for B (Bq and B (Bq have been recast into the same form as the results presented by HPQCD and FNAL/MILC, to guarantee that they are unity in the VSA.The required input for M Bs , M B d , m b (m b ), m d (m b ) and ms(m b ) were taken from Ref. [93].
Bq we notice a slight tension between the ETM and the other two results, whilst for B(3)

Fig. 7
Fig.7Using results from Refs.[92],[93] and[94] Results are quoted in MS(m b ) and using the BBGLN scheme.Magenta bands are the weighted averages between N f = 2 + 1 and N f = 2 + 1 + 1.In cases where χ 2 /dof > 1.25 the dashed magenta lines indicate the uncertainty including the PDG scale factor of χ 2 /dof.

Table 1 z
-expansion coefficients quoted by the three collaborations.