Inclusive semi-leptonic B meson decay structure functions from lattice QCD

TL;DR

This paper introduces a non-perturbative lattice QCD method to calculate structure functions relevant to inclusive semi-leptonic B meson decays, enabling better theoretical understanding of these processes.

Contribution

It presents a novel approach for computing inclusive decay structure functions via lattice QCD, including a numerical test for the B_s -> X_c l nu mode at zero recoil.

Findings

Successfully identified ground state and excited state contributions.

Demonstrated feasibility of analytic continuation from unphysical to physical kinematics.

Performed calculations with physical charm quark mass and varied bottom quark mass.

Abstract

We propose a method to non-perturbatively calculate the forward-scattering matrix elements relevant to inclusive semi-leptonic B meson decays. Corresponding hadronic structure functions at unphysical kinematics are accessible through lattice QCD calculation of four-point correlation functions. The unphysical kinematical point may be reached by analytic continuation from the physical differential decay rate. A numerical test is performed for the B_s -> X_c l nu mode in the zero-recoil limit. We use lattice ensembles generated with 2+1 dynamical quark flavors. The valence charm quark mass is tuned to its physical value, while the bottom quark mass is varied in the range (1.56-2.44)m_c. From the numerical results we can identify the contributions of the ground state D_s^(*) meson as well as those of excited states or continuum states.

Tables

Table 1. Table 1: Lattice ensembles used in this work. They are generated with 2+1 flavors of dynamical Möbius domain-wall fermions. a ​ m u ​ d 𝑎 subscript 𝑚 𝑢 𝑑 am_{ud} and a ​ m s 𝑎 subscript 𝑚 𝑠 am_{s} are the light and strange quark masses included as the sea quark. The heavy quark masses a ​ m c 𝑎 subscript 𝑚 𝑐 am_{c} and a ​ m b 𝑎 subscript 𝑚 𝑏 am_{b} are only in the valence sector. The charm quark mass is tuned to its physical value, while the bottom quark mass is chosen such that they are 1.56 or 2.44 times heavier than charm. The number of statistical samples N cfg subscript 𝑁 cfg N_{\mathrm{cfg}} for each ensemble is given in the last column.

$\beta$	$a$ [fm]	$L^3\times T$	$a{m}_{ud}$	$a{m}_s$	$a{m}_c$	$a{m}_b$	$N_{\mathrm{cfg}}$
4.17	0.080	32${}_{}{}^3\times$64	0.007	0.040	0.44037	0.68808	50
4.35	0.055	48${}_{}{}^3\times$96	0.0042	0.0250	0.27287	0.42636	25
						0.66619	50
4.47	0.044	64${}_{}{}^3\times$128	0.0030	0.0150	0.210476	0.328869	25
						0.5138574	50

Table 2. Table 2: D s ( ∗ ) superscript subscript 𝐷 𝑠 D_{s}^{(*)} and B s subscript 𝐵 𝑠 B_{s} masses measured for each heavy quark mass parameters. The numerical data are those from Fahy:2017enl , which has more statistics than Table tab:lattices. The numbers are in the lattice unit.

$\beta$	$a{m}_c$	$a{M}_{D_s}$	$a{M}_{D_s^{\ast }}$	$a{m}_b$	$a{M}_{B_s}$
4.17	0.44037	0.8142(5)	0.8685(11)	0.68808	1.0624(11)
4.35	0.27287	0.5498(3)	0.5876(10)	0.42636	0.7220(5)
				0.66619	0.9524(8)
4.47	0.210476	0.4355(3)	0.4649(9)	0.328869	0.5737(4)
				0.513857	0.7682(8)

Table 3. Table 3: B s subscript 𝐵 𝑠 B_{s} meson mass a ​ M B s 𝑎 subscript 𝑀 subscript 𝐵 𝑠 aM_{B_{s}} and matrix element a 3 ​ Z L ​ L superscript 𝑎 3 superscript 𝑍 𝐿 𝐿 a^{3}Z^{LL} extracted from the SL and SS two-point functions. The results are obtained assuming a saturation by the ground state in an interval [ t min , t max ] subscript 𝑡 min subscript 𝑡 max [t_{\mathrm{min}},t_{\mathrm{max}}] . The last column shows the inverse vector current renormalization constant Z V ​ b ​ b − 1 superscript subscript 𝑍 𝑉 𝑏 𝑏 1 Z_{Vbb}^{-1} for the heavy-heavy current.

$\beta$	$a{m}_b$	$[t_{\min },t_{\max }]$	$a{M}_{B_s}$	$a^3{Z}^{LL}$	$Z_{Vbb}^{-1}$
4.17	0.68808	[12,20]	1.0607(20)	0.01023(33)	0.8323(51)
4.35	0.42636	[18,30]	0.7234(22)	0.00464(15)	0.9896(18)
	0.66619		0.9554(22)	0.00386(14)	0.7986(35)
4.47	0.328869	[24,40]	0.5767(15)	0.00262(10)	1.0088(11)
	0.513857		0.7674(16)	0.00214(09)	0.9274(34)

Table 4. Table 4: Numerical results for | D μ ​ ν J ​ J ​ ( ω ) | subscript superscript 𝐷 𝐽 𝐽 𝜇 𝜈 𝜔 |D^{JJ}_{\mu\nu}(\omega)| at ω = 0 𝜔 0 \omega=0 .

$\beta$	$a{m}_b$	$\left|D_{00}^{VV}\left(0\right)\right|$	$\left|D_{kk}^{VV}\left(0\right)\right|$	$\left|D_{00}^{AA}\left(0\right)\right|$	$\left|D_{kk}^{AA}\left(0\right)\right|$
4.17	0.68808	0.926(66)	0.035(07)	0.048(10)	0.805(59)
4.35	0.42636	0.924(69)	0.036(05)	0.034(10)	0.774(57)
	0.66619	1.055(92)	0.070(11)	0.060(18)	0.881(90)
4.47	0.328869	1.021(66)	0.041(06)	0.056(12)	0.842(58)
	0.513857	1.068(94)	0.065(10)	0.053(18)	0.871(81)

Table 5. Table 5: Zero-recoil form factors | h + ​ ( 1 ) | subscript ℎ 1 |h_{+}(1)| and | h A 1 ​ ( 1 ) | subscript ℎ subscript 𝐴 1 1 |h_{A_{1}}(1)| extracted from | D μ ​ ν J ​ J ​ ( ω → m D s ( ∗ ) ) | subscript superscript 𝐷 𝐽 𝐽 𝜇 𝜈 → 𝜔 subscript 𝑚 superscript subscript 𝐷 𝑠 |D^{JJ}_{\mu\nu}(\omega\to m_{D_{s}^{(*)}})| .

$\beta$	$a{m}_b$	$\left|h_{+}\left(1\right)\right|$	$\left|h_{A_1}\left(1\right)\right|$
4.17	0.68808	0.969(28)	0.898(44)
4.35	0.42636	0.946(25)	0.843(45)
	0.66619	0.967(30)	0.840(45)
4.47	0.328869	0.987(30)	0.883(60)
	0.513857	0.981(33)	0.857(56)

Equations

$$H_W=V_{qQ}\frac{G_F}{\sqrt{2}}\bar{\ell }{\Gamma }^{\mu }{\nu }_{\ell }\cdot J_{\mu }.$$

$$|\mathcal{M}{|}^2=|V_{qQ}{|}^2{G}_F^2{M}_B{l}^{\mu \nu }{W}_{\mu \nu }.$$

$$W_{\mu \nu }=\sum _X(2\pi {)}^3{\delta }^4(p_B-q-p_X)\frac{1}{2M_B}⟨B(p_B)|J_{\mu }^{†}(0)|X⟩⟨X|J_{\nu }(0)|B(p_B)⟩.$$

$$W_{\mu \nu }=-W_1{g}_{\mu \nu }+W_2{v}_{\mu }{v}_{\nu }+i{W}_3{ϵ}_{\mu \nu \alpha \beta }{v}^{\alpha }{q}^{\beta }+W_4{q}_{\mu }{q}_{\nu }+W_5(q_{\nu }{v}_{\mu }+q_{\mu }{v}_{\nu }).$$

$$T_{\mu \nu }=i\int d^{4}x{e}^{-iqx}\frac{1}{2M_B}⟨B|T\{J_{\mu }^{†}(x)J_{\nu }(0)\}B⟩,$$

$$T_{\mu \nu }=-T_1{g}_{\mu \nu }+T_2{v}_{\mu }{v}_{\nu }+i{T}_3{ϵ}_{\mu \nu \alpha \beta }{v}^{\alpha }{q}^{\beta }+T_4{q}_{\mu }{q}_{\nu }+T_5(q_{\nu }{v}_{\mu }+q_{\mu }{v}_{\nu }).$$

$$-\frac{1}{\pi }\mathrm{Im}{T}_i=W_i.$$

