Inclusive semi-leptonic decays from lattice QCD

We develop a method to compute inclusive semi-leptonic decay rate of hadrons fully non-perturbatively using lattice QCD simulations. The sum over all possible final states is achieved by a calculation of the forward-scattering matrix elements on the lattice, and the phase-space integral is evaluated using their dependence on the time separation between two inserted currents. We perform a pilot lattice computation for the B_s ->X_c l nu decay with an unphysical bottom quark mass and compare the results with the corresponding OPE calculation. The method to treat the inclusive processes on the lattice can be applied to other processes, such as the lepton-nucleon inelastic scattering.


Shoji Hashimoto
Theory Center, Institute of Particle and Nuclear Studies, High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801, Japan and School of High Energy Accelerator Science, The Graduate University for Advanced Studies (SOKENDAI), Tsukuba 305-0801, Japan (Dated: May 29, 2020) We develop a method to compute inclusive semi-leptonic decay rate of hadrons fully nonperturbatively using lattice QCD simulations.The sum over all possible final states is achieved by a calculation of the forward-scattering matrix elements on the lattice, and the phase-space integral is evaluated using their dependence on the time separation between two inserted currents.We perform a pilot lattice computation for the Bs → Xc ν decay with an unphysical bottom quark mass and compare the results with the corresponding OPE calculation.The method to treat the inclusive processes on the lattice can be applied to other processes, such as the lepton-nucleon inelastic scattering.
Quark-hadron duality plays a key role in perturbative Quantum Chromodynamics (QCD) calculations of physical processes.It states that hadronic processes can be calculated taking quarks and gluons as final states, even though the actually observed final states are composed of hadrons.In order that the duality is satisfied, the processes must be summed or smeared over all possible hadronic final states in some kinematical range [1], such as a region of invariant mass squared, but it is not a priori known how large the smearing should be.A systematic approach to duality is based on the Operator Product Expansion (OPE) [2,3], which is constructed in the Euclidean domain and analytically continued to the Minkowski domain.In the context of heavy quark decays the OPE is an expansion in inverse powers of the heavy quark mass, or more precisely of the energy release, and the analytic continuation entails an inevitable violation of duality.While there are indications that duality violation plays a minor role in the analysis of inclusive semi-leptonic B meson decays to determine the Cabibbo-Kobayashi-Maskawa (CKM) matrix element |V cb | [4,5], full control of the systematic error can only be achieved by non-perturbative methods.The current tension between the inclusive and exclusive determinations of |V cb | [6,7] makes any contribution in this direction timely.
Lattice QCD simulation provides a means of nonperturbative QCD computation for various hadronic processes including heavy quark decays.It has been successfully applied to the calculation of exclusive decay form factors, which are essential for a precise determination of the CKM elements, e.g.K → π ν, B → D ( * ) ν, etc. (see [8] for their world averages), while the study of inclusive processes is scarce, except for recent attempts to formulate methods to introduce an analytic continuation [9] or a smearing [10].The inclusive processes are, on the other hand, difficult to treat on the lattice because they consist of many physical states often including multiple hadrons.To identify each amplitude and to sum over the phase space is nearly impossible due to the number of states involved.One may instead use analyticity and the optical theorem to relate the total rate to another quantity that is calculable on the lattice.This approach has been followed in simple cases, such as the e + e − → q q processes [11][12][13][14] and τ -lepton decays [15,16], while the application to B meson semi-leptonic decays is much more complicated [9].
In this work we develop a novel and general method to compute the inclusive semi-leptonic decay rate on the lattice.The method is based on a technique to calculate smeared spectral density of hadron correlators [17] (see also [10,18] for a slightly different strategy).The extraction of the spectral density ρ(ω) of hadronic correlation functions remains intractable, but once ρ(ω) is smeared over some energy range, one can construct a good approximation using the correlation functions calculated on the lattice.In semi-leptonic decays of hadrons, the phase-space integral plays the role of this smearing.The method is systematically improvable as more computational resources are made available.
In this paper we use the inclusive semi-leptonic decays Bs → X c ν to demonstrate how the method works.Here, X c stands for all possible charmed states which may occur with the quark-level decay process b → c ν.After describing the kinematics of the decay and the method to calculate the inclusive decay rate, we present a pilot lattice study.
For the analysis of the Bs → X c ν decay, we assign a momentum p µ to the initial Bs meson, and momenta p µ and p µ ν to the leptons and ν in the final state, respectively.Thus, the hadronic state X c has momentum r µ = (p−q) µ with q µ = (p +p ν ) µ .The differential decay rate is written as [19,20]

