Inclusive semi-leptonic decays from lattice QCD
KKEK-CP-376
Inclusive semi-leptonic decays from lattice QCD
Paolo Gambino
Dipartimento di Fisica, Universit`a di Torino and INFN TorinoVia P. Giuria 1, I-10125, Torino, Italy
Shoji Hashimoto
Theory Center, Institute of Particle and Nuclear Studies,High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801, Japan andSchool of High Energy Accelerator Science, The Graduate Universityfor 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 non-perturbatively using lattice QCD simulations. The sum over all possible final states is achievedby a calculation of the forward-scattering matrix elements on the lattice, and the phase-space inte-gral is evaluated using their dependence on the time separation between two inserted currents. Weperform a pilot lattice computation for the ¯ B s → X c (cid:96) ¯ ν decay with an unphysical bottom quark massand compare the results with the corresponding OPE calculation. The method to treat the inclu-sive processes on the lattice can be applied to other processes, such as the lepton-nucleon inelasticscattering. Quark-hadron duality plays a key role in perturbativeQuantum Chromodynamics (QCD) calculations of phys-ical processes. It states that hadronic processes can becalculated taking quarks and gluons as final states, eventhough the actually observed final states are composedof hadrons. In order that the duality is satisfied, theprocesses must be summed or smeared over all possiblehadronic final states in some kinematical range [1], suchas a region of invariant mass squared, but it is not apriori known how large the smearing should be. A sys-tematic approach to duality is based on the OperatorProduct Expansion (OPE) [2, 3], which is constructedin the Euclidean domain and analytically continued tothe Minkowski domain. In the context of heavy quarkdecays the OPE is an expansion in inverse powers of theheavy quark mass, or more precisely of the energy release,and the analytic continuation entails an inevitable viola-tion of duality. While there are indications that dualityviolation plays a minor role in the analysis of inclusivesemi-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 achievedby non-perturbative methods. The current tension be-tween the inclusive and exclusive determinations of | V cb | [6, 7] makes any contribution in this direction timely.Lattice QCD simulation provides a means of non-perturbative QCD computation for various hadronic pro-cesses including heavy quark decays. It has been success-fully applied to the calculation of exclusive decay formfactors, which are essential for a precise determinationof the CKM elements, e.g. K → π(cid:96) ¯ ν , B → D ( ∗ ) (cid:96) ¯ ν , etc. (see [8] for their world averages), while the study of in-clusive processes is scarce, except for recent attempts toformulate methods to introduce an analytic continuation[9] or a smearing [10]. The inclusive processes are, on theother hand, difficult to treat on the lattice because theyconsist of many physical states often including multiple hadrons. To identify each amplitude and to sum over thephase space is nearly impossible due to the number ofstates involved. One may instead use analyticity and theoptical theorem to relate the total rate to another quan-tity that is calculable on the lattice. This approach hasbeen followed in simple cases, such as the e + e − → q ¯ q pro-cesses [11–14] and τ -lepton decays [15, 16], while the ap-plication to B meson semi-leptonic decays is much morecomplicated [9].In this work we develop a novel and general methodto compute the inclusive semi-leptonic decay rate on thelattice. The method is based on a technique to calcu-late smeared spectral density of hadron correlators [17](see also [10, 18] for a slightly different strategy). Theextraction of the spectral density ρ ( ω ) of hadronic cor-relation functions remains intractable, but once ρ ( ω ) issmeared over some energy range, one can construct agood approximation using the correlation functions cal-culated 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 com-putational resources are made available.In this paper we use the inclusive semi-leptonic decays¯ B s → X c (cid:96) ¯ ν to demonstrate how the method works. Here, X c stands for all possible charmed states which may oc-cur with the quark-level decay process b → c(cid:96) ¯ ν . Afterdescribing the kinematics of the decay and the methodto calculate the inclusive decay rate, we present a pilotlattice study.For the analysis of the ¯ B s → X c (cid:96) ¯ ν decay, we assigna momentum p µ to the initial ¯ B s meson, and momenta p µ(cid:96) and p µ ¯ ν to the leptons (cid:96) and ¯ ν in the final state, re-spectively. Thus, the hadronic state X c has momentum r µ = ( p − q ) µ with q µ = ( p (cid:96) + p ¯ ν ) µ . The differential decayrate is written as [19, 20] d Γ dq dq dE (cid:96) = G F | V cb | π L µν W µν , (1) a r X i v : . [ h e p - l a t ] M a y where G F is the Fermi constant. The momentum trans-fer q µ and the lepton energy E (cid:96) are evaluated in the restframe of the initial ¯ B s meson. The leptonic tensor L µν is explicitly written as L µν = p µ(cid:96) p ν ¯ ν − p (cid:96) · p ¯ ν g µν + p ν(cid:96) p µ ¯ ν − i(cid:15) µανβ p (cid:96),α p ¯ ν,β for massless neutrinos. The hadronic ten-sor W µν ( p, q ) is defined through W µν ( p, q ) = (cid:88) X c (2 π ) δ (4) ( p − q − r ) × E B s (cid:104) ¯ B s ( p ) | J µ † | X c ( r ) (cid:105)(cid:104) X c ( r ) | J ν | ¯ B s ( p ) (cid:105) . (2)It is summed over all possible final states X c to representthe inclusive decay. The electroweak current relevant forthis decay mode is J µ = ( V − A ) µ = ¯ cγ µ (1 − γ ) b .One can perform an integral over the lepton energy E (cid:96) in (1), and the remaining integrals over q and q canbe rewritten in terms of ω and q , energy and spatialmomentum squared of the final hadrons X c , respectively.Thus, the total decay rate can be calculated asΓ = G F | V cb | π (cid:90) q d q (cid:112) q (cid:88) l =0 ¯ X ( l ) , (3)where q = (( m B s − m D s ) / m B s ) and¯ X ( l ) ≡ (cid:90) m Bs − √ q √ m Ds + q dω X ( l ) (4)with X (0) = q ( W − W ii ) , (5) X (1) = − ( m B s − ω ) q k ( W k + W k ) , (6) X (2) = ( m B s − ω ) ( W kk + 2 W ii ) . (7)Here, we take the momentum q in the k -th direction,while the i -th direction is assumed to be perpendicularto that. The repeated indices in (5)–(7) are not summed.The integral with respect to ω in (4) represents the sumover states that could appear for a given momentum q .On the lattice, as a counterpart of the hadronic ten-sor W µν , one can calculate the forward-scattering matrixelements of the form [9] C JJµν ( t ; q ) = (cid:88) x e i q · x m B s (cid:104) ¯ B s ( ) | J † µ ( x ,t ) J ν ( , | ¯ B s ( ) (cid:105) (8)from four-point functions including the interpolating op-erators for the ¯ B s meson state | ¯ B s ( ) (cid:105) . Now we intro-duce the transfer matrix on the lattice e − ˆ Ht to expressthe time dependence of the matrix element in (8) as1 V m B s (cid:104) ¯ B s ( ) | ˜ J † µ ( − q ) e − ˆ Ht ˜ J ν ( q ) | ¯ B s ( ) (cid:105) , (9)where ˜ J ν ( q ) denotes a Fourier transform of the insertedcurrent: ˜ J ν ( q ) = (cid:80) x e i q · x J ν ( x ). On the other hand, the integral over ω in (4) can be rewritten in the form (cid:90) ∞ dω K ( ω, q ) (cid:104) ¯ B s ( ) | ˜ J † µ ( − q ) δ ( ˆ H − ω ) ˜ J ν ( q ) | ¯ B s ( ) (cid:105) = (cid:104) ¯ B s ( ) | ˜ J † µ ( − q ) K ( ˆ H, q ) ˜ J ν ( q ) | ¯ B s ( ) (cid:105) . (10)Here K ( ω, q ) represents an integral kernel determinedby the explicit form of the integrands (5)–(7). The ω -integral is implicit on the right hand side; all the inter-mediate states may exist between the currents. Compar-ing the right hand side with (9), we find that the integral(10) can be evaluated if the kernel operator is well ap-proximated by a polynomial of the form K ( ˆ H, q ) = k ( q ) + k ( q ) e − ˆ H + · · · + k N ( q ) e − N ˆ H (11)with some coefficients k j ( q ), since the matrix elementsof the individual term on the right hand side are nothingbut C JJµν ( t ; q )’s.The best approximation of K ( ˆ H, q ) can be achievedusing the Chebyshev polynomials. We define a state | ψ µ ( q ) (cid:105) on which the kernel operator is evaluated as | ψ µ ( q ) (cid:105) = e − ˆ Ht ˜ J µ ( q ) | ¯ B s ( ) (cid:105) . A small time evolution e − ˆ Ht with a constant time t is introduced to avoid anypotential divergence in (cid:104) ψ µ ( q ) | ψ ν ( q ) (cid:105) . We can then con-struct an approximation as (cid:104) ψ µ | K ( ˆ H ) | ψ ν (cid:105)(cid:104) ψ µ | ψ ν (cid:105) (cid:39) c ∗ N (cid:88) j =1 c ∗ j (cid:104) ψ µ | T ∗ j ( e − ˆ H ) | ψ ν (cid:105)(cid:104) ψ µ | ψ ν (cid:105) . (12)(The dependence on q is omitted for simplicity.) T ∗ j ( x )stands for the shifted Chebyshev polynomials, whichare derived from the standard Chebyshev polynomials T j ( x ) as T ∗ j ( x ) ≡ T j (2 x − ≤ x ≤
1. Their first few termsare T ∗ ( x ) = 1, T ∗ ( x ) = 2 x − T ∗ ( x ) = 8 x − x + 1, and the others can be obtained recursivelyby T ∗ j +1 ( x ) = (4 x − T ∗ j ( x ) − T ∗ j − ( x ). Each termof (cid:104) ψ µ | T ∗ j ( e − ˆ H ) | ψ ν (cid:105) / (cid:104) ψ µ | ψ ν (cid:105) can be constructed from C JJµν ( t + 2 t ) /C JJµν (2 t ) = (cid:104) ψ µ | e − ˆ Ht | ψ ν (cid:105) / (cid:104) ψ µ | ψ ν (cid:105) .The coefficients c ∗ j in (12) are obtained from c ∗ j = 2 π (cid:90) π dθ K (cid:18) − ln 1 + cos θ (cid:19) cos( jθ ) , (13)according to the general formula of the Chebyshev ap-proximation. The Chebyshev approximation is the best in the sense that its maximum deviation in x ∈ [0 ,
1] isminimized among polynomials of order N .The integral kernel K ( ω, q ) is chosen as K ( l ) σ ( ω ) = e ωt ( − (cid:112) q ) − l ( m B s − ω ) l × θ σ ( m B s − (cid:112) q − ω ) (14)for l = 0, 1, or 2 corresponding to X ( l ) , (5)–(7). An ap-proximate Heaviside step function θ σ ( x ) is introduced to . . . . . ω K ( ) σ ( ω ) σ = 0 . true N = 5 N = 10 N = 20 . . . . . ω K ( ) σ ( ω ) σ = 0 . true N = 5 N = 10 N = 20 . . . . . ω K ( ) σ ( ω ) σ = 0 . true N = 5 N = 10 N = 20 FIG. 1. Approximation of the weight function K ( l =0) σ ( ω ) withthe Chebyshev polynomials of e − ω . For each value of thesmearing width σ (= 0.2 (top), 0.1 (middle), 0.05 (bottom)),the approximations with the polynomial order N = 5 (dot-ted), 10 (dot-dashed), 20 (dashed) are plotted as well as thetrue curve (solid curve). realize the upper limit of the ω -integral. In order to sta-bilize the Chebyshev approximation, we smear the stepfunction over a small width σ . For an explicit form, wechose θ σ ( x ) = 1 / (1 + exp( − x/σ )). The extra factor e ωt in (14) cancels the short time evolution e − ˆ Ht in | ψ µ ( q ) (cid:105) .Fig. 1 demonstrates how well K ( l ) σ ( ω ) is approximatedwith certain orders of the polynomials, i.e. N = 5, 10and 20. An example for l = 0 is shown. Here we takethree representative values of σ : 0.2, 0.1 and 0.05 in lat-tice units. The comparison is made for parameters that roughly correspond to our lattice setup: the inverse lat-tice spacing 1 /a (cid:39) am B s (cid:39) . t /a = 1.The momentum insertion q is set to zero. The kernelfunction is well approximated with relatively low ordersof the polynomials, such as N = 10, when sufficientlysmeared, e.g. σ = 0.2. For smaller σ ’s, the function ex-hibits a more rapid change near the threshold ω = 1 . N = 20. Eventually wehave to take the limit σ →
0, and the error due to finite N has to be estimated. For l = 1 and 2 the polynomialapproximations are better than those for l = 0.We perform a pilot study of the method describedabove using lattice data computed on an ensemble with2+1 flavors of M¨obius domain-wall fermions (the ensem-ble “M- ud s a” in [21], which has 1 /a = 3.610(9) GeV).For the charm and bottom quarks in the valence sec-tor, the same lattice formulation is used. The charmquark 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, respec-tively. The bottom quark mass is taken as 2 . m c , whichis substantially smaller than the physical b quark mass.The corresponding B s meson mass is 3.45 GeV. In thissetup, the maximum possible spatial momentum in the B s → D s (cid:96) ¯ ν decay is ( m B s − m D s ) / m B s (cid:39) L × L t = 48 ×
96, and we calcu-late the forward-scattering matrix elements with spatialmomenta q of (0,0,0), (0,0,1), (0,0,2) and (0,0,3) in unitsof 2 π/La . The number of lattice configurations averagedis 100, and the measurement is performed with four dif-ferent source time-slices.For a fixed spatial momentum q , we compute a four-point function to extract C JJµν ( t ; q ) (more details of thelattice calculation are presented in [9]). We perform the ω -integral (4) using the representation (12). Matrix ele-ments of the shifted Chebyshev polynomials are obtainedfrom C JJµν ( t + 2 t ; q ) /C JJµν (2 t ; q ) at various t ’s (and t =1) by a fit with constraints |(cid:104) ψ µ | T ∗ j ( e − ˆ H ) | ψ ν (cid:105) / (cid:104) ψ µ | ψ ν (cid:105)| <
1, which is a necessary condition for the Chebyshev poly-nomials.First, we inspect how well the Chebyshev approxima-tion works by comparing the results for ¯ X (2) obtainedwith the polynomial order N = 5, 10, 15 at various val-ues of σ , the width of the smearing. Fig. 2 shows that thedependence on σ is mild and the limit of σ = 0 is alreadyreached at around σ = 0 .
05. The dependence on N isnot significant, which indicates that the approximationis already saturated at N (cid:39)
10. This is crucial becausethe error of the lattice data is too large to constrain thematrix elements (cid:104) ψ µ | T ∗ j ( e − ˆ H ) | ψ ν (cid:105) / (cid:104) ψ µ | ψ ν (cid:105) at j (cid:39)
10 orlarger. The results for ¯ X (0) and ¯ X (1) show the similartendency. We take σ = 0 .
05 in the following analysis; theresults are within statistical error even if we extrapolateto σ = 0.The lattice results for ¯ X = (cid:80) l =0 ¯ X ( l ) are comparedwith the OPE predictions in Fig. 3 as a function of q .Here, the results for different polarizations, i.e. longi-tudinal ( (cid:107) : µ , ν = 0 and 3) and perpendicular ( ⊥ : µ , .
00 0 .
05 0 .
10 0 .
15 0 . σ . . . . . ¯ X ( ) N = 5 N = 10 N = 15 FIG. 2. ¯ X (2) at q = 2 π/La (0 , ,
1) plotted as a function ofthe smearing width σ . Results with polynomial orders N =5, 10, 15 are shown. . . . . . . q [GeV ] ¯ X [ G e V ] X V V k X V V ⊥ X AA ⊥ X AA k FIG. 3. ¯ X as a function of q plotted in the physical unit.Longitudinal ( (cid:107) ) and perpendicular ( ⊥ ) polarizations are plot-ted for vector ( V V ) and axial-vector ( AA ) channels. Dot-ted and dashed curves show the lowest order and O (1 /m )OPE estimates for each channel of corresponding color, re-spectively. ν = 1 and 2) directions to q , are separately plotted forvector ( V V , squares) and axial-vector ( AA , circles) cur-rent contributions. The lowest order and O (1 /m ) OPEestimates [20] are shown in the same plot. The OPEpredictions are sensitive to the heavy quark masses. Wetake the MS mass for the charm quark, ¯ m c (3 GeV) =1.00 GeV, and the kinetic mass for the fictitious b quark, m kinb (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 ofthe semi-leptonic fit of [5], although they refer to a lightspectator and to the physical b mass. The dashed linesinclude O (1 /m ) power corrections, which are large andtend to improve the agreement with the lattice data com-pared to the free quark decay (dotted lines). . . . . . . q [GeV ] ¯ X p q [ G e V ] LO OPE O (1 /m ) OPE O (1 /m , α s ) OPElattice data
FIG. 4. Integrand of the q -integral plotted in the physicalunit. The dot-dashed curve is an interpolation of the latticedata, and the O (1 /m , α s ) OPE calculation is shown by thered curve. To obtain the total decay rate, we integrate ¯ X (cid:112) q over q as in (3). The vector and axial-vector contributions ofdifferent polarizations are added. The integrand is shownin Fig. 4. We fit ¯ X ( l ) / (cid:112) q − l by a polynomial of q tointerpolate the data points. The fit curve (dot-dashed) isterminated at q . We compare the lattice results withthe corresponding OPE prediction (red curve) including O (1 /m ) [23] and O ( α s ) [24] terms with α s = 0 .
27. Thepower corrections are controlled here by powers of thepartonic energy (cid:112) m c + q which ranges between 1 and1.5 GeV, significantly less than that for a physical b .They are singular at the partonic endpoint, where themaximum energy hits the mass-shell of charm quark andthe perturbative corrections show an integrable singular-ity.Integrating the fit to lattice data we obtain Γ / | V cb | =4 . × − GeV, where only the statistical error isshown. We note that the total decay rate is about fivetimes smaller than that of the physical B s meson, becauseof the smaller phase space for the artificially small b quarkmass. On the OPE side, several higher order correctionsare available for the total width, including the complete O ( α s ) [25, 26] and the O ( α s /m ) [27, 28] corrections. Weimplement them in the kinetic scheme using the same in-puts as above and obtain Γ / | V cb | = 5 . × − GeV.The dominant uncertainty is due to the value of the b quark mass, but missing higher order corrections and un-certainties on the matrix elements would also induce an O (10%) uncertainty. Despite these limitations, the agree-ment between the lattice and the OPE is remarkable.An immediate extension of this work is of course thecalculation of the inclusive semi-leptonic decay rate of B mesons and b baryons ( b → c(cid:96) ¯ ν and b → u(cid:96) ¯ ν ). Momentsof kinematical variables, such as the lepton energy mo-ments and hadronic invariant mass moments, can also becalculated by a slight modification of the method. A nu-merical challenge for the lattice calculation is the largerecoil momentum up to ∼ b → u transitions, the experimental analysis involvesvarious momentum cuts to veto unnecessary b → c back-grounds. Our method allows to apply arbitrary kine-matical cuts, and a fully non-perturbative calculation ispossible according to the experimental setup. A compar-ison to the OPE calculation at or closer to the physical b mass would provide a valuable test of the OPE, includingthe assumption of quark-hadron duality. It may also beused to determine the hadronic parameters appearing inthe heavy quark expansion. The fully non-perturbativelattice calculation can also be applied to D meson de-cays, for which the energy release is not sufficiently largeto yield reliable OPE calculations, and where one couldobserve the onset of quark-hadron duality.The possible applications of the framework are not lim-ited to heavy quark decays. Lepton-nucleon ( (cid:96)N ) scat-tering is another large area of application. Traditionally,it has been analyzed combining perturbation theory andnon-perturbative inputs, such as the parton distributionfunctions (PDFs). Instead, the method described in thiswork allows to directly compute the cross sections with- out recourse to intermediate quantities like PDFs, andit opens a new strategy to study the inelastic scatter-ings. Moreover, it will make it possible to perform non-perturbative calculation of low-energy scatterings, whichcannot be treated with the presently available techniques. ACKNOWLEDGMENTS
We thank the members of the JLQCD collaboration fordiscussions and for providing the computational frame-work and lattice data. Numerical calculations are per-formed on SX-Aurora TSUBASA at High Energy Accel-erator Research Organization (KEK) under its Particle,Nuclear and Astro Physics Simulation Program, as wellas on Oakforest-PACS supercomputer operated by JointCenter for Advanced High Performance Computing (JC-AHPC). This work is supported in part by JSPS KAK-ENHI Grant Number JP26247043 and by the Post-K andFugaku supercomputer project through the Joint Insti-tute for Computational Fundamental Science (JICFuS).PG is supported in part by the Italian Ministry of Re-search (MIUR) under grant PRIN 20172LNEEZ. [1] E. Poggio, H. R. Quinn, and S. Weinberg, Phys. Rev. D , 1958 (1976).[2] M. A. Shifman, in , Vol. 3 (World Scientific, Singapore, 2000)pp. 1447–1494, arXiv:hep-ph/0009131.[3] I. I. Bigi and N. Uraltsev, Int. J. Mod. Phys. A , 5201(2001), arXiv:hep-ph/0106346.[4] A. Alberti, P. Gambino, K. J. Healey, and S. Nandi,Phys. Rev. Lett. , 061802 (2015), arXiv:1411.6560[hep-ph].[5] P. Gambino, K. J. Healey, and S. Turczyk, Phys. Lett.B , 60 (2016), arXiv:1606.06174 [hep-ph].[6] M. Tanabashi et al. (Particle Data Group), Phys. Rev.D , 030001 (2018).[7] P. Gambino, M. Jung, and S. Schacht, Phys. Lett. B , 386 (2019), arXiv:1905.08209 [hep-ph].[8] S. Aoki et al. (Flavour Lattice Averaging Group), Eur.Phys. J. C , 113 (2020), arXiv:1902.08191 [hep-lat].[9] S. Hashimoto, PTEP , 053B03 (2017),arXiv:1703.01881 [hep-lat].[10] M. T. Hansen, H. B. Meyer, and D. Robaina, Phys. Rev.D , 094513 (2017), arXiv:1704.08993 [hep-lat].[11] D. Bernecker and H. B. Meyer, Eur. Phys. J. A , 148(2011), arXiv:1107.4388 [hep-lat].[12] X. Feng, S. Hashimoto, G. Hotzel, K. Jansen,M. Petschlies, and D. B. Renner, Phys. Rev. D ,034505 (2013), arXiv:1305.5878 [hep-lat].[13] A. Francis, B. Jaeger, H. B. Meyer, and H. Wittig, Phys.Rev. D , 054502 (2013), arXiv:1306.2532 [hep-lat].[14] C. Lehner and A. S. Meyer, Phys. Rev. D , 074515(2020), arXiv:2003.04177 [hep-lat]. [15] M. Tomii, G. Cossu, B. Fahy, H. Fukaya, S. Hashimoto,T. Kaneko, and J. Noaki (JLQCD), Phys. Rev. D ,054511 (2017), arXiv:1703.06249 [hep-lat].[16] P. Boyle, R. J. Hudspith, T. Izubuchi, A. Jttner,C. Lehner, R. Lewis, K. Maltman, H. Ohki, A. Portelli,and M. Spraggs (RBC, UKQCD), Phys. Rev. Lett. ,202003 (2018), arXiv:1803.07228 [hep-lat].[17] G. Bailas, S. Hashimoto, and T. Ishikawa, PTEP ,043B07 (2020), arXiv:2001.11779 [hep-lat].[18] J. Bulava and M. T. Hansen, Phys. Rev. D , 034521(2019), arXiv:1903.11735 [hep-lat].[19] A. V. Manohar and M. B. Wise, Phys. Rev. D , 1310(1994), arXiv:hep-ph/9308246.[20] B. Blok, L. Koyrakh, M. A. Shifman, and A. Vainshtein,Phys. Rev. D , 3356 (1994), [Erratum: Phys. Rev. D50, 3572 (1994)], arXiv:hep-ph/9307247.[21] K. Nakayama, B. Fahy, and S. Hashimoto, Phys. Rev.D , 054507 (2016), arXiv:1606.01002 [hep-lat].[22] P. Gambino, A. Melis, and S. Simula, Phys. Rev. D ,014511 (2017), arXiv:1704.06105 [hep-lat].[23] M. Gremm and A. Kapustin, Phys. Rev. D , 6924(1997), arXiv:hep-ph/9603448.[24] V. Aquila, P. Gambino, G. Ridolfi, and N. Uraltsev,Nucl. Phys. B , 77 (2005), arXiv:hep-ph/0503083.[25] A. Pak and A. Czarnecki, Phys. Rev. Lett. , 241807(2008), arXiv:0803.0960 [hep-ph].[26] K. Melnikov, Phys. Lett. B , 336 (2008),arXiv:0803.0951 [hep-ph].[27] A. Alberti, P. Gambino, and S. Nandi, JHEP , 147(2014), arXiv:1311.7381 [hep-ph].[28] T. Mannel, A. A. Pivovarov, and D. Rosenthal, Phys.Rev. D92