$$-\infty <v\cdot q<\frac{1}{2M_B}(M_B^2+q^2-m_X^2)$$

$$\frac{1}{2M_B}((2M_B+M_X{)}^2-q^2-M_B^2)<v\cdot q<\infty .$$

$$\frac{1}{2M_B}(M_B^2+q^2-m_X^2)<v\cdot q.$$

$$T(v\cdot q)=\frac{1}{\pi }{\int }_{-\infty }^{(v\cdot q{)}_{\max }}d(v\cdot q^{\prime })\frac{\mathrm{Im}T(v\cdot q^{\prime })}{v\cdot q^{\prime }-v\cdot q}.$$

$${⟨\frac{1}{q_0-(m_B-\omega )}⟩}_{\mathrm{fixed}{q}^2}\mathrm{or}{⟨\frac{1}{q_0-(m_B-\omega )}⟩}_{\mathrm{fixed}{\mathbf{q}}^2}$$

$${⟨{\left(\frac{1}{q_0-(m_B-\omega )}\right)}^2⟩}_{\mathrm{fixed}{q}^2}\mathrm{or}{⟨{\left(\frac{1}{q_0-(m_B-\omega )}\right)}^2⟩}_{\mathrm{fixed}{\mathbf{q}}^2}.$$

$$C_{\mu \nu }^{SJJS}(t_{\mathrm{snk}},t_1,t_2,t_{\mathrm{src}})=\sum _{\mathbf{x}}⟨P^S(\mathbf{x},t_{\mathrm{snk}}){\tilde{J}}_{\mu }^{†}(\mathbf{q},t_1){\tilde{J}}_{\nu }(\mathbf{q},t_2)P^{S†}(0,t_{\mathrm{src}})⟩,$$

$$C^{XY}(t_{\mathrm{snk}},t_{\mathrm{src}})=\sum _{\mathbf{x}}⟨P^X(\mathbf{x},t_{\mathrm{snk}})P^{Y†}(0,t_{\mathrm{src}})⟩,$$

$$\frac{C_{\mu \nu }^{SJJS}(t_{\mathrm{snk}},t_1,t_2,t_{\mathrm{src}})}{C^{SL}(t_{\mathrm{snk}},t_1)C^{LS}(t_2,t_{\mathrm{src}})}⟶\frac{\frac{1}{2M_B}⟨B(0)|{\tilde{J}}_{\mu }(\mathbf{q},t_1{)}^{†}{\tilde{J}}_{\nu }(\mathbf{q},t_2)|B(0)⟩}{\frac{1}{2M_B}|⟨0|P^L|B(0)⟩{|}^2,}$$

$$C_{\mu \nu }^{JJ}(t;\mathbf{q})=\int d^3\mathbf{x}{e}^{i\mathbf{q}\cdot \mathbf{x}}\frac{1}{2M_B}⟨B(0)|J_{\mu }^{†}(\mathbf{x},t)J_{\nu }(0)|B(0)⟩,$$

$$T_{\mu \nu }^{JJ}(\omega ,\mathbf{q})={\int }_0^{\infty }dt{e}^{\omega t}{C}_{\mu \nu }^{JJ}(t;\mathbf{q}).$$

$$\frac{d}{d\omega }{T}_{\mu \nu }^{JJ}(\omega ,\mathbf{q})={\int }_0^{\infty }dtt{e}^{\omega t}{C}_{\mu \nu }^{JJ}(t;\mathbf{q}),$$

$$C_{\mu \nu }^{JJ}(t;0)=\frac{1}{2M_B}\sum _X⟨B(0)|J_{\mu }^{†}|X(0)⟩\frac{e^{-E_{X}t}}{2E_X}⟨X(0)|J_{\nu }|B(0)⟩,$$

$$⟨D(0)|V^0|B(0)⟩$$

$$⟨D^{\ast }(0)|A^k|B(0)⟩$$

$$⟨D(0)|V^0|B(0)⟩$$

$$⟨D^{\ast }(0)|A^k|B(0)⟩$$

$$T_{00}^{VV}(\omega ,0)$$

$$T_{kk}^{AA}(\omega ,0)$$

$${\mu }_{\pi }^2$$

$${\mu }_G^2$$