arXiv:2005.13730v1 [hep-lat] 28 May 2020
where G F is the Fermi constant.The momentum transfer q µ and the lepton energy E are evaluated in the rest frame of the initial Bs meson.The leptonic tensor L µν is explicitly written as L µν = p µ p ν ν − p • p ν g µν + p ν p µ ν − i µανβ p ,α p ν,β for massless neutrinos.The hadronic tensor W µν (p, q) is defined through It is summed over all possible final states X c to represent the inclusive decay.The electroweak current relevant for this decay mode is One can perform an integral over the lepton energy E in (1), and the remaining integrals over q 2 and q 0 can be rewritten in terms of ω and q 2 , energy and spatial momentum squared of the final hadrons X c , respectively.Thus, the total decay rate can be calculated as where with Here, we take the momentum q in the k-th direction, while the i-th direction is assumed to be perpendicular to that.The repeated indices in ( 5)-( 7) are not summed.
The integral with respect to ω in (4) represents the sum over states that could appear for a given momentum q.
On the lattice, as a counterpart of the hadronic tensor W µν , one can calculate the forward-scattering matrix elements of the form [9] from four-point functions including the interpolating operators for the Bs meson state | Bs (0) .Now we introduce the transfer matrix on the lattice e − Ĥt to express the time dependence of the matrix element in (8) as where Jν (q) denotes a Fourier transform of the inserted current: Jν (q) = x e iq•x J ν (x).On the other hand, the integral over ω in (4) can be rewritten in the form Here K(ω, q) represents an integral kernel determined by the explicit form of the integrands ( 5)- (7).The ωintegral is implicit on the right hand side; all the intermediate states may exist between the currents.Comparing the right hand side with (9), we find that the integral (10) can be evaluated if the kernel operator is well approximated by a polynomial of the form (11) with some coefficients k j (q), since the matrix elements of the individual term on the right hand side are nothing but C JJ µν (t; q)'s.The best approximation of K( Ĥ, q) can be achieved using the Chebyshev polynomials.We define a state |ψ µ (q) on which the kernel operator is evaluated as Jµ (q)| Bs (0) .A small time evolution e − Ĥt0 with a constant time t 0 is introduced to avoid any potential divergence in ψ µ (q)|ψ ν (q) .We can then construct an approximation as (The dependence on q is omitted for simplicity.)T * j (x) stands for the shifted Chebyshev polynomials, which are derived from the standard Chebyshev polynomials T j (x) as T * j (x) ≡ T j (2x − 1), so that they are defined in the range 0 ≤ x ≤ 1.Their first few terms are and the others can be obtained recursively by according to the general formula of the Chebyshev approximation.The Chebyshev approximation is the best in the sense that its maximum deviation in x ∈ [0, 1] is minimized among polynomials of order N .The integral kernel K(ω, q) is chosen as for l = 0, 1, or 2 corresponding to X (l) , ( 5)-( 7).An approximate Heaviside step function θ σ (x) is introduced to (ω) with the Chebyshev polynomials of e −ω .For each value of the smearing width σ (= 0.2 (top), 0.1 (middle), 0.05 (bottom)), the approximations with the polynomial order N = 5 (dotted), 10 (dot-dashed), 20 (dashed) are plotted as well as the true curve (solid curve).
realize the upper limit of the ω-integral.In order to stabilize the Chebyshev approximation, we smear the step function over a small width σ.For an explicit form, we chose θ σ (x) = 1/(1 + exp(−x/σ)).The extra factor e 2ωt0 in (14) cancels the short time evolution e − Ĥt0 in |ψ µ (q) .Fig. 1 demonstrates how well K (l) σ (ω) is approximated with certain orders of the polynomials, i.e.N = 5, 10 and 20.An example for l = 0 is shown.Here we take three representative values of σ: 0.2, 0.1 and 0.05 in lattice units.The comparison is made for parameters that roughly correspond to our lattice setup: the inverse lattice spacing 1/a 3.61 GeV, am Bs 1.0, t 0 /a = 1.The momentum insertion q is set to zero.The kernel function is well approximated with relatively low orders of the polynomials, such as N = 10, when sufficiently smeared, e.g.σ = 0.2.For smaller σ's, the function exhibits a more rapid change near the threshold ω = 1.0, and one needs higher orders, like N = 20.Eventually we have to take the limit σ → 0, and the error due to finite N has to be estimated.For l = 1 and 2 the polynomial approximations are better than those for l = 0.
We perform a pilot study of the method described above using lattice data computed on an ensemble with 2+1 flavors of Möbius domain-wall fermions (the ensemble "M-ud3-sa" in [21], which has 1/a = 3.610(9) GeV).For the charm and bottom quarks in the valence sector, the same lattice formulation is used.The charm quark mass m c is tuned to its physical value and the D s and D * s meson masses are 1.98 and 2.12 GeV, respectively.The bottom quark mass is taken as 2.44m c , which is substantially smaller than the physical b quark mass.The corresponding B s meson mass is 3.45 GeV.In this setup, the maximum possible spatial momentum in the Bs − m 2 Ds )/2m Bs 1.16 GeV.The lattice volume is L 3 × L t = 48 3 × 96, and we calculate the forward-scattering matrix elements with spatial momenta q of (0,0,0), (0,0,1), (0,0,2) and (0,0,3) in units of 2π/La.The number of lattice configurations averaged is 100, and the measurement is performed with four different source time-slices.
First, we inspect how well the Chebyshev approximation works by comparing the results for X(2) obtained with the polynomial order N = 5, 10, 15 at various values of σ, the width of the smearing.Fig. 2 shows that the dependence on σ is mild and the limit of σ = 0 is already reached at around σ = 0.05.The dependence on N is not significant, which indicates that the approximation is already saturated at N 10.This is crucial because the error of the lattice data is too large to constrain the matrix elements ψ µ |T * j (e − Ĥ )|ψ ν / ψ µ |ψ ν at j 10 or larger.The results for X(0) and X(1) show the similar tendency.We take σ = 0.05 in the following analysis; the results are within statistical error even if we extrapolate to σ = 0.
The lattice results for X = 2 l=0 X(l) are compared with the OPE predictions in Fig. 3 as a function of q 2 .Here, the results for different polarizations, i.e. longitudinal ( : µ, ν = 0 and 3) and perpendicular (⊥: µ, 3. X as a function of q 2 plotted in the physical unit.Longitudinal ( ) and perpendicular (⊥) polarizations are plotted for vector (V V ) and axial-vector (AA) channels.Dotted and dashed curves show the lowest order and O(1/m 2 ) OPE estimates for each channel of corresponding color, respectively.
ν = 1 and 2) directions to q, are separately plotted for vector (V V , squares) and axial-vector (AA, circles) current contributions.The lowest order and O(1/m 2 ) OPE estimates [20] are shown in the same plot.The OPE predictions are sensitive to the heavy quark masses.We take the MS mass for the charm quark, mc (3 GeV) = 1.00 GeV, and the kinetic mass for the fictitious b quark, m kin b (1 GeV) = 2.70(4) GeV, tuned to reproduce the B s meson mass in the simulation using the results of [22].For the OPE matrix elements we employ the results of the semi-leptonic fit of [5], although they refer to a light spectator and to the physical b mass.The dashed lines include O(1/m 2 ) power corrections, which are large and tend to improve the agreement with the lattice data compared to the free quark decay (dotted lines).To obtain the total decay rate, we integrate X q 2 over q 2 as in (3).The vector and axial-vector contributions of different polarizations are added.The integrand is shown in Fig. 4. We fit X(l) / q 2 2−l by a polynomial of q 2 to interpolate the data points.The fit curve (dot-dashed) is terminated at q 2 max .We compare the lattice results with the corresponding OPE prediction (red curve) including O(1/m 3 ) [23] and O(α s ) [24] terms with α s = 0.27.The power corrections are controlled here by powers of the partonic energy m 2 c + q 2 which ranges between 1 and 1.5 GeV, significantly less than that for a physical b.They are singular at the partonic endpoint, where the maximum energy hits the mass-shell of charm quark and the perturbative corrections show an integrable singularity.
Integrating the fit to lattice data we obtain Γ/|V cb | 2 = 4.9(6) × 10 −13 GeV, where only the statistical error is shown.We note that the total decay rate is about five times smaller than that of the physical B s meson, because of the smaller phase space for the artificially small b quark mass.On the OPE side, several higher order corrections are available for the total width, including the complete O(α 2 s ) [25,26] and the O(α s /m 2 ) [27,28] corrections.We implement them in the kinetic scheme using the same inputs as above and obtain Γ/|V cb | 2 = 5.4(8) × 10 −13 GeV.The dominant uncertainty is due to the value of the b quark mass, but missing higher order corrections and uncertainties on the matrix elements would also induce an O(10%) uncertainty.Despite these limitations, the agreement between the lattice and the OPE is remarkable.
An immediate extension of this work is of course the calculation of the inclusive semi-leptonic decay rate of B mesons and b baryons (b → c ν and b → u ν).Moments of kinematical variables, such as the lepton energy moments and hadronic invariant mass moments, can also be calculated by a slight modification of the method.A nu-merical challenge for the lattice calculation is the large recoil momentum up to ∼ 2.3 GeV, which requires fine lattices to keep the discretization effects under control.For b → u transitions, the experimental analysis involves various momentum cuts to veto unnecessary b → c backgrounds.Our method allows to apply arbitrary kinematical cuts, and a fully non-perturbative calculation is possible according to the experimental setup.A comparison to the OPE calculation at or closer to the physical b mass would provide a valuable test of the OPE, including the assumption of quark-hadron duality.It may also be used to determine the hadronic parameters appearing in the heavy quark expansion.The fully non-perturbative lattice calculation can also be applied to D meson decays, for which the energy release is not sufficiently large to yield reliable OPE calculations, and where one could observe the onset of quark-hadron duality.
The possible applications of the framework are not limited to heavy quark decays.Lepton-nucleon ( N ) scattering is another large area of application.Traditionally, it has been analyzed combining perturbation theory and non-perturbative inputs, such as the parton distribution functions (PDFs).Instead, the method described in this work allows to directly compute the cross sections with-out recourse to intermediate quantities like PDFs, and it opens a new strategy to study the inelastic scatterings.Moreover, it will make it possible to perform nonperturbative calculation of low-energy scatterings, which cannot be treated with the presently available techniques.

FIG. 4 .
FIG.4.Integrand of the q 2 -integral plotted in the physical unit.The dot-dashed curve is an interpolation of the lattice data, and the O(1/m 3 , αs) OPE calculation is shown by the red curve.