$$\left\{\begin{array}[c]{l}T_{00}=-T_{1}+T_{2}+T_{4}q_{0}^{2}+2T_{5}q_{0},\\ T_{0k}=T_{k0}=0,\\ T_{kk}=T_{1},\\ T_{kl}=0\;\;\;(k\not=l).\end{array}\right.$$

$$T_1^{VV}$$

$$T_1^{AA}$$

$$T_2^{VV}=T_2^{AA}$$

$$T_5^{VV}=T_5^{AA}$$

$$T_{kk}^{AA}=T_{00}^{VV}$$

Full text

1 Introduction

There is a long-standing tension between the inclusive and exclusive determinations of the Cabbibo-Kobayashi-Maskawa (CKM) matrix elements $|V_{cb}|$ and $|V_{ub}|$. For instance, the Particle Data Book (2016) quotes $|V_{cb}|$ = $(42.2\pm 0.8)\times 10^{-3}$ for the inclusive determination, while the exclusive determination from $B\to D^{(*)}\ell\nu$ modes gives $(39.2\pm 0.7)\times 10^{-3}$. They are separated by 2.8 standard deviations. Reflecting this tension, the quoted average for $|V_{cb}|$ has a larger uncertainty, $(40.5\pm 1.5)\times 10^{-3}$ Olive:2016xmw , than the individual inclusive and exclusive determinations. Here, the dominant source of error is “unknown” until the reason for the tension is understood. The situation is similar for $|V_{ub}|$. (See, for instance Dingfelder:2016twb for a summary of the recent theoretical and experimental status.)

The exclusive determination uses the differential decay rates of $B$ meson to the ground state $D$ or $D^{*}$ meson, $B\to D^{(*)}\ell\nu$, combined with the corresponding form factors calculated using lattice Quantum Chromodynamics (LQCD). The precision of the lattice calculation will be improved as more calculations on finer lattices and with more statistics become available in the near future. The inclusive determination, on the other hand, is obtained from experimental measurements of all possible semi-leptonic final states compared to the quark level calculation of the corresponding decay rate in QCD. This theoretical calculation involves perturbative expansion, heavy quark expansion and other techniques, each of which would need non-trivial steps to improve. Furthermore, the potential error due to the quark-hadron duality, which is a common assumption in the inclusive calculations, is hard to quantify. This situation makes the understanding of the possible cause of the tension rather difficult.

In this paper, we propose a method to non-perturbatively calculate the amplitudes related to the inclusive decay rate using the LQCD technique. Our goal is to identify the source of the tension by comparing the inclusive and exclusive determinations on a single set of lattice calculations. As a first step, this paper describes the framework to perform such a comparison.

The theoretical formulation for the inclusive semi-leptonic $B$ decay Chay:1990da is analogous to that of nucleon deep inelastic scattering, where the cross section is written in terms of forward-scattering matrix elements of the nucleon. For the inclusive $B$ decay, the differential decay rate is related to the imaginary part of the forward matrix elements through the optical theorem. For the $b\to c\ell\nu$ channel, one considers the matrix elements of bi-local operators $J_{\mu}^{\dagger}(x)J_{\nu}(y)$ with $J_{\mu}$ the electroweak $V-A$ current to describe the $b\to c$ transition. In the kinematical region where hadronic resonances appear in the final state, perturbation theory is not applicable to directly calculate the differential decay rate. Rather, under the assumption of “quark-hadron duality” the rate smeared over some kinematical variables are evaluated using perturbation theory Poggio:1975af . For semi-leptonic decays, such smearing is naturally provided by the integral over the phase space to obtain the total decay rate. The question is, then, whether the amount of smearing is enough to ensure the use of the perturbation theory. In this work, on the other hand, we propose to focus on different kinematical region that is off-shell and no resonances directly appear.

We consider unphysical kinematical setting where the final hadronic system does not have enough energy to be on-shell. Using the analyticity of the decay amplitude, the forward-scattering matrix element of the $B$ meson in this kinematical region can be related to an integral of the physical decay amplitude, which is an imaginary part of the matrix element. In this unphysical kinematical region, perturbation theory may be used, and more importantly a direct lattice QCD calculation is possible. We are thus able to apply the different approaches for the same physical quantity, which may give some clue to understand the source of the tension between the inclusive and exclusive analyses.

In this paper we focus on the zero-recoil limit of the $b\to c\ell\nu$ inclusive decays to demonstrate how the method works, although the theoretical framework is general and applies to non-zero recoil processes including both $b\to c\ell\nu$ and $b\to u\ell\nu$ decay modes. The zero-recoil limit has its own interest because the analysis of the inclusive decay rate has been used to derive a bound and an estimate of the zero-recoil form factor $h_{A_{1}}(1)$ of the $B\to D^{(*)}\ell\nu$ decay Bigi:1994ga ; Shifman:1994jh ; Kapustin:1996dy ; Gambino:2010bp ; Gambino:2012rd . This so-called zero-recoil sum rule relates the perturbatively calculated forward-scattering matrix elements (with heavy quark expansion) to the sum of resonance contributions including the lowest-lying $D^{(*)}$ meson. Since the contributions of each hadronic state are positive, it may give an upper bound on the $B\to D^{(*)}$ semi-leptonic form factor. The analysis of Shifman:1994jh implies ${\cal F}(1)<0.94$, while a recent analysis Gambino:2012rd gives ${\cal F}(1)<0.92$; there is also a caution that with perturbative corrections the bound is much weaker ${\cal F}(1)<0.98$ Kapustin:1996dy . Estimating the size of the excited state contributions, the central value could also be evaluated: ${\cal F}(1)\simeq 0.89(3)$ Shifman:1994jh , ${\cal F}(1)\approx 0.86$ Gambino:2012rd . The recent lattice calculation by the Fermilab-MILC collaboration Bailey:2014tva , on the other hand, gives ${\cal F}(1)=0.906(4)(12)$, which is on the upper side of the sum rule estimate. Since the lower ${\cal F}(1)$ implies larger $|V_{cb}|$ extracted, the numbers obtained by the sum rule approach make $|V_{cb}|$ consistent with its inclusive determination. Full analysis of the unphysical matrix elements to be proposed in this paper will allow us to directly compare the sum rule and lattice approaches and may give some hints to understand the tension between the inclusive and exclusive determinations.

Lattice QCD has most widely been used for the calculation of lowest-lying hadron masses and matrix elements. This is achieved by extracting the contribution of the ground state from (Euclidean) long-distance correlation functions. For the calculation of the inclusive decays, this is insufficient as the excited states and scattering states are discarded. The method we apply in this work is based on the analytic continuation from the space-like momentum region to the time-like. The first such proposal was made by Ji and Jung for a calculation of the hadronic structure function of photon Ji:2001wha , and the method has been applied to the two-photon decays of Charmonium Dudek:2006ut as well as pion Feng:2012ck ; Gerardin:2016cqj . An application to the hadronic vacuum polarization function for muon $g-2$ has also been developed Feng:2013xsa . These are all related to the processes of non-QCD external states with its four-momentum specified. The inclusive $B$ meson decay is also in this class, i.e. the emitted virtual $W(\to\ell\nu)$ has a specified four-momentum $q_{\mu}$.

We compute four-point functions on the lattice. Two operators create and annihilate the $B$ meson, while the other two represent the weak currents with the momentum transfer $q_{\mu}=(q_{0},\mathbf{q})$. The matrix elements obtained from such four-point functions are integrated over time separation $t$ between the two currents with a weight factor $e^{\omega t}$. Then, we obtain the matrix element for a fixed value of the momentum transfer $(m_{B}-\omega,\mathbf{q})$, as long as $\omega$ is lower than the lowest energy of the final state, i.e. the $D^{(*)}$ meson. In other words, the momentum flowing into the hadronic final state is $(\omega,-\mathbf{q})$, and the above condition for $\omega$ implies unphysical kinematics. This is the kinematical region one can access in this method. For the purpose of the comparison between the lattice calculation and the continuum perturbation theory, this unphysical setup provides a convenient testing ground.

The matrix elements have two independent kinematical variables $q^{2}$ and $v\cdot q$. Here, $v_{\mu}$ denotes the four-velocity of the initial $B$ meson. Connection to the differential decay rate may be provided by an analytic continuation in the complex plane of $v\cdot q$ while $q^{2}$ is fixed. One may also take $(v\cdot q)^{2}-q^{2}$ and $v\cdot q$ as two independent variables. In the rest frame of the $B$ meson they correspond to $\mathbf{q}^{2}$ (three-momentum squared) and $m_{B}-\omega$, which have direct correspondences to the momentum and energy of the final hadronic system.

The method introduced in this work is analogous to those for the lattice calculation of deep inelastic scattering structure functions Liu:1993cv ; Liu:1998um ; Liu:1999ak . The difference is in the kinematical setup, i.e. by considering the unphysical kinematical region we are able to obtain the relevant amplitude unambiguously. (Our setup corresponds to the calculation of the structure function $W(Q^{2},\nu)$ at small $\nu$, or equivalently the Bjorken $x$ greater than 1, which is unphysical.) More closely related, a method to calculate the $B$ meson shape function or light-cone wave function was proposed some years ago Aglietti:1998mz ; Aglietti:1998ur . The main emphasis was in the region of large recoil momentum, where the lattice calculation is substantially more difficult, and to our knowledge there has been no successful numerical calculation performed using this method. The proposal in this paper includes the lattice calculation in such kinematical range but is not limited to it. We are able to proceed step by step from a less expensive case of the zero-recoil momentum.

This paper describes the method of lattice calculation as well as a numerical test, which serves as a proof of concept. The lattice ensembles used for this purpose are among those of state-of-the-art. They are generated with 2+1 flavors of dynamical Möbius domain-wall fermions. The lattice spacing covered corresponds to the inverse lattice spacing $1/a\simeq$ 2.45, 3.61 and 4.50 GeV. Up and down quarks are slightly heavier than their physical mass, while the strange quark mass is approximately set to the physical value. In this initial study, we choose the strange quark as the spectator quark, so that the initial state is $B_{s}$ and the final states have a $c\bar{s}$ quantum number with its ground states being $D_{s}^{(*)}$. With the relatively large lattice cutoff of our ensemble, the charm quark can be treated as in the same manner as for light quarks. Our recent calculation of the $D_{(s)}$ meson decay constant Fahy:2015xka suggests that the discretization effect is well under control.

For the bottom quark, on the other hand, the lattice cutoff is not sufficiently large to simply ignore the discretization effects of $O((am_{b})^{2})$ originating from large $b$ quark mass $m_{b}$. It appears, for instance, as a large deviation of the renormalization constant $Z_{V}$ of vector current $\bar{Q}\gamma_{\mu}Q$ from being 1. (Without using the conserved current on the lattice, $Z_{V}$ can deviate from 1 in general. Our analysis of light-quark current correlator indicates that $Z_{V}$ for massless quarks is close to 1 within 5%.) We renormalize the heavy quark current partially using the measured matrix element $\langle B|\bar{Q}\gamma_{0}Q|B\rangle$. Remaining discretization effects are estimated using the calculations at three lattice spacings. In order to keep them under control, we limit the value of $m_{b}$ so that $am_{b}$ is below 0.6–0.7. The resulting “$b$ quark” is only 1.56–2.44 times heavier than charm, and still lighter than the physical $b$ quark. Still, because of the heavy quark symmetry, the motion of the $b$ quark would be well approximated by the lighter “$b$ quark”. A more realistic calculation is left for future studies.

The rest of the paper is organized as follows. Section 2 summarizes the kinematics of the semi-leptonic $B$ meson decays, and Section 3 introduces the unphysical kinematics we treat. Our strategy for the lattice calculation is given in Section 4, and some estimates in the zero-recoil limit are discussed in Section 5. Lattice results in Section 6 demonstrate how the method works by numerical calculations. Future directions are discussed in Section 7.

2 Inclusive $B$ meson decays

We consider the semi-leptonic decay process $B(p_{B})\to X(p_{X})\ell(p_{\ell})\nu_{\ell}(p_{\nu})$. The lepton $\ell$ and neutrino $\nu_{\ell}$ are generated through a virtual $W$ which has a four momentum $q=p_{\ell}+p_{\nu}$. $X$ represents the hadronic final state, which may either be a single hadron like $D^{(*)}$ meson or multibody states. The momentum $p_{X}$ stands for the total four-momentum of the whole hadronic final state. The initial $B$ meson could either be $B^{\pm,0}$ or $B_{s}$; the final state $X$ carries strangeness when the initial is $B_{s}$. The discussions in the following closely follow those of Chay:1990da ; Manohar:1993qn ; Blok:1993va .

The momentum of the initial and final states are specified above. The momentum of the lepton system $q=p_{\ell}+p_{\nu}$ also plays the role of the momentum transfer $q=p_{B}-p_{X}$. The electroweak Hamiltonian for this process is given by

[TABLE]

Here the hadronic current $J_{\mu}$ also has the $V-A$ structure $J_{\mu}=\bar{q}\Gamma_{\mu}Q$ with $\Gamma_{\mu}\equiv\gamma_{\mu}(1-\gamma_{5})$. ($q$ is either $u$ or $c$; $Q$ denotes $b$.) $V_{qQ}$ denotes a Cabibbo-Kobayashi-Maskawa matrix element, either $V_{cb}$ or $V_{ub}$.

The decay rate is given by an absolute value squared of the amplitude, and has a decomposition analogous to the deep inelastic scattering analysis:

[TABLE]

The leptonic tensor $l^{\mu\nu}$ is known, and the hadronic tensor is written as a sum of matrix elements

[TABLE]

Here the sum over the final state $X$ includes an integral over its four-momentum $p_{X}$. One can introduce the structure functions $W_{i}$ to parametrize the hadronic tensor:

[TABLE]

Here, $v_{\mu}=(p_{B})_{\mu}/M_{B}$ is the four-velocity of the initial $B$ meson. The $W_{i}$’s are functions of two invariant variables $v\cdot q$ and $q^{2}$, and have mass dimension $-1$ ($W_{1}$, $W_{2}$), $-2$ ($W_{3}$, $W_{5}$) or $-3$ ($W_{4}$). Among five structure functions, three ($W_{1}$, $W_{2}$ and $W_{3}$) contribute to the decay amplitude for light leptons ($\ell$ = $e$ or $\mu$), and others are relevant only for $\ell=\tau$.

One may rewrite the hadronic tensor $W_{\mu\nu}$ using the forward scattering matrix element

[TABLE]

and the corresponding structure functions $T_{i}$ defined similarly as

[TABLE]

By inserting the complete set of states between the currents, one can see that the imaginary part of $T_{i}$ gives the hadronic tensor,

[TABLE]

Namely, $T_{i}$’s are defined for more general external momenta and analytically continued from $W_{i}$’s.

For a given value of $q^{2}$, the cut for $T_{i}$’s in the complex plane of $v\cdot q$ runs on the real axis in the region

[TABLE]

as well as

[TABLE]

This analytic structure is depicted in Figure 1, where the cuts in the complex $v\cdot q$ plane are shown by thick lines.

The process $B\to X\ell\nu_{\ell}$ occurs only in part of (8), i.e. $\sqrt{q^{2}}<v\cdot q<\frac{1}{2M_{B}}(M_{B}^{2}+q^{2}-m_{X}^{2})$. The upper cut (9) is also unphysical, corresponding to the process of $b\to\bar{c}bb$, which is always far apart from the physical cut for the $b\to c$ process. When we consider the $b\to u$ process near $q^{2}\approx 0$, the gap between the two cuts becomes narrow.

3 Structure functions at unphysical kinematics

Not the entire kinematical region of $T_{i}$’s is accessible by LQCD, which is defined on the Euclidean space. In particular, the imaginary part never shows up. Instead, we consider the region on the real axis between (8) and (9). Namely,

[TABLE]

(The upper cut is too far to be relevant in this work.) At the rest frame of the initial $B$ meson, this condition implies $p_{X}^{2}=(M_{B}-q_{0})^{2}<m_{X}^{2}$. The lowest $m_{X}$ is the mass of the ground state $D_{(s)}$ meson for $b\to c$ (or pion for $b\to u$). The region (10) corresponds to the case where the final state energy is not sufficient to become on-shell.

For the connection to the physical decay amplitude, we need to perform a contour integral of the form

[TABLE]

Here, $v\cdot q$ on the left hand side corresponds to the unphysical kinematics specified in (10), and the integral over $(v\cdot q^{\prime})$ is along the physical cut. The upper limit of the integral is $(v\cdot q)_{\mathrm{max}}\equiv(M_{B}^{2}+q^{2}-m_{X}^{2})/2M_{B}$. In the integrand, the experimental result may be inserted for $\mathrm{Im}T$, but only in the kinematically accessible region $v\cdot q>\sqrt{q^{2}}$. For the inaccessible kinematical region, one needs to use perturbation theory, which should be well-behaved as it is far away from the hadronic resonances.

Figure 2 shows the phase space of the $b\to c$ semi-leptonic decay on the plane of invariant variables $v\cdot q$ and $q^{2}$. The axes are normalized by $m_{b}$ as $v\cdot\hat{q}=v\cdot q/m_{b}$ and $\hat{q}^{2}=q^{2}/m_{b}^{2}$. The gray region represents the physically allowed range $\sqrt{\hat{q}^{2}}\leq v\cdot\hat{q}\leq(1+\hat{q}^{2}-m_{c}^{2}/m_{b}^{2})/2$ with $0\leq\hat{q}^{2}\leq(1-m_{c}/m_{b})^{2}$ Falk:1995me . The integral in (11) is along a horizontal line at a fixed $q^{2}$. It is clear that other than the region of small $q^{2}$ the physically available region is quite limited and one must rely on the perturbative calculation for the region below $v\cdot\hat{q}\leq\sqrt{\hat{q}^{2}}$. One may instead consider a plane of $v\cdot q$ and $(v\cdot q)^{2}-q^{2}$, which is shown in Figure 3. In the rest frame of the initial $B$ meson, these variables correspond to $q_{0}$ and $\mathbf{q}^{2}$, respectively. For small $\mathbf{q}^{2}$, such as the zero-recoil limit considered later in this work, the contour integral along the fixed $\mathbf{q}^{2}$ would be more useful.

Strictly speaking, the contribution from the $\bar{c}bb$ cut must be taken into account in the contour integral (11), which amounts to add another integral starting from $((2M_{B}+M_{X})^{2}-q^{2}-M_{B}^{2})/2M_{B}$. Such contribution must be negligible because of the large suppression factor appearing in the denominator of the integrand. We ignore the corresponding contribution also in the lattice calculation, so that the treatment is consistent.

The integral (11) of the experimentally observed differential decay rate corresponds to an inverse hadronic energy moment of the form

[TABLE]

in the rest frame of the initial $B$ meson. The variable $\omega$ is an arbitrary parameter between zero and the lowest energy of the final hadronic system at a fixed $q^{2}$ (or at a fixed $\mathbf{q}^{2}$). The integral is more inclusive when $\omega$ is lower. Near the upper limit of $\omega$, the moment is dominated by the ground state of the final hadron. Later in this paper, we also introduce a derivative of the structure functions with respect to $\omega$, which is obtained as a second moment

[TABLE]

In (12) and (13) the moments are considered on the double differential decay rate $d^{2}\Gamma/dq^{2}dq_{0}$, which is non-zero in the phase space shown in Figure 2 or 3. (For the definition, see Falk:1995me ; Falk:1997jq ; Ligeti:2014kia for instance.) So far, in the literature, the moments of hadron energy and invariant mass as well as the lepton energy have been considered; our proposal is to analyze the inverse moments (12) and (13) at sufficiently small $\omega$, instead, to extract $|V_{cb}|$ or $|V_{ub}|$. To actually extract the moments from the experimental data is beyond the scope of this work.

The structure functions $T_{i}$ have been calculated within the heavy quark expansion approach. At the tree-level, the explicit form is given in the appendix of Blok:1993va . One-loop or even two-loop calculations have also been carried out Trott:2004xc ; Aquila:2005hq ; Biswas:2009rb , but they only concern the differential decay rates (or the imaginary part of the structure functions), and one needs to perform the contour integral to relate them to the unphysical kinematical region.

4 Lattice calculation strategy

In this section, we describe the method to extract $T_{i}$’s from a four-point function calculated on the lattice. Although we take the $B\to D^{(*)}\ell\nu$ channel to be specific, the extension to other related channels is straightforward.

We consider the four-point function of the form

[TABLE]

where $P^{S}$ is a smeared pseudo-scalar density operator to create/annihilate the initial $B$ meson at rest. The inserted currents $\tilde{J}_{\mu}$ are either vector or axial-vector $b\to c$ current and assumed to carry the spatial momentum projection $\sum_{\mathbf{x}_{1}}e^{i\mathbf{q}\cdot\mathbf{x}_{1}}J(\mathbf{x}_{1},t_{1})$. Thus, the mass dimension of $\tilde{J}_{\mu}$ is zero. The quark-line diagram representing (14) is shown in Figure 4.

We also define the $B$ meson two-point functions

[TABLE]

with $XY$ denoting a combination of local ($L$) or smeared ($S$) operators at the source (src) and sink (snk). We then consider a ratio

[TABLE]

where the asymptotic form on the right hand side is reached at long time separations $t_{\mathrm{snk}}-t_{1}$, $t_{2}-t_{\mathrm{src}}\gg 1$ by inserting the complete set of states for the initial and final states. The amplitude to create/annihilate $B$ meson is canceled by the denominator, as well as the exponential dumping of the correlators due to the $b$ quark mass. In practice, we may assume the translational invariance and time reversal to rewrite the two-point functions in the denominator by $C^{LS}(t,t_{\mathrm{src}})$.

The matrix element $|\langle 0|P^{L}|B(\mathbf{0})\rangle|$ appearing in (16) may be obtained by combining the $LS$ and $SS$ two-point correlators as usual. After multiplying this factor to (16), we may extract

[TABLE]

which has the form close to the structure function defined in (5).

The remaining integral over time with $e^{-iq_{0}x_{0}}$ in (5) has to be carried out. Since the “time” on the lattice is the Eudlidean time, this factor actually corresponds to $e^{\omega t}$. We then arrive at

[TABLE]

Absorbing the time dependence of the $B$ meson states, $\omega$ is related to $q_{0}$ as $\omega=m_{B}-q_{0}$. The integral region is limited to $t>0$, because the negative $t$ region corresponds to the $\bar{c}bb$ state, which we ignore.

Since the function $C_{\mu\nu}^{JJ}(t;\mathbf{q})$ falls off exponentially $e^{-E_{X}(\mathbf{q})t}$ with $E_{X}(\mathbf{q})$ the ground state energy at large time separations, $\omega$ must be smaller than $E_{X}(\mathbf{q})$ for the integral to be convergent. Indeed, this is the same as the condition $p_{X}^{2}=(M_{B}-q_{0})^{2}<m_{X}^{2}$, which defines the unphysical kinematical region (10).

This kind of “Fourier transform” can be justified by using the analytic continuation when there are no other singularities. A good example is the vacuum polarization function $\langle J_{\mu}J_{\nu}\rangle$, which appears in the analysis of muon $g-2$. The analytic continuation is possible for $q^{2}$ below the threshold of $m_{\rho}^{2}$ (or the lowest energy state $\pi\pi$). Another way to understand it is to consider $n$-th derivatives with respect to $\omega$ to construct a Taylor expansion around $\omega=0$. It gives temporal moments of the correlators that is defined without the factor $e^{\omega t}$. Then, as far as the Taylor series converges, it can be extended to real and imaginary $\omega$; on the real axis, it can be extended until the function hits the first singularity, which is the ground state pole. Some details are discussed in Feng:2013xsa for an application to the vacuum polarization function.

For the inclusive $B$ meson decays, there is another cut representing the $\bar{c}bb$ state. This state corresponds to the case where the time-ordering of the current insertions are reversed, i.e. $\hat{J}_{\mathbf{p}}\hat{J}_{\mathbf{p}}^{\dagger}$ rather than $\hat{J}_{\mathbf{p}}^{\dagger}\hat{J}_{\mathbf{p}}$. Since the corresponding intermediate states are far more suppressed by the factor $e^{-2m_{b}(t_{1}-t_{2})}$ for heavy $b$ quark, their contribution is negligible. We will neglect the contribution from this state consistently for both the lattice computation and the contour integral.

The integral (18) contains the point $t=0$, where the two currents sit on the same time-slice and generate a contact term. Possible ultraviolet divergence due to the contact term does not survive in $T^{JJ}_{\mu\nu}$ after the four-dimensional integral over space-time, since the divergence scales as $1/|x|^{3}$ for small separation $x$. Still, we try to avoid the $t=0$ time slice, as it contains the contribution from the distant cut $\bar{c}bb$ which we ignore. Since the fractions of the $t=0$ contribution to the $c$ states in $t\geq 0$ and to the $\bar{c}bb$ states in $t\leq 0$ are unknown, this would be unavoidable. The $t=0$ time-slice may be avoided by introducing a derivative

[TABLE]

which we mainly analyze in this work. The derivative with respect to $\omega$ can be extracted from the experimental data through the second moment defined in (13).

5 Zero recoil limit

In this section we summarize some useful relations obtained in the limit of zero recoil $\mathbf{q}=0$. They are from both the hadronic (Section 5.1) and quark (Section 5.2) pictures. On the hadronic side, we consider the limit where the ground state saturates the final states. The structure functions can then be written in terms of the corresponding decay form factors. On the quark side, the estimates from the heavy quark expansion are summarized.

5.1 Contribution of the ground state $D^{(*)}$ meson

Consider the function $C^{JJ}_{\mu\nu}(t;\mathbf{q})$ defined in (17). By setting $\mathbf{q}=0$ and inserting the complete set of states between the currents, it is written as

[TABLE]

where the state $X(\mathbf{0})$ has zero total spatial momentum. Among the final states the ground state contribution would be most significant. The ground state is $D$ for the temporal vector current $V^{0}$ and $D^{*}$ for the spatial axial-vector current $A^{k}$. The corresponding matrix elements are expressed by the form factors $f^{0}(q^{2})$ and $A_{1}(q^{2})$ as

[TABLE]

in the conventional definition of the form factors. Here, $\varepsilon^{*}$ denotes the polarization vector of the vector particle. The maximum momentum transfer $q^{2}_{\mathrm{max}}$ is $q^{2}_{\mathrm{max}}=(M_{B}-M_{D^{(*)}})^{2}$. With the notation of heavy quark effective theory, they are

[TABLE]

The argument $v\cdot v^{\prime}=1$ of $h_{+}$ and $h_{A_{1}}$ means the zero-recoil limit $v=v^{\prime}$, i.e. the initial and final heavy mesons have the same four-velocity $v$ and $v^{\prime}$. The form factors $h_{+}(v\cdot v^{\prime})$ and $h_{A_{1}}(v\cdot v^{\prime})$ become identical in the heavy quark limit, which is called the Isgur-Wise function. It is normalized to 1 in the zero-recoil limit.

Using them, the ground state contribution to $C^{JJ}_{\mu\nu}(t;\mathbf{0})$ in (20) can simply be written as $|h(1)|^{2}\exp(-M_{D^{(*)}}t)$. ($h(1)$ represents either $h_{+}(1)$ or $h_{A_{1}}(1)$ depending on the channel.) Then, the ground-state contribution to the structure function $T^{JJ}_{\mu\nu}(\omega,\mathbf{0})$ defined in (18) are written as

[TABLE]

The excited state contributions have similar structure with higher masses such as $M_{D^{(*)\prime}}$ and different form factors. Note that they are all constructive for $\omega$ smaller than $M_{D^{(*)}}$.

5.2 Quark level estimate with heavy quark expansion

The same set of structure functions can also be calculated using the heavy quark expansion. It is a form of the operator product expansion applied for the heavy hadron decays. The product of currents $J_{\mu}^{\dagger}(x)J_{\nu}(0)$ may be written in terms of local operators $\bar{Q}Q$, $\bar{Q}\mathbf{D}^{2}Q$, $\bar{Q}\sigma_{\mu\nu}D_{\mu\nu}Q$, and so on; the higher dimensional operators are suppressed by a power of $1/m_{b}$. The leading order corresponds to the free quark decay, and the next-to-leading order operators represent the motion of the $b$ quark inside the hadron. They are conventionally parametrized by two quantities:

[TABLE]

or equivalently $\lambda_{1}$ and $\lambda_{2}$ depending on the literature. (Precise definition of the operators including their renormalization scheme becomes important once the higher order corrections are included. In this paper we do not discuss the one-loop corrections, and the $1/m_{b}$ terms associated with these operators are included only to estimate their potential size.)

The results for the structure functions $T_{i}$ at the tree level are available in Blok:1993va including the first non-trivial correction in the heavy quark expansion. In the zero-recoil limit, the hadronic tensors reduce to

[TABLE]

At the tree-level and at the leading order of the heavy quark expansion, these functions are written as

[TABLE]

Others are zero at this order. Here $m_{q}$ denotes the quark mass in the final state, which is the charm quark mass in our case. We introduce a notation $\omega\equiv m_{Q}-q_{0}$, which is consistent with the one adopted in our lattice calculation. Then, for the components of the structure functions we obtain

[TABLE]

at this order.

One finds a pole at $\omega=m_{q}$ or $q_{0}=m_{Q}-m_{q}$ for $T_{00}^{VV}$ and $T_{kk}^{AA}$. These are also the channels for which the ground state dominates as shown in the previous subsection. For the real QCD, this pole extends to a cut towards lower values of $v\cdot q$.

Including the corrections of $O(1/m_{b})$ the expressions are lengthy, but given in the appendix of Blok:1993va .

6 Numerical tests

In this section we present a numerical test of the method to extract the structure functions from lattice calculations of the four-point functions. Lattice calculations are still unrealistic as the initial $b$ quark mass is lower than its physical value. The main purpose of this section is rather to demonstrate how the method works.

6.1 Setup

We use the lattice ensembles generated by the JLQCD collaboration. They contain the dynamical quark effects of 2+1 flavors of the Möbius domain-wall fermion formulation. For this particular analysis, we use the ensembles listed in Table 1. The light quark masses ($am_{ud}$ and $am_{s}$) are those in the ensemble generation, while heavy quarks ($am_{c}$ and $am_{b}$) are only in the valence sector.

In this work, we only consider the strange quark for the spectator quark; an extension to the up and down quarks is technically straightforward. For the heavy quarks, we use the same domain-wall fermion formulation as those for the light quarks. Discretization effects for the charm quark are under good control, as the previous calculation of the charm quark mass suggests Nakayama:2016atf . The extrapolation towards the bottom quark still needs extra care Fahy:2015xka ; Fahy:2017enl , where we introduced the reinterpretation of the heavy quark propagator according to ElKhadra:1996mp . In this work, since the propagation of the $b$ quark cancels between the numerator and denominator of (16), such reinterpretation does not affect the result.

The lattice spacing corresponds to the cutoff $a^{-1}$ of 2.453(4), 3.610(9), and 4.496(9) GeV for $\beta$ = 4.17, 4.35, and 4.47, respectively. They are determined through the Wilson flow scale $t_{0}^{1/2}$ with an input from Borsanyi:2012zs . The same ensembles have been used for various physical quantities including the charm quark mass determination through the charmonium correlator Nakayama:2016atf , calculation of heavy-light decay constants Fahy:2017enl , $D$ meson semi-leptonic form factors Kaneko:2017sct , chiral condensate through the Dirac operator spectrum Cossu:2016eqs , as well as the extraction of the $\eta^{\prime}$ mass from topological charge density correlators Fukaya:2015ara . The code set Iroiro++ has been used for writing the simulation codes Cossu:2013ola .

We take the charm quark mass such that it reproduces the physical value of the spin-averaged 1S mass of charmonium. The “bottom” quark masses are taken to be $(1.25)^{2}\simeq 1.56$ or $(1.25)^{4}\simeq 2.44$ times heavier than the charm. In principle, they will allow us to extrapolate to the physical bottom quark mass, which is 4.5 times heavier Chakraborty:2014aca . Given the exploratory nature of this work, we do not attempt to extrapolate.

The $D_{s}^{(*)}$ and $B_{s}$ meson masses calculated on the same ensembles but with higher statistics are listed in Table 2. They come from the study of $D$ and $B$ meson decay constants Fahy:2017enl .

We compute the four-point function with $t_{\mathrm{src}}$ = 0 and $t_{\mathrm{snk}}$ = 28, 42, or 56 on the $32^{3}\times 64$, $48^{3}\times 96$, and $64^{3}\times 128$ lattices, respectively. They roughly correspond to the same physical separation of $\sim$ 2.3 fm.

We calculate the four-point function as depicted in Figure 4. We first calculate the spectator $s$ quark propagator starting from a smeared source at $t_{\mathrm{src}}$ and switch to $b$ at $t_{\mathrm{snk}}$ using the sequential source method after applying the same smearing as the source. The $b$ propagator is terminated at $t_{2}$, and the sequential source method is applied again with the current $J_{\nu}$ ($V_{\nu}$ or $A_{\nu}$) to switch to $c$. We finally contract with the initial $b$ quark propagator, which is calculated with a point source at $t_{\mathrm{src}}$, with an insertion of another current $J_{\mu}$ ($V_{\mu}$ or $A_{\mu}$) at $t_{1}$.

The vector and axial-vector currents defined on the lattice are local currents constructed from the Möbius domain-wall fermion fields. They receive finite renormalization since the lattice currents are non-conserving. Precise chiral symmetry achieved with Möbius domain-wall fermion ensures that the same renormalization constant is shared by the vector and axial-vector currents. The renormalization constant is determined in the massless limit using the position-space method Tomii:2016xiv .

In this work, we concentrate on the zero-recoil limit, i.e. $\mathbf{q}$ = (0,0,0). With zero spatial momentum insertion, a symmetry emerges under an exchange of three spatial directions. We may average over the equivalent correlators under this symmetry, i.e. $C^{JJ}_{kk}(t)$ with $k$ = 1, 2 and 3. Also in this limit, the parity conservation forbids the crossing channels $C^{VA}_{\mu\nu}(t)$ and $C^{AV}_{\mu\nu}(t)$. (In the following we suppress the argument $\mathbf{q}=\mathbf{0}$ for brevity.)

6.2 Two-point and three-point functions

First of all, let us show the effective mass for the initial $B_{s}$ meson in Figure 5. They correspond to the two-point functions with smeared source and local sink (LS). We identify the plateau as shown in the plots by solid lines. From the data in the time interval $[t_{\mathrm{min}},t_{\mathrm{max}}]$ listed in Table 3, we extract the meson mass. The data for the smeared sink (SS) are noisier but gives a consistent result.

By combining the SL and SS correlators, we obtain the matrix element of the local pseudo-scalar density $a^{3}Z^{LL}\equiv(1/2M_{B})|\langle 0|P^{L}|B(\mathbf{0})\rangle|^{2}$. The numerical results for the masses and matrix elements are listed in Table 3. These correlators play the role of denominator in (16).

We also analyze the lattice data for a three-point function corresponding to the matrix element $\langle B|V_{0}|B\rangle$, which is normalized to 1 in the continuum theory, for which the vector current conserves. Since the current we used on the lattice is local, there exits a finite renormalization. Such finite renormalization factor $Z_{V}=Z_{A}$ on the same lattice ensemble was calculated in the massless limit Tomii:2016xiv . The matrix element $\langle B|V_{0}|B\rangle$ provides another way to determine the renormalization constant $Z_{V}$. Here, the current is made of heavy $b$ (and $\bar{b}$) quark field, and the discretization effect due to large $am_{b}$ could become significant.

We calculate a three-point function

[TABLE]

and divide it by the two-point function $C^{SS}(t_{\mathrm{snk}},t_{\mathrm{src}})$ to define the inverse renormalization constant $Z_{Vbb}^{-1}$.

We plot the ratio $C^{SVS}_{\mu\nu}(t_{\mathrm{snk}},t,t_{\mathrm{src}})/C^{SS}(t_{\mathrm{snk}},t_{\mathrm{src}})$ in Figure 6. We obtain a plateau for each $\beta$ and $b$ quark mass. When $am_{b}$ is small ($am_{b}\lesssim$ 0.4), the plateau starts from early time slices and resulting $Z_{Vbb}^{-1}$ is close to 1. For larger $am_{b}$ ($\sim$ 0.5–0.7) we find a larger deviation from 1 and that the plateau is narrower. This indicates larger discretization effects for finite $am_{b}$.

The results for $Z_{Vbb}^{-1}$ are listed in Table 3. The corresponding renormalization factor $Z_{V}^{-1}$ determined in the massless limit Tomii:2016xiv is 1.047(10), 1.038(6), 1.031(5) at $\beta$ = 4.17, 4.35, 4.47, respectively.

Such a strong effect due to finite quark mass can be partly understood within the framework of heavy quark effective theory implemented for lattice fermions ElKhadra:1996mp . For the domain-wall fermion, the inverse of the wave function normalization is given at the tree level as $(Z_{Vbb}^{(0)})^{-1}=1-(am_{b})^{2}/2-7(am_{b})^{4}/8-5(am_{b})^{6}/16+\cdots$, which shows significant deviation from 1 for large $am_{b}$.

In principle, the discretization effect as seen in the renormalization constant is irrelevant after taking the continuum limit. In the practical calculation, however, $am_{b}$ is not small enough to ignore the terms of $(am_{b})^{2}$ and higher, and some error will remain after the extrapolation with only three lattice spacings. To reduce the error as much as we can, we utilize the mass dependent renormalization factor $Z_{Vbb}^{-1}$ determined from the three-point function. For the $b\to c$ current, we combine the factor for the massless and massive currents as $\sqrt{Z_{V}Z_{Vbb}}$, assuming that the renormalization constant is mainly given by the wave function renormalization and the vertex correction plays only a minor role. Whether or not this prescription gives a good approximation can be checked by inspecting the continuum extrapolation of the final results.

6.3 Four-point functions

Now, let us show the main results for the four-point functions.

We set $t_{2}$ so that it is separated from $t_{\mathrm{snk}}$ by 8, 12, 16 at $\beta$ = 4.17, 4.35, 4.47 to allow for the ground state saturation of the final $B_{s}$ meson. According to the effective mass plots of the two-point functions, see Figure 5, the correlator may be slightly affected by excited state contributions at these time separations. But it is within two standard deviations, and is sufficient for this exploratory study.

The position of $t_{1}$ is then varied between 0 and $t_{2}$. The integral over $t_{1}$ to obtain the structure function uses the data with $t_{1}$ greater than the separation between $t_{2}$ and $t_{\mathrm{snk}}$.

The ratio (16) is plotted in Figure 7 for various channels. The data on the coarsest lattice ($\beta$ = 4.17 and $am_{b}$ = 0.68808) are shown in this plot. The $VV$ channel in the temporal direction ($\mu\nu=00$) and the $AA$ channel in the spatial direction ($\mu\nu=kk$) give the largest contributions, while the others (temporal $AA$ and spatial $VV$) are sub-dominant. The sign is flipped for the temporal channels as they acquire $i^{2}$ by the Wick rotation.

One can see that the ratio decays as $t_{1}$ is more separated from $t_{2}$. The spatial $AA$ channel corresponds to the $D_{s}^{*}$ meson in the final state, while the temporal $VV$ channel probes the $D$ meson. In the plot, we also draw lines showing the $D_{s}^{(*)}$ contribution using their masses separately measured and given in Table 2. We find that this ground state contribution actually explains the four-point function up to $t_{1}=t_{2}-1$, which implies that the ground state indeed gives the dominant contribution for these channels.

The other channels (temporal $AA$ and spatial $VV$) are smaller in magnitude and may have slightly steeper decays suggesting the contributions of excited states.

Similar plots are shown in Figures 8 and 9 for the data at $\beta$ = 4.35 and 4.47, respectively. The data at lighter $b$ quark $m_{b}=1.25^{2}m_{c}$ look very similar to those of the coarser lattice as they should be since the $b$ quark mass is matched. The results at heavier $b$ mass $am_{b}=1.25^{4}m_{c}$, on the other hand, show a prominent deviation from the ground state dominance of the temporal $VV$ and spatial $AA$ channels. Such a dependence on the initial $b$ quark mass is expected, since the similarity of the ground-state wave function between the initial ($B_{s}$) and final ($D_{s}^{(*)}$) states is worse. In other words, the excited state effects become more significant as heavy quark symmetry, a symmetry under the exchange of $c$ and $b$, is violated more strongly.

6.4 Structure functions

The structure function $T^{JJ}_{\mu\nu}(\omega)$ is obtained by integrating $C^{JJ}_{\mu\nu}(t_{2}-t_{1})$ over $t=t_{2}-t_{1}$ as defined in (18). The value of $\omega$ can be chosen arbitrarily in the region below the lowest hadronic energy state.

We first show the results for (18), which contains a contribution from the contact term at $t=0$. The contact term gives a constant contribution to $T^{JJ}_{\mu\nu}(\omega)$ without changing the overall shape of the function. For the final numerical results, we use the derivative form (19).

The integral in (18) is replaced by a sum over $t$ from 0 to $\infty$. We truncate the sum at $t_{1}=t_{\mathrm{snk}}-t_{2}$ to ensure the ground state saturation of the initial $B_{s}$ meson. For larger $t$, we assume that the correlator $C^{JJ}_{\mu\nu}(t)$ is dominated by the ground state $D_{s}^{(*)}$ and replace the correlator by the corresponding exponential function $\exp(-m_{D_{s}^{(*)}}t)$ to evaluate the rest of the summation. This assumption seems to be well satisfied for the two main channels (temporal $VV$ and spatial $AA$) as shown in Figures 7–9. For the mass $m_{D_{s}^{(*)}}$ we use the value measured from two-point functions, which is statistically more precise than the four-point function. When the lowest energy state is unknown as in the case for two other channels (temporal $AA$ and spatial $VV$), we need to rely on the exponential fall-off of the four-point functions shown in Figures 7–9. Because of the noisy signal at larger separations, $\omega$ would be limited to $\omega\approx 0$ such that the sum is well saturated by the short distance contributions.

The results for the four channels obtained on the coarsest lattice are shown in Figure 10. Note that the channels of $VA$ and those of $kl(k\not=l)$ vanish for zero spatial momentum.

As expected, in Figure 10 one finds that the dominant channels show an increase of $T^{JJ}_{\mu\nu}(\omega)$ toward the singularity at the physical $D_{s}$ (temporal $VV$) and $D_{s}^{*}$ (spatial $AA$) poles. In the plot, the pole location is shown by vertical dashed lines. For the channels without the corresponding single particle in the final state, the structure function is much smaller and no sign of singularity associated with $D_{s}^{(*)}$ is found.

The contribution from the $D_{s}^{(*)}$ meson pole, $1/(m_{D_{s}^{(*)}}-\omega)$, is shown in the plot (dotted curve) with an arbitrary normalization. As expected from the exponential fall-off of the four-point correlator, the temporal $VV$ and spatial $AA$ channels are well described by the pole contribution.

Figure 11 shows similar plots for the data at $\beta$ = 4.35. For the lighter $b$ quark (top), the structure functions look very similar to the corresponding plot at $\beta$ = 4.17. The result for heavier $b$ quark, on the other hand, shows some deviation from the pole dominance.

We now discuss the physics interpretation of the results. In order to avoid the contamination from the contact term, we define a normalized derivative of the structure functions:

[TABLE]

The prefactor $(m_{D^{(*)}}-\omega)^{2}$ is multiplied to make it dimensionless. It also plays the role of absorbing the leading singularity due to the ground state $D_{s}^{(*)}$ meson. $D^{JJ}_{\mu\nu}(\omega)$ approaches a constant $|h(1)|^{2}$ in the limit of $\omega\to m_{D_{s}^{(*)}}$ as (25) and (26) indicate. Therefore, in principle, one can extract the zero-recoil form factor $|h(1)|$ for the exclusive modes through $D^{JJ}_{\mu\nu}(\omega)$ in the limit $\omega\to m_{D_{s}^{(*)}}$. It would be statistically noisier than the standard calculation from three-point functions, though.

The results for a “$b$ quark” which is $1.25^{2}=1.56$ times heavier than charm is shown in Figure 12. The prefactor is $(m_{D_{s}}-\omega)^{2}$ or $(m_{D_{s}^{*}}-\omega)^{2}$ for the temporal $VV$ or spatial $AA$ channel, respectively. The results at three lattice spacings agree within their statistical errors, showing that the discretization effects are under good control.

We find that the function $D^{JJ}_{\mu\nu}(\omega)$ is nearly constant in the whole range of $\omega$ between 0 and $m_{D_{s}^{(*)}}$, which indicates the dominance of the ground state. This is natural, as one expects that the insertion of the temporal vector current $V_{0}$ does not disturb the meson when the initial and final quark masses are degenerate. Namely, charge conservation ensures $D^{VV}_{00}(\omega)=1$ in the limit of degenerate heavy quark masses. This relation is slightly violated for non-degenerate masses. For the axial-vector channel, this argument does not apply, but the heavy-quark spin symmetry between $D_{s}$ and $D_{s}^{(*)}$ explains the similarity of the two channels.

The results for a heavier “$b$ quark”, whose mass is $1.25^{4}=2.44$ times more than charm, are shown in Figure 13. Lattice data on the two fine lattices show a good agreement.

With a larger mass difference between charm and bottom, we expect more contributions from excited states since the heavy quark symmetry is more violated and the wave function overlaps less. The lattice results support this expectation, showing a slight increase of $|D^{JJ}_{\mu\nu}(\omega)|$ towards lower $\omega$. This implies the contribution from the excited states.

Our proposal is to use $|D^{JJ}_{\mu\nu}(\omega)|$ at small $\omega$ for the determination of $|V_{cb}|$. It is inclusive in the sense that all possible excited states contribute, and at the same time non-perturbative lattice calculation is possible. Numerical results at $\omega=0$ obtained in this work for the channels available are summarized in Table 4. As Figures 10 and 11 show, the channels $D^{VV}_{00}$ and $D^{AA}_{kk}$ that contain large contributions from the ground states are large $\sim O(1)$, while others are an order of magnitude smaller. This is not obvious from the tree-level analysis (34)–(35), in which all channels would be $O(1)$.

As a by-product of this analysis, we may also extract the zero-recoil form factors of $B_{s}\to D_{s}^{(*)}$ transitions from the limit $\omega\to m_{D_{s}^{(*)}}$ of $|D^{JJ}_{\mu\nu}(\omega)|$:

[TABLE]

We take the temporal $VV$ and spatial $AA$ channels for $h_{+}(1)$ and $h_{A_{1}}(1)$, respectively. Numerical results are listed in Table 5. These numbers are consistent with the recent lattice calculations $h_{+}(1)$ = 1.02–1.03 Bailey:2015rga and $h_{A_{1}}(1)$ = 0.906(4)(12) Bailey:2014tva , but have significantly larger errors. It should be noted that our results are obtained for lighter $b$ quark masses and the spectator quark is strange rather than up or down.

Even though the $b$ quark mass is different from its physical value, one can gain some information about the physical form factor using the $1/m$ expansion of the form factors, In the zero-recoil limit, it is written as

[TABLE]

with numerical coefficients calculated on quenched lattices as $c_{+}^{(2)}=(0.20(4)\mathrm{~{}GeV})^{2}$ and $c_{+}^{(3)}=(0.26(3)\mathrm{~{}GeV})^{3}$ Hashimoto:1999yp . For our parameters of $m_{b}/m_{c}$ = 2.45, the formula (39) gives 0.99(1), which is consistent with our results. A similar analysis for $|h_{A_{1}}(1)|$ with a quenched coefficients Hashimoto:2001nb gives $|h_{A_{1}}(1)|\approx$ 0.85, which again shows a good agreement.

6.5 Comparison to heavy quark expansion

Analytic calculation of the structure functions is possible using Heavy Quark Expansion (HQE), which is a form of the operator product expansion applied for heavy hadron decays. In the region away from the resonances, the perturbative calculation would be reliable. So far, the calculation is available only at the tree level but including the leading $1/m_{b}$ corrections. The leading order results are summarized in Section 5.2. They are essentially $\pm 1$ for $D^{JJ}_{\mu\nu}(\omega)$, and the $1/m_{b}$ corrections add some non-trivial $\omega$ dependence.

In order to compare the lattice data with analytic estimates, we plot the tree-level estimate including the $1/m_{b}$ corrections Blok:1993va in Figure 14. For the input parameters we took $\mu_{G}^{2}$ = 0.37 GeV2, $\mu_{\pi}^{2}$ = 0.5 GeV2, as well as $m_{c}$ = 1.5 GeV, as representative values. At this level of accuracy, i.e. zeroth order of $\alpha_{s}$, details of the input values are not very important. Results are however in good agreement near $\omega\approx 0$, and the HQE prediction rapidly diverges near the resonance region. It should be possible to extend this comparison to one-loop using existing calculations Trott:2004xc ; Aquila:2005hq . In order to do that we need to carry out the contour integral to convert the results for the decay rate to the structure functions at unphysical kinematics.

7 Conclusion and discussions

The inverse hadronic energy moments of the inclusive semi-leptonic $B$ meson decay rate can be calculated from the first principles of QCD using lattice QCD techniques. The formulation is based on the analytic continuation of the relevant structure functions from the imaginary cut to the unphysical kinematical region where the lattice QCD calculation can be applied. Performing an exploratory numerical computation, we have demonstrated that the strategy works well at least in the zero-recoil limit. The parameter $\omega$ controls the degree of inclusiveness. On the one end, one can extract the form factors of the exclusive modes, while the other end contains significant contributions from excited and continuum states.

In order to use this method for the determination of $|V_{cb}|$, the lattice calculation has to be extended toward the physical quark masses for both light and heavy quarks. Such calculation is technically straightforward, but computationally much more costly. Since the statistical signal would then become noisier, we need to investigate how much precision one can eventually achieve.

One interesting application in the zero-recoil momentum limit is to test the consistency with the heavy quark expansion approach that lead to the bound and estimates for $|h_{A_{1}}(1)|$, e.g. Gambino:2010bp ; Gambino:2012rd . These estimates appeared to be slightly lower than the lattice calculation and the source of the inconsistency is not understood so far. More direct comparison would be possible by calculating the structure function at the unphysical kinematics by the both methods. Since the lattice calculation can be performed with various heavy quark masses, it also provides a qualitative test of the heavy quark expansion approach.

An extension to the non-zero recoil momentum is another important direction required to compare with the experimental data. It will enable us to disentangle the individual structure functions $T_{i}$ from various components $T^{JJ}_{\mu\nu}$.

The experimental data taken at BaBar Aubert:2004aw and Belle Urquijo:2006wd ; Schwanda:2006nf contain the necessary information of the differential decay rate, but they need to be reanalyzed to extract the inverse hadron energy moment of the form (12) and (13) at each $q^{2}$ or $\mathbf{q}^{2}$. (The experimental analysis has been focused on the lepton energy moments and hadron invariant mass moments, so far.) We have not considered any effect due to experimental cuts and backgrounds that may arise in the moment analyses. These would be a main challenge on the experimental side.

The contribution from the unphysical (negative $q^{2}$) region needs to be estimated theoretically. It should be possible because this region is away from the resonances and may be treated by perturbation theory. One may also minimize the contribution from this region by considering the higher moments or other tailored moments. Such attempts has so far been made for the analysis of inclusive $\tau$ lepton decays to strange final state Maltman:2016ziu , which introduces a multiple pole function in the Cauchy integral. A similar idea should also work for the $B$ meson inclusive decay.

A natural extension of the strategy proposed in this paper would be to apply for the $b\to u$ process, for which the tension between the inclusive and exclusive $|V_{ub}|$ determination is more critical. (According to the summary in Particle Data Book (2016) Olive:2016xmw , the inclusive determination gives $|V_{ub}|=(4.49\pm 0.16^{+0.16}_{-0.18})\times 10^{-3}$, while the exclusive is $|V_{ub}|=(3.72\pm 0.19)\times 10^{-3}$.) In principle, the same technique is applicable by changing the final quark mass from $m_{c}$ to $m_{u}$, unless one approaches the region of $q^{2}\approx 0$ where the gap between the two cuts in the complex $v\cdot q$ plane becomes narrow. Since the statistical signal of lattice calculations is worse for large momentum in general, the small $q^{2}$ region would be more difficult also technically. On the experimental side, the $b\to u$ transition is more challenging because of the $b\to c$ backgrounds; the kinematical region where $b\to c$ could also occur is essentially not available. Whether or not there is a region that can be used for the inverse moment analysis proposed in this work needs to be investigated.

The rare decays $B\to X_{s}\ell^{+}\ell^{-}$ may be treated similarly, except that one has to avoid the $c\bar{c}$ resonance region in the invariant mass of $q^{2}$. The related process $B\to X_{s}\gamma$ would be difficult because of the same problem as the $q^{2}\approx 0$ region of the $b\to u$ transition.

The semi-leptonic decays of $D$ meson may play a role to validate the method. The heavy quark expansion would not work best for charm and theoretical calculation might not be reliable. The lattice calculation is, on the other hand, easier for charm than for bottom, as the conventional relativistic fermion formulation is applicable. The branching ratio of $D\to X_{s}\mu\nu$ is about 17%, among which the contributions of $K\mu\nu$ and $K\pi\mu\nu$ are 9% and 4%, respectively, and the bulk of $K\pi\mu\nu$ is actually $K^{*}(892)\mu\nu$. As there is still some room for other excited state contributions, it would become a good testing ground, provided that the experimental measurement of the differential decay rate is available. The situation for $B\to X_{c}\mu\nu$ is similar, i.e. the branching ratio for $B\to X_{c}\ell\nu$ (11%) is dominated by $D\ell\nu$ (2.3%) and $D^{*}\ell\nu$ (5.7%), and excited states $D\mathrm{n}\pi\ell\nu$ are small but still significant ($\sim$ 2%).

Apart from the $B$ meson decays, the method outlined in this work may be applied for inelastic scattering processes such as $e^{-}p^{+}\to e^{-}X$, for which the formulation of the structure functions was originally developed. The relevant structure functions, generally written as $W(Q^{2},\nu)$ with momentum transfer $q^{2}=-Q^{2}$ and electron energy loss $\nu$, may be analytically continued in the complex $\nu$ plane from a physical to unphysical kinematics at a fixed $Q^{2}$. In the kinematical setup where the energy of the final state $X$ is not sufficient to make it on-shell, the lattice calculation can be applied. If successful, it will allow us to analyze the inelastic processes without recourse to perturbation theory. In particular, relatively low-energy processes may be treated.

As an application of the lattice QCD calculation for particle phenomenology, the method proposed in this work opens new possibilities. By treating the inclusive processes, it probes relatively high energy processes compared to the standard applications, such as the ground state masses, decay constants and form factors. Their energy scale is, on the other hand, not sufficiently high to apply perturbation theory. Such intermediate energy region of QCD between perturbative and non-perturbative regimes has not been fully explored theoretically. One such attempt is the lattice calculation of the vacuum polarization function and its comparison with the $\tau$ decay experiment Tomii:2017xuk , and the analysis of Maltman:2016ziu is also in the similar direction. The present work extends these previous works from the simple two-point function to the semi-leptonic topology. The method will enable us to investigate how the quark-hadron duality works quantitatively for such delicate cases. The lattice QCD calculation provides non-perturbative calculation in the entire region covering both the perturbative and non-perturbative regimes.

Among other potential applications, the tension between the inclusive and exclusive determinations of $|V_{cb}|$ and $|V_{ub}|$ is the most important problem to be understood. Through the method discussed in this paper, it will become possible to test the various theoretical tools in the intermediate energy regimes by fully non-perturbative lattice calculations. Our goal is to finally reach the understanding of the problem and to obtain more precise determinations of $|V_{cb}|$ and $|V_{ub}|$.

Acknowledgments

I thank the members of the JLQCD collaboration for discussions and for providing the computational framework and lattice data. I am also grateful to Taku Izubuchi, Luchang Jin, Christoph Lehner, Zoltan Ligeti, Hiroshi Ohki, Amarjit Soni for useful discussions, and to Brian Colquhoun for carefully reading the manuscript. Numerical calculations are performed on Blue Gene/Q at KEK under its Large Scale Simulation Program (No. 16/17-14), on Oakforest-PACS supercomputer operated by Joint Center for Advanced High Performance Computing (JCAHPC), on COMA at Center for Computational Sciences, University of Tsukuba, and on Camphor 2 at Institute for Information Management and Communication, Kyoto University. This work is supported in part by JSPS KAKENHI Grant Number JP26247043 and by the Post-K supercomputer project through the Joint Institute for Computational Fundamental Science (JICFuS).