B c →J/ψ Form Factors for the full q 2 range from Lattice QCD
BB c → J/ψ
Form Factors for the full q range from Lattice QCD Judd Harrison,
1, a
Christine T. H. Davies,
1, b and Andrew Lytle (HPQCD Collaboration) c SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK INFN, Sezione di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Roma RM, Italy
We present the first lattice QCD determination of the B c → J/ψ vector and axial-vector formfactors. These will enable experimental information on the rate for B c semileptonic decays to J/ψ to be converted into a value for V cb . Our calculation covers the full physical q range of the decayand uses non-perturbatively renormalised lattice currents. We use the Highly Improved StaggeredQuark (HISQ) action for all valence quarks on the second generation MILC ensembles of gluonfield configurations including u , d , s and c HISQ sea quarks. Our HISQ heavy quarks have massesranging upwards from that of c ; we are able to reach that of the b on our finest lattices. This enablesus to map out the dependence on heavy quark mass and determine results in the continuum limitat the b . We use our form factors to construct the differential rates for B − c → J/ψµ − ν µ and obtaina total rate with 6% uncertainty: Γ( B − c → J/ψµ − ν µ ) / | η EW V cb | = 1 . × s − . Includingvalues for V cb , η EW and τ B c yields a branching fraction for this decay mode of 0.0151(10)(10)(3) withuncertainties from lattice QCD, η EW V cb and τ B c respectively. PACS numbers: 12.38.Gc, 13.20.Gd, 13.40.Hq, 14.40.Pq
I. INTRODUCTION
Accurate calculations of hadronic parameters areneeded for the determination of Cabibbo-Kobayashi-Maskawa (CKM) matrix elements from the compari-son of Standard Model theory and experimental resultsfor exclusive flavour-changing decay rates. Leptonic de-cay rates require decay constants to be determined andsemileptonic decay rates require form factors. LatticeQCD is now established as the method of choice for thecalculation of these hadronic parameters and efforts areongoing both to improve the accuracy of the results andto expand the range of processes for which calculationalresults are available. Here we report on the first lat-tice QCD calculation of the form factors for B c to J/ψ semileptonic decay, a process under active study by theLHCb experiment [1]. Because the valence quarks in-volved in this process are all heavy, the calculation isunder very good control in lattice QCD as we will show.This opens the prospect of a new exclusive process thatcan be used for the determination of | V cb | that has a re-duced theoretical uncertainty.Currently exclusive determinations of | V cb | are focussedon B → D ∗ and B → D semileptonic decays. Theemphasis is on the former pseudoscalar to vector mesontransition because of more favourable kinematic factorsfor the differential decay rate towards the zero recoil re-gion. Lattice QCD calculations are generally more accu-rate in this region because the daughter meson has smallspatial momentum in the lattice frame (where the par-ent meson is usually taken to be at rest). The emphasis a [email protected] b [email protected] c on this region led to the early lattice QCD B → D ∗ form factor calculations being purely done at zero recoil,where there is a single form factor [2–4]. Comparisonwas then made to results derived from experimental dif-ferential rates in this same limit. More recently (see, forexample, [5, 6]) it has become clear that extrapolatingthe experimental results to the zero recoil point comeswith significant systematic uncertainties, associated withthe underlying model dependence of such extrapolations,that were previously being underestimated. The way for-ward requires a much more direct comparison q -bin by q -bin of the experimental decay rate and that from the-ory, determined from lattice QCD form factors. For thiswe need lattice QCD form factors as a function of q and preferably covering the full q range of the decay.This allows the uncertainty in the extracted value of | V cb | to be optimised between the changing experimental andtheoretical uncertainties as a function of q . We showhere that it is possible to calculate the form factors for B c → J/ψ semileptonic decay across the full q range inlattice QCD.The B c → J/ψ form factor calculation that we de-scribe here acts as a prototype for upcoming calculationsof form factors for B → D ∗ and B s → D ∗ s . B c → J/ψ is amore attractive starting point for lattice QCD, however,because the mesons are ‘gold-plated’ (with tiny widths)and the valence quarks involved are all relatively heavy.This means that the valence quark propagators, fromwhich appropriate correlation functions are constructed,are inexpensive to calculate. The correlation functionsthen have small statistical errors, even when the daugh-ter has maximum spatial momentum. The absence ofvalence light quarks means that finite-volume effects arenegligible and sensitivity to the u/d quark mass in thesea should be small.The main obstacle for the calculation of B c → J/ψ form factors is that of the discretisation effects associ- a r X i v : . [ h e p - l a t ] J u l ated with the heavy quarks. The c quark is handled veryaccurately in lattice QCD as long as improved discreti-sations of the Dirac equation are used. A particularlyaccurate approach is that of HPQCD’s Highly ImprovedStaggered Quark (HISQ) action [7] and it is the one thatwe will use here for all quarks. For c quarks a recentcalculation using HISQ [8] obtained a 0.4% uncertaintyin the J/ψ decay constant and showed good control oflattice discretisation effects all the way to very coarselattices with a spacing of 0.15 fm. Discretisation effectsare larger for the b quark. Indeed, since we expect dis-cretisation effects to grow as a power of the heavy quarkmass in lattice units, am h , we need to work on fine lat-tices to approach the b quark mass, and this is what wewill do here. The dominant discretisation effects in theHISQ action behave as a , since tree-level a errors areremoved and those that include powers of α s are heavilysuppressed. This means that discretisation effects canbe controlled for quarks with masses m h between thatof the c and the b on lattices with lattice spacing below0.1 fm. By working with a range of quark masses reach-ing up to that of the b on a range of lattice spacings, wecan map out both the dependence on m h and the de-pendence on am h , both of which are smooth functions,and obtain a physical value in the continuum limit forthe b quark. This ‘heavy-HISQ’ approach was developedby HPQCD for B meson decay constants [9, 10] and isnow being extended to form factors [11]. Its efficacy hasbeen demonstrated for B s → D s form factors for the full q range in [12] and here we apply it to B c → J/ψ . Abig advantage of this approach is that the lattice currentoperators that couple to the W boson can be normalisedfully nonperturbatively, e.g.[12–15], avoiding the system-atic errors associated with the perturbative normalisa-tion needed if a nonrelativistic approach is used for the b quark in lattice QCD [3, 4].A further motivation for studying B c → J/ψ semilep-tonic decay in detail is to calculate the ratio, R ( J/ψ ), ofthe partial widths for the outgoing lepton to be a τ com-pared to that for it to be an e or µ . The analogous resultsfor B decays, R ( D ) and R ( D ∗ ), e.g.[16–19], have beena source of tension between experiment and the Stan-dard Model [20], implying lepton-universality violation.Recent results from Belle, on the contrary, show goodagreement with the Standard Model [21]. This makesit very important to test lepton-universality violation inother processes and we can obtain R ( J/ψ ) from our formfactors for comparison with ongoing LHCb analyses. Wewill present those results and analyses of other lepton-universality violation tests separately [22]; here we focuson the form factors and differential rates for W decay to µ or e .The subsequent sections are organised as follows: • In section II we begin by outlining the relevant ex-perimental observables and relate them to the in-variant form factors coming from QCD matrix ele-ments.
J/ µ + µ ` ⌫ ` B c ✓ W ✓ J/ W
Here we give the partial rates for B − c → J/ψ ( → µ + µ − ) (cid:96) − ν (cid:96) where (cid:96) is the final state lepton as differ-entials with respect to q and angular variables definedin the standard way in Fig. 1. In this work we consideronly the cases (cid:96) = µ and (cid:96) = e . The partial rates are ob-tained from the full differential decay rate assuming the J/ψ decay is purely electromagnetic and summing over µ + µ − helicities (assuming that the µ + µ − are masslessand hence pure helicity eigenstates) [23]. This gives d Γ dq = G F (2 π ) | η EW V cb | ( q − m (cid:96) ) | (cid:126)p (cid:48) | M B c q × (cid:104) (cid:0) H − + H + H +2 (cid:1) + m (cid:96) q (cid:0) H − + H + H +2 + 3 H t (cid:1)(cid:105) , (1) d Γ dq d cos( θ W ) = G F (2 π ) | η EW V cb | ( q − m (cid:96) ) | (cid:126)p (cid:48) | M B c q B ( J/ψ → µ + µ − ) (cid:104) H −
12 (1 + cos( θ W )) + H
12 (1 − cos( θ W )) + H sin ( θ W )+ m (cid:96) q (cid:16) (cid:0) H − + H (cid:1) sin ( θ W )+2( H t − H cos( θ W )) (cid:17)(cid:105) , (2) d Γ dq d cos( θ J/ψ ) = G F (2 π ) | η EW V cb | ( q − m (cid:96) ) | (cid:126)p (cid:48) | M B c q B ( J/ψ → µ + µ − ) (cid:104) ( H − + H ) 12 (cid:0) ( θ J/ψ ) (cid:1) + H (cid:0) − cos ( θ J/ψ ) (cid:1) + m (cid:96) q (cid:16) ( H − + H ) 12 (cid:0) ( θ J/ψ ) (cid:1) +( H + 3 H t ) (cid:0) − cos ( θ J/ψ ) (cid:1)(cid:17)(cid:105) , (3)and d Γ dq dχ = G F (2 π ) | η EW V cb | ( q − m (cid:96) ) | (cid:126)p (cid:48) | M B c q B ( J/ψ → µ + µ − ) 23 π (cid:104) H − + H + H + 12 H − H + cos(2 χ )+ m (cid:96) q (cid:16) H − + H + H + 3 H t − H − H + cos(2 χ ) (cid:17)(cid:105) . (4)Here p (cid:48) is the momentum of the J/ψ , p is the momen-tum of the B − c , q = p − p (cid:48) , | (cid:126)p (cid:48) | is the magnitude of the J/ψ spatial momentum in the B − c rest frame and m (cid:96) isthe mass of the final state lepton (cid:96) . G is the Fermi con-stant defined from muon lifetime, V cb is the appropriateCabibbo-Kobayashi-Maskawa matrix element and η EW is a structure-independent electroweak correction. Forthe electrically neutral J/ψ final state this is small [24],1 + α log( M Z /µ ) /π . Taking µ = M B c , the mass of the B − c meson, but allowing for variation by a factor of 2 [3],gives η EW = 1.0062(16). Note that we include in theexpressions for the differential rate terms with factors of m (cid:96) /q , which are negligible for (cid:96) = e but have a verysmall visible effect near q = 0 for (cid:96) = µ . The helicity amplitudes are defined as H ± ( q ) =( M B c + M J/ψ ) A ( q ) ∓ M B c | (cid:126)p (cid:48) | M B c + M J/ψ V ( q ) ,H ( q ) = 12 M J/ψ (cid:112) q (cid:16) − M B c (cid:126)p (cid:48) M B c + M J/ψ A ( q )+ ( M B c + M J/ψ )( M B c − M J/ψ − q ) A ( q ) (cid:17) ,H t ( q ) = 2 M B c | (cid:126)p (cid:48) | (cid:112) q A ( q ) (5)and correspond to the nonzero values of (cid:15) ∗ µ ( q, λ (cid:48) ) (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ cγ µ (1 − γ ) b | B − c ( p ) (cid:105) for the dif-ferent combinations of the J/ψ and W − polarisations λ and λ (cid:48) respectively. The form factors in Eq. (5) are thestandard Lorentz invariant ones, their relations to thematrix elements are given by [25] (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ cγ µ b | B − c ( p ) (cid:105) =2 iV ( q ) M B c + M J/ψ ε µνρσ (cid:15) ∗ ν ( p (cid:48) , λ ) p (cid:48) ρ p σ (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ cγ µ γ b | B − c ( p ) (cid:105) =2 M J/ψ A ( q ) (cid:15) ∗ ( p (cid:48) , λ ) · qq q µ + ( M B c + M J/ψ ) A ( q ) (cid:104) (cid:15) ∗ µ ( p (cid:48) , λ ) − (cid:15) ∗ ( p (cid:48) , λ ) · qq q µ (cid:105) − A ( q ) (cid:15) ∗ ( p (cid:48) , λ ) · qM B c + M J/ψ (cid:104) p µ + p (cid:48) µ − M B c − M J/ψ q q µ (cid:105) . (6)We also have (cid:104) | ¯ cγ ν c | J/ψ ( p (cid:48) , λ ) (cid:105) = N J/ψ (cid:15) ν ( p (cid:48) , λ ) , (7) (cid:104) B − c ( p ) | ¯ bγ c | (cid:105) = N B c , (8)and (cid:88) λ (cid:15) ν ( p (cid:48) , λ ) (cid:15) ∗ µ ( p (cid:48) , λ ) = − g νµ + p (cid:48) ν p (cid:48) µ M (9)where N J/ψ and N B c are amplitudes proportional to thedecay constant of the corresponding meson and (cid:15) is the J/ψ polarisation vector. We will make use of these whenwe come to extract the form factors in Eq. (6) from ourlattice correlation functions.
III. COMPUTATIONAL STRATEGY
We follow the strategy of previous heavy-HISQ calcula-tions of form-factors [11, 12, 15]. We work with multipleheavy quarks with masses m h between the physical c and b masses on lattices with a range of fine lattice spacings.Most of our m h masses are below the b mass but we areable to reach a mass very close to the physical b quarkmass on our finest lattices. An important point is that,at every lattice spacing, we are able to cover the full q range of the decay for the heavy quark masses thatwe can reach at that lattice spacing, i.e. the q rangeexpands on finer lattices in step with the heavy quarkmass range [12]. We compute form factors by extract-ing combinations of the relevant matrix elements definedin Eq. (6) from correlation functions computed on thelattice. We then fit the form factors as a function oflattice spacing and heavy quark mass to determine theirfunctional form in the continuum limit at the physical b mass.The correlation functions that we calculate, for generalchoices ν and Γ of J/ψ polarisation and current operatorrespectively, are C J/ψ ( t,
0) = (cid:104) | ¯ cγ ν c ( t ) (¯ cγ ν c (0)) † | (cid:105) ,C H c ( t,
0) = (cid:104) | (cid:0) ¯ hγ c ( t ) (cid:1) † ¯ hγ c (0) | (cid:105) ,C ( T, t,
0) = (cid:104) | ¯ cγ ν c ( T ) ¯ c Γ h ( t ) ¯ hγ c (0) | (cid:105) . (10)Note that in this section, for notational simplicity, weconsider the matrix elements in terms of continuum cur-rent operators. The nonperturbative renormalisation ofour lattice current operators is discussed in Section IV A.By considering the insertion of complete sets of stateswe may express these correlation functions in terms ofgeneral fit forms (folding the two-point functions aboutthe mid-point of the lattice so that t extends from 0 to N t / C J/ψ ( t,
0) = (cid:88) n (cid:16) ( A n ) e − tE n + ( − t ( A no ) e − tE on (cid:17) C H c ( t,
0) = (cid:88) n (cid:16) ( B n ) e − tM n + ( − t ( B no ) e − tM on (cid:17) and C ( T, t,
0) = (cid:88) n,m (cid:16) A n B m J nm e − ( T − t ) E n − tM m +( − T + t A no B m J nmoe e − ( T − t ) E on − tM m +( − t A n B mo J nmeo e − ( T − t ) E n − tM om +( − T A no B mo J nmoo e − ( T − t ) E on − tM om (cid:17) . (11)Here n , m correspond to on shell particle states withquantum numbers resulting in nonzero amplitudes andthe o labels indicate energies and amplitudes corre-sponding to the time-doubled states typical of staggeredquarks. The lowest energy, non-oscillating states arethose corresponding to the J/ψ and H − c . We work withthe H − c at rest, choosing to access the full q range bygiving momentum only to the J/ψ , and extract the ma-trix elements of these states from our lattice correlators. This gives: A = N J/ψ (cid:112) E J/ψ (cid:32) (cid:126)p (cid:48) ν M J/ψ (cid:33) / , (12) B = N H c (cid:112) M H c (13)and J ν, Γ) = (cid:88) λ (cid:15) ν ( p (cid:48) , λ ) (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ c Γ b | H − c (cid:105) (cid:114) E J/ψ M H c (cid:16) (cid:126)p (cid:48) ν /M J/ψ (cid:17) (14)where (cid:126)p (cid:48) ν is the ν component of the J/ψ spatial momen-tum, with ν corresponding to the choice of polarisationin Eq. (10), with current c Γ h . In the subsequent subsec-tions we discuss the combinations of ν and Γ for whichwe must compute correlation functions in order to extractthe full set of correlation functions defined in Eq. (6). A. Extracting V ( q ) The choice of operators used to extract V ( q ) is themost straightforward since the matrix element of the vec-tor current cγ µ b involves only V . With p = ( M H c , , , (cid:88) λ (cid:15) ν ( p (cid:48) , λ ) (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ cγ µ h | H − c (cid:105) = (cid:88) λ iV ( q ) M H c + M J/ψ ε µκρσ (cid:15) ν ( p (cid:48) , λ ) (cid:15) ∗ κ ( p (cid:48) , λ ) p (cid:48) ρ p σ = − i(cid:15) µνρ p (cid:48) ρ M H c M H c + M J/ψ V ( q ) (15)In this calculation we give the J/ψ spatial momentum (cid:126)p (cid:48) = ( k, k, (cid:126)p (cid:48) to be zero. Keeping both ofthe others non-zero minmises the discretisation errors fora given magnitude of p (cid:48) . Here we choose µ = 3 and ν = 1and find V ( q ) =Φ( k ) M H c + M J/ψ ikM H c J ,γ ) (16)where we have defined the relativistic normalisationΦ( k ) = (cid:114) E J/ψ M H c (cid:16) k /M J/ψ (cid:17) (17)with k the ν component of p (cid:48) . B. Extracting A ( q ) In order to isolate A ( q ), following [26], we make useof the partially conserved axial-vector current (PCAC)relation (cid:104) ∂A (cid:105) = ( m c + m h ) (cid:104) P (cid:105) where A = cγ γ ν h and P = cγ h . From Eq. (6) we have (cid:88) λ (cid:15) ν ( p (cid:48) , λ ) q µ (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ cγ µ γ h | H − c (cid:105) = (cid:88) λ (cid:15) ν ( p (cid:48) , λ )2 M J/ψ A ( q ) (cid:15) ∗ ( p (cid:48) , λ ) · q = 2 kE J/ψ M H c M J/ψ A ( q ) . (18)Taking Γ µ = γ and ν = 1 in Eq. (14) and multiplyingby m c + m b we then have A ( q ) = Φ( k ) ( m c + m b ) M J/ψ ikE J/ψ M H c J ,γ ) (19) C. Extracting A ( q ) In order to isolate A we use the axial-vector currentΓ = ¯ cγ µ γ h and J/ψ operator ¯ cγ ν c and choose µ = ν = 3along the spatial direction with zero J/ψ momentum.Using Eq. (6) this gives (cid:88) λ (cid:15) ( p (cid:48) , λ ) (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ cγ γ h | H − c (cid:105) =( M J/ψ + M H c ) A ( q ) (20)which gives, in terms of J A ( q ) = Φ(0) J ,γ γ ) M J/ψ + M H c . (21) D. Extracting A ( q ) The extraction of A is more complicated than theextraction of the other form factors since no trivial choiceof directions in axial-vector and J/ψ operators isolatesthe contribution of A relative to A or A . We use axial-vector current operator Γ = ¯ cγ γ h and J/ψ operator¯ cγ c . This yields contributions from each form factor,Φ( k ) J ,γ γ ) = (cid:88) λ (cid:15) ( p (cid:48) , λ ) (cid:104) J/ψ ( p (cid:48) , λ ) | ¯ cγ γ b | H − c (cid:105) = − k E J/ψ M H c q M J/ψ A ( q )+ ( M H c + M J/ψ ) (cid:32) k M J/ψ + E J/ψ M H c k M J/ψ q (cid:33) A ( q ) − A ( q ) k E J/ψ M H c M J/ψ ( M H c + M J/ψ ) (cid:32) M H c − M J/ψ q (cid:33) . (22)We must then subtract the A and A contributions toobtain A . TABLE I. Details of the gauge field configurations used inour calculation [27, 28]. Set 1 is referred to as ‘fine’, set 2 as‘superfine’, set 3 as ‘ultrafine’ and set 4 as ‘physical fine’. Thephysical value of ω was determined in [29] to be 0.1715(9)fmand the values of ω /a were taken from [12, 30, 31].Set ω /a N x × N t am l am s am c n configs . ×
96 0 . .
037 0 .
440 9802 2 . ×
144 0 . .
024 0 .
286 5003 3 . ×
192 0 . . .
188 3744 1 . ×
96 0 . . .
188 300 Tt J J/ H c h c c FIG. 2. Arrangement of propagators in the three point func-tion, we refer to c and c as the ‘spectator’ and ‘active’ charmquarks respectively and h as the ‘extended’ heavy quark. J represents the insertion of a vector or axial-vector current. IV. LATTICE CALCULATION
In this section we give details of the gauge field con-figurations used in this calculation, as well as the valuesof masses, momenta and staggered spin-taste operatorsused in the calculation. We use the the second genera-tion MILC gluon ensembles including light, strange andcharm sea quarks [27, 28] and compute HISQ charm andheavy quark propagators. The details of these gauge con-figurations are given in table I. As stated earlier, we workwith the H c at rest on the lattice. The arrangement of op-erators in the three point correlation functions is shownin figure 2. We refer to c as the ‘spectator’ charm quark, c as the ‘active’ quark and h as the ‘extended’ heavyquark. We use a single random wall source at time T for both the spectator and active quark, with a phasepatterning the source for the active quark to achieve anappropriate staggered spin-taste for the meson. We usethe spectator propagator at time 0 as the source for theextended heavy quark propagator with a spin-taste oper-ator corresponding to the H c quantum numbers and thentie the heavy quark propagator together with the activecharm quark propagator at the current at time t .Following the notation of [7] for staggered spin-tasteoperators, we choose lattice meson and current operatorswhich have the same quantum numbers as those detailedin section III. We always take the currents to be the lo-cal ones, γ i ⊗ γ i and γ i γ ⊗ γ i γ . This means that thecurrent insertions only need to be renormalised once. In TABLE II. Spin-taste operators used to isolate form factors.The first column is the operator used for the H c , the secondfor the J/ψ and the third column is the operator used at thecurrent. O H c O J/ψ O J V γ γ ⊗ γ γ γ ⊗ γ γ γ ⊗ γ A γ ⊗ γ γ ⊗ γ ⊗ γ A γ ⊗ γ γ ⊗ γ γ γ ⊗ γ γ A γ ⊗ γ γ ⊗ γ γ γ ⊗ γ γ TABLE III. Details of the charm and heavy valence masses.Set am val h am val c . , . , . . . , . , . , . . . , . , . . . , . , . . order to make the desired matrix element an overall tastesinglet, this requires that we use different tastes of J/ψ and H c . The specific combinations we use are given intable II. This means that we have three different tastesof J/ψ and two different tastes of H c . Different tastes ofa meson differ in mass by discretisation effects that aresmall for HISQ quarks [7, 27] and particularly small forheavy mesons. When fitting two- and three-point corre-lators the energies corresponding to different spin-tasteoperators use separate priors (i.e. we assume that theyare different) except for γ ⊗ γ and γ ⊗ γ whose ener-gies we fix to be equal. We will show later how close thedifferent taste meson masses are.We also calculate pseudoscalar heavyonium two-pointcorrelation functions with spin-taste γ ⊗ γ . This allowsus to determine the mass of the ground-state heavyoniummeson that we denote η h . The mass of the η h is a usefulphysical proxy for the heavy quark mass when we mapout the dependence on heavy quark mass of our results.We choose values of am h for the heavy quark massesat each value of the lattice spacing following [11]. Theyrange from above the c mass up to a value of 0.8. Thiscorresponds to a growing range in m h as the lattice spac-ing becomes finer. The valence masses used to computepropagators are given in table III. We also choose thevalues of T for our 3-point correlation functions follow-ing [11].Since we must compute heavy quark propagators atmultiple masses, from multiple T , the choice to insertmomentum via the active charm quark minimises thenumber of inversions. We put momentum on the c quarkpropagator via a momentum twist [32, 33]. The twistswere chosen to span evenly the physical q range for thelargest value of am h , approximated from aM H c valuesgiven in [11] and the physical J/ψ mass. The values oftwists and T are given in table IV. Note that because we TABLE IV. Details of the twists and three-point ranges, T ,used in our calculation. The twists are given in units of π/L and are applied along both the x and y direction. The mo-mentum component of the J/ψ in each of these directions inlattice units is then πθ/N x .Set θ T , . , . , . , . , .
807 14 , ,
202 0 , . , . , . , . , .
128 22 , ,
283 0 , . , . , . , . , .
207 31 , ,
414 0 , . , . , . , . , .
807 14 , , choose twists to evenly span the physical q range for themaximum value of am val h on a given lattice, the greatestvalues of twist, θ , will correspond to negative values of q for the smaller values of am val h used. We include thesenegative q points in our analysis, since the form factorsare analytic for −∞ < q < M minpole where M minpole is thelowest mass subthreshold resonance. A. Non-Perturbative Current Renormalisation
The renormalisation factors relating our lattice cur-rent operators to the partially conserved lattice currentswere computed in [11] for the axial-vector and in [12]for the vector current using the PCAC and PCVC re-lations respectively. These factors were not computedfor am h = 0 . amh = 0 .
65 on set 4 andfor these we interpolate between the other values givenin [11, 12], noting that the differences between results onset 1 and set 4 are very small. We tabulate the renor-malisation factors in table V, where we also include the ma -dependent discretisation correction terms, Z disc , forthe HISQ-quark tree level wavefunction renormalisationcomputed in [34]. Z disc starts at ( am h ) because of thehigh level of improvement in the HISQ action, so is veryclose to 1. The total renormalisation factor is given by Z V,A Z disc . Note that the pseudoscalar current used todetermine A is absolutely normalised and needs no Z factor.The statistical uncertainty of each lattice matrix ele-ment is typically at least an order of magnitude greaterthan that of the corresponding current renormalisationfactor. We therefore ignore correlations between differ-ent Z factors and between Z factors and our correlatordata. V. RESULTS
The correlator fits discussed in Section III were donesimultaneously for all correlators on each set using the corrfitter python package [35]. This allows us to deter-mine the ground-state fit parameters that we need alongwith a complete covariance matrix for them that can be
TABLE V. Z factors from [11] and [12] for the axial-vector andvector operators used in this work, together with the discreti-sation corrections, Z disc . Z A and Z V values for am h = 0 . am h = 0 .
65 on set 4 were obtained by interpolationfrom the other values for those sets.Set am h Z A Z V Z disc E i = E i +1 − E i , i (cid:54) = 0,Ω J/ψ = (3 . + p (cid:48) ) and Ω H c = M max H c (cid:0) am h . (cid:1) where M max H c is the value of M H c corresponding to the largest am h , takenfrom [11]. While Ω J/ψ was chosen to follow the relativisticdispersion relation, Ω H c was chosen heuristically to give priorvalues approximately following the observed H c masses oneach set while remaining suitably loose so as not to constrainthe fit results.Prior η b η c J/ψ ( p (cid:48) ) H c E / GeV m h . .
5) 3 . .
6) Ω
J/ψ . .
2) Ω H c . . E i / GeV m h . .
9) 1 . .
25) Ω
J/ψ . .
75) Ω H c . . E o / GeV − − Ω J/ψ . .
5) Ω H c . . E oi / GeV − − Ω J/ψ . .
75) Ω H c . . A ( B ) n ( o ) . .
0) 0 . .
0) 0 . .
0) 1 . . fed into our continuum extrapolation fits. The priorsfor our correlator fits are listed in table VI and take thesame form on all sets. Our correlator fits include an SVDcut on the correlation matrix following the tools includedin [35] (see Appendix D of [36] for a discussion of this).We also truncate the time-length of the three point andtwo point correlators that we fit to remove excited-statedominated data points at early times and increase fitstability. The specifics of these fit choices are given intable VII, along with variations used to test the stabilityof our analysis in Section V B.As discussed in section III we isolate form factors fromthe matrix elements extracted from correlation functionsand combine these with current renormalisation factorsto obtain A , A , A and V across a range of q forseveral different heavy quark masses. The form factor TABLE VII. Details of fit parameters, together with varia-tions used in section V B to check stability. ∆ T indicates thenumber of data points at the extremities of correlation func-tions not included in the fit. Bold values are those used toproduce our final results. χ / dof is estimated by introducingsvd and prior noise as in [35]. We do not compute χ valuesincluding prior and svd noise for those fits with n exp = 4.Set n exp ∆ T ∆ T J/ψ ∆ T H c SVD cut χ / dof δ
03 2 5 5 0.05 0.99 13 3 7 7 0.025 0.96 24 2 5 5 0.025 −
03 3 6 6 0.1 0.98 13 4 9 9 0.05 1.0 24 3 6 6 0.05 −
03 2 5 5 0.1 0.94 13 3 7 7 0.025 0.97 24 3 6 6 0.025 −
03 2 5 5 0.1 0.97 13 3 7 7 0.075 0.97 24 2 5 5 0.075 − ak here isthe value of the x and y components of the lattice momentumfor the J/ψ . ak is calculated from the corresponding twist inTable IV. am h ak A A A V − − − − − − − − − TABLE IX. Lattice form factor results for set 2. ak here isthe value of the x and y components of the lattice momentumfor the J/ψ . ak is calculated from the corresponding twist inTable IV. am h ak A A A V − − − − − − − − − − − − results are given in Tables VIII, IX, X and XI and plottedin Figure 3. We see that the least well-determined formfactor is A , especially close to zero recoil. This can betraced back to the division by k needed to determineit after the subtraction of the O (1) A contribution inEq. (22). In Figure 4 we plot A ( q ) | (cid:126)p (cid:48) | /M J/ψ , whichis the combination of A and | (cid:126)p (cid:48) | entering the helicityamplitudes in Eq. (5). We see that this has much clearerbehaviour as a function of q .We also give values for the meson masses in tables XIIfor the J/ψ and table XIII for the H c and η h . Fromtable II we have 3 different J/ψ tastes, using one localand two different 1-link point-split operators. The differ-ence in mass of these different taste
J/ψ mesons is tiny,barely visible above statistical uncertainties (at a levelof 1 MeV) even on our coarsest lattices. For the H c weuse two spin-taste operators, both local. As expectedfor pseudoscalars the taste-differences are slightly largerthan for the vector J/ψ and they are visible above thesmaller statistical errors in this case. The mass differ-
TABLE X. Lattice form factor results for set 3. ak here is thevalue of the x and y components of the lattice momentumfor the J/ψ . ak is calculated from the corresponding twist inTable IV. am h ak A A A V − − − − − − − − − ak here isthe value of the x and y components of the lattice momentumfor the J/ψ . ak is calculated from the corresponding twist inTable IV. am h ak A A A V − − − − − − − − − TABLE XII.
J/ψ masses for the local spin-taste operator γ ⊗ γ and 1 − link operators γ ⊗ γ ⊗ γ γ used in ourcalculation, see table II. aM J/ψ
Set γ ⊗ γ γ ⊗ γ ⊗ γ γ η h masses and H c masses for the local spin-taste operators γ ⊗ γ and γ γ ⊗ γ γ that we use in ourcalculation, see table II.Set am h aM H c ( γ ⊗ γ ) aM H c ( γ γ ⊗ γ γ ) aM η h ences are still only a few MeV, however, for a meson ofa mass of several GeV, and they fall on the finer latticessince they are a discretisation effect.We calculate the value of q for each form-factor basedon the meson masses of the corresponding meson tastes,along with the lattice momentum of the J/ψ determinedfrom the twists in Table IV.
A. Extrapolation to the Physical Point
In order to extract physical continuum form factorswe must fit our lattice results to a fit function includingdiscretisation effects arising from the use of large valuesof am h as well as the physical dependence upon M H c and q . In order to fit the q dependence we make use of the z -expansion (see, for example, [46–48]) which has becomestandard for the description of form factors. We map thephysical q range to within the unit circle via the changeof variables z ( q , t ) = (cid:112) t + − q − √ t + − t (cid:112) t + − q + √ t + − t . (23) TABLE XIV. Values used in our fits for the physical massesof relevant mesons, in GeV. These are from the Particle DataGroup [37] except for the unphysical η s which we take fromlattice calculations of the mass of the pion and kaon [38].The η s mass is used to set mass mistuning terms in our fitand so include an uncertainty. The other masses are usedas kinematic parameters in setting up our fit in z -space andused without uncertainties. The J/ψ mass is also used totune the c quark mass, and there we do include its (negligible)uncertainty. meson M phys [GeV] η b B c B J/ψ D ∗ η s B c pseudoscalar, vector and axialvector masses below BD ∗ threshold that we use in our polefactor, Eq. (26). Pseudoscalar values for the ground-state andfirst radial excitation are taken from experiment [37, 39–41];the other values are taken from [4] and are derived from latticeQCD calculations [42] and model estimates [43–45].0 − / GeV 1 − / GeV 1 + / GeV6.275 6.335 6.7456.872 6.926 6.757.25 7.02 7.157.28 7.15
The physical q range is from 0 to t − where t − = ( M H c − M J/ψ ) . (24)In Eq. 23 t + is the pair production threshold for a bc current, and the start of a cut in the complex q plane.This point is mapped to z = −
1. We take t + = ( M H + M D ∗ ) . (25)where H is an h ¯ u pseudoscalar meson. Note there is asizeable gap between t / and t / − and it is independentof m h up to binding energy effects. In this calculationwe do not compute M H entering t + . Instead we use anestimate of it derived from the M H c values that we have: M latt H = M H c − ( M phys B c − M phys B ). This ensures thatwhen we evaluate our fit function at the physical heavyquark mass then the correct physical continuum value of t + is recovered. We take the D ∗ mass from experimentsince our valence c masses are tuned to the physical valuethroughout and the mistuning of the light quark massesin the sea is taken care of in our fit function.Table XIV lists the physical masses we use for the B c , B and D ∗ . These are used as numerical constants with0 − q / GeV . . . . . . . A a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 . − q / GeV . . . . . . . . . A a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 . − q / GeV − − A a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 . − q / GeV . . . . . . V a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 . FIG. 3. The points show our lattice QCD results for each form factor as given in Tables VIII, IX, X and XI as a function ofsquared 4-momentum transfer, q . The legend gives the mapping between symbol colour and shape and the set of gluon fieldconfigurations used, as given by the lattice spacing, and the heavy quark in lattice units. The blue curve with error band isthe result of our fit in lattice spacing and heavy quark mass, evaluating the result in the continuum limit and for the b quarkmass, to give the physical form factor for B c → J/ψ . no uncertainty to simplify the necessary covariances re-quired to reconstruct our fit function. Their uncertaintiesare negligibly small in any case.The choice of t in Eq. (23) determines the value of q at which z = 0. Here we choose t = t − . This choicemeans that z is simply related to the variable w which isthe dot product of the 4-velocities of the B c and J/ψ andknown to be a good variable for a description of b → c form factors in the Heavy Quark Effective Theory ap-proach [46, 47] to B → D ∗ decays. We expect that still tobe a useful consideration here because the B c meson hassimilar properties to that of a heavy-light meson, interpo-lating between heavy-light and heavyonium [10]. Heavy-onium mesons also display similar mass-suppressed dif-ferences in amplitudes between ground-state vector andpseudoscalar mesons, so that we expect a B c → J/ψ de-1 − q / GeV − . − . − . . . . . A × | ~ p | / M J / ψ a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 . FIG. 4. Plot showing A ( q ) | (cid:126)p (cid:48) | /M J/ψ against q . The com-bination of A and | (cid:126)p (cid:48) | entering the helicity amplitudes inEq. (5). This combination does not exhibit as much noise as A alone, seen in Figure 3. cay to behave rather similarly to B → D ∗ .Physical particles with bc quark content, masses be-tween t / and t / − (i.e. above the physical q range forthe decay) and the appropriate quantum numbers to cou-ple to the current operator result in the appearance of asub-threshold pole in the corresponding form factor. Fol-lowing [46, 47] we include these poles in our fit using theform (Blaschke factor) P ( q ) = (cid:89) M pole z ( q , M ) . (26) P ( q ) is zero when q = M and has modulus 1 alongthe cut in the q plane for real q > t + where multi-particle production can occur. For the pole masses cor-responding to B c states we use the values listed in [4],which we reproduce here in table XV. We use 4 massesin each of the vector and axial channels and 3 in the pseu-doscalar channel. These values are derived from experi-ment and from lattice QCD calculations for the low-lying B c masses and from models for the higher-lying states.Since these masses are again simply numbers used in set-ting up our fit in z -space the values are used withoutuncertainties. As we will show in Section V B our resultsare not sensitive to the exact form of the pole expressionabove (Eq. (26)). The pole masses are shifted from thoseof the B c to account for our use of unphysical heavy-quark masses. We use M pole = M H c + M physpole − M phys B c again ensuring that when we extrapolate to the physi-cal heavy quark mass the physical continuum values of Table XV are recovered.We fit each form factor, F ( q ), to the fit function F ( q ) = 1 P ( q ) (cid:88) n =0 a n z n N n (27)where P ( q ) is the appropriate pole form for that formfactor (constructed using 1 − states for V ( q ), 1 + statesfor A ( q ) and A ( q ) or 0 − states for A ( q )) as inEq. (26). The remainder of the fit function is a poly-nomial in z with separate coefficients, a n , for each formfactor that take the form a n = (cid:88) j,k,l =0 b jkln ∆ ( j ) h (cid:18) am val c π (cid:19) k (cid:18) am val h π (cid:19) l . (28)The ∆ ( j ) h allow for the dependence on the heavy quarkmass using the η h mass as a physical proxy for this. Wehave ∆ (0) h = 1 and∆ ( j (cid:54) =0) h = (cid:18) M η h (cid:19) j − (cid:32) M phys η b (cid:33) j . (29)We take Λ = 0.5 GeV.The polynomials in am c /π and am h /π in Eq. (28) al-low for discretisation effects that are either constant withthe heavy quark mass (for example coming from the c -quark action) or that come from the heavy-quark actionand so depend on am h . Because of the form of the HISQaction only even powers of a can appear here.The remainder of Eq. (27), N n , takes into account theeffect of mistuning the valence and sea quark masses foreach form factor. N n = 1 + A n δ val m c + B n δ sea m c + C n δ sea m s + D n δ sea m l (30)with δ val m c = ( am val c − am tuned c ) /am tuned c ,δ sea m c = ( am sea c − am tuned c ) /am tuned c ,δ sea m s ( l ) = ( am sea s ( l ) − am tuned s ( l ) ) / (10 am tuned s ) . (31)We use the J/ψ meson mass to tune the c quark massfor the reasons discussed in [8]. We take am tuned c = am val c (cid:32) M phys J/ψ M J/ψ (cid:33) . (32)We use the η s mass [38] to determine the mistuning ofthe s quark mass in the sea, taking am tuned s = am val s (cid:32) M phys η s M η s (cid:33) . (33)We use the η s masses in lattice units from [11] to do this.The values of the physical J/ψ and η s masses that we2 .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . z . . . . . . A × P ( q ) a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 .
800 0 .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . z . . . . . A × P ( q ) a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 . .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . z − . − . . . . . . A × P ( q ) a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 .
800 0 .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . z . . . . . . V × P ( q ) a \ amh . . . . .
427 0 .
500 0 .
525 0 .
600 0 .
650 0 . FIG. 5. The points show our lattice QCD results for each form factor as given in Tables VIII, IX, X and XI multiplied by thepole function of Eq. (26) and plotted in z -space. The legend gives the mapping between symbol colour and shape and the setof gluon field configurations used, as given by the lattice spacing, and the heavy quark in lattice units. The blue curve witherror band is the result of our polynomial fit in z with lattice spacing and heavy quark mass dependence (Eq. (27)), evaluatingthe result in the continuum limit and for the b quark mass, to give the physical form factor for B c → J/ψ . use are given, with their uncertainties, in Table XIV. Todetermine the mistuning of the u/d = l quark mass inthe sea we take am tuned l = am tuned s / ∆[ m s /m l ] , (34)with ∆[ m s /m l ] = 27 . b n for each form factor,multiplying terms of order O ( a ) by 0 . a errors are removed in the HISQ action at tree-level [7]. We alsouse priors of 0 . .
5) for each B n , C n and D n for eachform-factor since sensitivity to sea quark masses enteronly at 1-loop. All remaining priors are taken as 0(1). Wehave checked that the prior width is conservative usingthe empirical Bayes criterion [50].In doing our fit to Eq. (27) we impose the kinematical3 TABLE XVI. Physical z -expansion coefficients for the pseu-doscalar, axial-vector and vector form factors for B c → J/ψ decay. The correlation matrices between the parameters aregiven in Appendix A. a a a a A A A V constraint2 M J/ψ A (0) = ( M J/ψ + M H c ) A (0) − ( M J/ψ − M H c ) A (0) . (35)This condition holds in the continuum limit with phys-ical sea quark masses for each heavy quark mass. Weimpose it using our lattice meson masses at each valueof am h and allowing for discretisation and quark massmistuning effects effects in this implementation. We dothis by imposing A (0) − ( aM J/ψ + aM H c ) / (2 aM J/ψ ) A (0) + (36)( aM J/ψ − aM H c ) / (2 aM J/ψ ) A (0) = ∆ kin . ∆ kin here is a nuisance term made up of leading orderdiscretisation and mistuning effects to account for theuse of lattice masses rather than values in the physicalcontinuum limit. We take∆ kin = (cid:88) i =1 α c,i ( am val c /π ) i + α h,i ( am h /π ) i + β c δ val m c + β (cid:48) c δ sea m c + β l δ sea m l + β s δ sea m s (37)where α and β are priors taken as 0(1). We find thatthe fit returns values for α and β well within their priorwidths.The fit was done simultaneously across all form fac-tors in order to preserve correlations important for con-structing helicity amplitudes. We use the nonlinear leastsquares fitting routine from the python package lsq-fit [51]. It returns a χ / dof of 0.18. The fit curvesevaluated at the b quark mass in the continuum limitare plotted along with the raw lattice results in Figure 3.In Figure 5 we show the lattice results and the fitfunction in z -space in a way that makes clearer how thefit is operating. The figure shows that the form factormultiplied by the pole function takes a very simple, lin-ear, form in z -space. The constant term in the z -spacepolynomial varies with heavy quark mass but the slopechanges very little.The extrapolated continuum values of a n = b n aregiven in Table XVI. In the continuum limit at the phys-ical b mass Eq. (27) becomes, for each form factor, F ( q ) = 1 P ( q ) (cid:88) n =0 a n z n . (38) N . . Γ / | η E W V c b | s − FIG. 6. Plot showing the stability of the total rate for B − c → J/ψµ − ν µ under variations of the correlator fits. The x axisvalue corresponds to N = δ +4 δ +16 δ +64 δ where δ n is thevalue of δ corresponding to the fit given in Table VII for set n .The black horizontal line and red error band correspond to ourfinal result and the blue points and blue error band correspondto the combination of fit variations associated to N . Ourresult for the total rate is very stable to these variations. P ( q ) is the appropriate pole factor for that case(Eq. (26)) using pole masses from Table XV. To recon-struct the form factors in q space both the values of the z -expansion coefficients from Table XVI and their corre-lations are needed; the correlation matrices are given inAppendix A. Table XVI shows that our fit is able to pindown the constant and linear pieces of the z -expansionbut is not able to pin down higher order terms in z . Thisis consistent with what we see in Figure 5. Notice thatthe values obtained for a and a in each case are wellwithin the prior width of 1.0 allowed for them in the fit. B. Tests of the Stability of the Analysis
Our results for the form factors are dependent to someextent upon choices made in the correlator fits as wellas choices made in the continuum/heavy-quark mass fit.Here we test the impact of those choices on both the con-tinuum fit coefficients and the final outcome of our calcu-lation to make sure that it is robust. It is convenenient inthat test to have a single quantity which we can compareand we take that to be the total rate of B − c → J/ψ(cid:96) − ν (cid:96) decay, i.e. Γ( B − c → J/ψµ − ν µ ) / | η EW V cb | . This is ob-tained by determining the helicity amplitudes from ourform factors and then integrating in q over the differen-tial rate they give (see Eqs (5) and (1)). The results forthe differential rates and total rate will be discussed inmore detail in Sec. VI; here we focus on the stability ofthe final result under variations of fit choice.We first look at the choices of the correlator fit param-eters: ∆ T , ∆ T J/ψ , ∆ T H c , the value of SVD cut andthe number of exponentials used in the fit. In order toverify that our results are independent of such choiceswe repeat the full analysis using all combinations of thevariations listed in Table VII. The total rate computed4 N . . . a A N − . − . − . a A N . . a V FIG. 7. As for Figure 6 showing the stability of the coeffi-cients of the z -expansion for the form factors under variationsof the correlator fits. We include a subset of coefficients here;other plots look very similar. using each of these fit variations is plotted in Figure 6,where we see that our final result is insensitive to suchvariations.Figure 7 shows a similar stability plot for a subsetof the z -expansion coefficients for separate form factors(Eq. 38). We again see that the z -expansion coefficientsare stable, within their uncertainty band, under changesin the correlator fits. Plots for other form factors notshown here are very similar.We now turn to the effect of choices in the extrapo-lation to the physical point. The pole form of Eq. (26)removes q dependence from poles in the q plane abovethe physical region for the decay. For b → c decays thesepoles are substantially above the physical region. For ex-ample, here q is (3 .
18 GeV) whereas the lowest pole, N poles012345 | a | | a A || a A || a A || a V | Γ 1 . . . . . . . . Γ / | η E W V c b | s − FIG. 8. Magnitude of the O ( z ) coefficient, a , for eachform factor plotted against the number of poles included inEq. (26). The prior widths on the b ijkn are scaled accordingto the number of poles, see text. Note that the maximumnumber of poles included for A is 3. The black crosses anderror bars give the total width for the (cid:96) = µ case, Γ / | η EW V cb | ,determined from that fit, using the right-hand y axis. Thegrey band corresponds to our final result for the total widthusing N poles = 4, and prior values for b ijkn of 0(1). This showshow the different coefficients as a function of N poles give avery stable result for the total width. for the pseudoscalar form factor is at the B c mass of6.275 GeV. Hence we do not expect the exact positionsor number of poles to have a big effect on the fits. Themagnitude of P ( q ) does depend on how many poles areincluded [4], however, and this affects the normalisationof the quantity form factor times P ( q ) that is expandedin z (Eq. (27)). This in turn affects the prior width thatmust be allowed on the z -expansion coefficients to achievea good fit.We have investigated the effect of including fewer polesin Eq. (26) by repeating our analysis including only thefirst N poles resonances listed in Table XV. We then take aprior width on the z -expansion coefficients of 5 . − N poles .We are able to obtain a good fit, with χ / dof ≈ . A expected below t + , we include only 3 poles for that form factor even inthe N poles case.Figure 8 shows these results, plotting against the left-hand y -axis the magnitude of the coefficient correspond-ing to the order z term, a , coming from the fits as afunction of the number of poles included. Results aregiven for each form factor. We see that as we includefewer poles, increasingly large z -expansion coefficients areneeded partly in order to account for the removal of phys-ical q dependence from missing poles but also becauseof the normalisation change.5 . . . . . . . . / | η EW V cb | s − O (2 , , , O (3 , , , O (3 , , , O (3 , , , O (2 , , , FIG. 9. Plot showing the stability of the total ratefor B − c → J/ψµ − ν µ considering lower order truncationsof z -expansion, discretisation and heavy mass dependentterms in Eq. (27) and Eq. (28). O ( n , n , n , n ) cor-responds to the result including terms of highest order O ((2Λ /M η h ) n , ( am c ) n , ( am h ) n , z n ). The vertical blackline is our final result, corresponding to O (3 , , , Figure 8 also gives, marked as black crosses and usingthe right-hand y axis, the total width for B c → J/ψ obtained from that fit. We see that the total width isvery stable as a function of N poles as the change in P ( q )is compensated by the different coefficients obtained forthe z -expansion.We have also investigated the effects of including fewer z terms in Eq. (27) as well as fewer am c , am h and2Λ QCD /M η h terms in Eq. (28). Figure 9 gives the to-tal width obtained using these variations, where we seethat our result is insensitive to the removal of the highestorder terms. VI. DISCUSSION
Using the form factors from Section V A we constructthe helicity amplitudes using Eq. (5). The helicity ampli-tudes are plotted in Figure 10, where we see that H and H t are singular at q = 0 as we would expect from the1 / (cid:112) q factors appearing in their definitions. This singu-lar behaviour is cancelled by the ( q − m (cid:96) ) factor ap-pearing in the differential decay rate (Eq. (1)). Togetherwith the factor of | (cid:126)p (cid:48) | this results in physical rates whichvanish at the extremes of the physical range, q = m (cid:96) and q = q Using these helicity amplitudes we construct the differ-ential decay rates using Eqs. (1), (2), (3) and (4). We take m µ = 105 . m e = 0 . d Γ /dq is plotted in Figure 12. In each case we plotthe rate for l = µ normalising each differential rate bythe total integrated decay rate. This means that any con-stant factors, such as | V cb | or B ( J/ψ → µ + µ − ) cancelout. Where an integration over q is necessary we usea simple trapezoidal interpolation in order to ensure co-variaces are carried through correctly, taking sufficientlymany points that the results are insensitive to using ad- q / GeV H / G e V H t H − H + H FIG. 10. Helicity amplitudes for B − c → J/ψ(cid:96) − ν (cid:96) plotted asa function of q . ditional points.The angular differential rates show the expected qual-itative picture for a pseudoscalar to vector meson de-cay [25]. The dependence on θ J/ψ is that for the casewhere the vector meson is seen in final states that carryspin (i.e. here in
J/ψ → µ + µ − ). This differs from, forexample, the case of B → D ∗ if the D ∗ is seen in Dπ .The θ W dependence is that of a b quark decay via the V − A weak interaction. This produces a virtual W − and charged lepton with, preferentially, a (1 + cos θ W ) distribution in the W − rest frame. The picture switchesfor a c decay producing a virtual W + [25].We also compute the forward-backward asymmetry, A F B , the lepton polarisation asymmetry (for the lepton (cid:96) produced from the virtual W ), A λ (cid:96) , and the fractionof J/ψ produced with longitudinal polarisation, F J/ψL .These are defined as A F B ( q ) = 1 d Γ /dq π (cid:90) π d Γ dq d cos( θ W ) cos( π − θ W ) dθ W , A λ (cid:96) ( q ) = d Γ λ (cid:96) = − / /dq − d Γ λ (cid:96) =+1 / /dq d Γ /dq ,F J/ψL ( q ) = d Γ λ J/ψ =0 /dq d Γ /dq (39)where we have chosen the forward direction for the pur-pose of A F B as being in the direction of the
J/ψ mo-mentum in the B c rest frame. These quantities are plot-ted in Figure 13 for the cases (cid:96) = e and (cid:96) = µ . Wesee in the top plot that A F B is negative except near q = 0 for (cid:96) = µ . This may be be understood fromFigure 10 where we see that H − H − , which is the dom-6 − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . θ J/ Ψ )0 . . . . . . d Γ ( B − c → J / ψ ( → µ + µ − ) µ − ν ) d c o s ( θ J / Ψ ) / Γ − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . θ W )0 . . . . . . . . . d Γ ( B − c → J / ψ ( → µ + µ − ) µ − ν ) d c o s ( θ W ) / Γ − − − χ . . . . . . . d Γ ( B − c → J / ψ ( → µ + µ − ) µ − ν ) d χ / Γ FIG. 11. Angular differential decay rates for B − c → J/ψ(cid:96) − ν (cid:96) for the (cid:96) = µ case. From the top down, d Γ /dq d cos( θ J/ψ ), d Γ /dq d cos( θ W ) and d Γ /dq dχ . Each rate is normalised bythe total decay rate Γ( B − c → J/ψ ( → µ + µ − ) µ − ν µ ). inant contribution to A F B for m (cid:96) << q , is less than orequal to zero across the full physical q range. The be-haviour of A F B near q = 0 in the (cid:96) = µ case originatesfrom the − m (cid:96) /q H t H cos( θ W ) term in Eq. 2. When m (cid:96) /q ≈ O (1) it is apparent from Figure 10 that thisterm will dominate over the H − H − contribution. Inthe middle plot of Figure 13 we see that A λ e is equal tounity across the full q range, in line with the expecta-tion that in the massless limit the lepton is produced ina purely left handed helicity eigenstate. In the bottom q / GeV . . . . . . . . . Γ d Γ d q / G e V − d Γ dq ( B − c → J/ψµ − ¯ ν µ ) / Γ( B − c → J/ψµ − ¯ ν µ ) FIG. 12. The differential rate d Γ /dq for B − c → J/ψ(cid:96) − ν (cid:96) for the (cid:96) = µ case, normalised by the total decay rate Γ. plot of Figure 13 we see that the longitudinal polarisa-tion fraction, F J/ψL , approaches unity near q = 0 where H and H t dominate the total rate, and goes to 1 / q where H = H + = H − and H t = 0.We also compute the total decay rate for the cases (cid:96) = e and (cid:96) = µ . We findΓ( B − c → J/ψµ − ν µ ) | η EW V cb | = 1 . × s − (40)= 11 . × − GeV . and Γ( B − c → J/ψe − ν e ) | η EW V cb | = 1 . × s − (41)= 11 . × − GeVwith the ratio Γ e / Γ µ = 1 . m µ amount to 0 .
4% of the rate. The errorbudget for the total rate for the (cid:96) = µ case is given inTable XVII. Not surprisingly we see that the largest con-tributions come from the discretisation effects from theheavy quark mass, the statistical uncertainties on thefinest (ultrafine) lattices and the quark mass mistuningterm, predominantly that from the sea quarks. Both ofthese uncertainties can be reduced in future by increas-ing statistics on ultrafine lattices, adding more results onlattices that include physical u/d quarks and obtainingresults on exafine ( a =0.03 fm) lattices.We can compare our results for the total rate to thosefrom earlier, non-lattice QCD approaches such as po-tential models and QCD sum rules. In these other ap-proaches it is much harder to quantify sources of uncer-tainty and derive an error budget. One way to obtain a7 q / GeV − . − . − . . . A F B SM ‘ = e SM ‘ = µ q / GeV . . . . . . A λ ‘ SM ‘ = e SM ‘ = µ q / GeV . . . . . . . . F J / ψ L SM ‘ = e SM ‘ = µ FIG. 13. Angular asymmetry variables defined in Eq. (39)for the cases (cid:96) = e and (cid:96) = µ . TABLE XVII. Error budget for the total rate for B c → J/ψµν µ . Errors are given as a percentage of the final answer.The top half gives sources of systematic errors: Extrapolationof data in M η h to the physical value of M η b , errors propor-tional to powers of am c , am h , mass mistuning effects, anduncertainties in determination of the lattice spacing. The sec-ond section gives a breakdown of the statistical uncertaintiescorresponding to each of our data sets. Finally ‘Other priors’includes all of the remaining sources of uncertainty such as∆ kin and the current renormalisation factors Z .Source % errorM η h → M η b am c → am h → δ m w /a , w picture of possible systematic errors is simply to compareresults from different model calculations. Ref. [53] usesa relativistic quark model to obtain a value 10 . × − GeV for Γ / | V cb | and gives a table of comparison to othermodel approaches. Results vary over a range up to dou-ble this value, which implies a worrying lack of control ofsystematic effects in at least some of these calculations.In future it may be possible for such calculations to betuned against our lattice QCD result for B c → J/ψ andthen yield useful information on B c decay to excited char-monium which is much harder to determine accurately inlattice QCD.From our value for Γ / | η EW V cb | in Eq. (40) we canderive a result for the total width for B − c → J/ψµ − ν µ using a value for η EW (from Section II) and for | V cb | . Wetake | V cb | = 41 . . × − from an average of inclusiveand exclusive determinations and with the error scaledby 2.4 to allow for their inconsistency [52]. This givesΓ( B − c → J/ψµ − ν µ ) = 2 . × s − . (42)Here the first uncertainty is from our lattice QCD calcu-lation and the second from the uncertainty in | η EW V cb | .These uncertainties are approximately equal.Using the experimental average value for the B c life-time of 0 . × − s [52, 54, 55] then gives thebranching fractionBr( B − c → J/ψµ − ν µ ) = 0 . . (43)The uncertainties are from lattice QCD, | η EW V cb | andthe lifetime respectively.8The accuracy of our result means that it can be usedto normalise other B c decay modes. In the Particle DataTables [52] B c branching fractions are given in the mod-ified form B ( b → B + c )Γ i / Γ. Here the first factor is theprobability for a produced b quark to hadronise as a B c or one of its excitations. Using the value of the modi-fied branching fraction for B − c → J/ψµ − ν µ obtained byCDF [56] of 8 . . × − [52] we infer B ( b → B + c ) = 0 . . (44)The first uncertainty is the combined uncertainty fromEq. (43) and the second uncertainty is from the exper-imental value of the modified branching fraction. Thisvalue can be compared to the value B ( b → B + ) =0 . b quark will form a B c .We have focussed here on the rates for production of a µ or e from the virtual W in the final state. We will dis-cuss results for B − c → J/ψτ − ν τ in a follow-up paper [22]. VII. CONCLUSIONS
We have carried out the first lattice QCD calcula-tion of the form factors within the Standard Model for B − c → J/ψ(cid:96) − ν (cid:96) decay. This is a process under activestudy by LHCb [1, 57] and more accurate form factorsare needed (and crucially with quantified errors) to re-duce uncertainties in the experimental determination ofthe lepton-universality violating ratio R ( J/ψ ).We give all four form factors across the full q rangeof the decay. The calculation is done using the HISQaction for all quarks, enabling us to normalise the lat-tice current operators that couple to the W fully non-perturbatively. We have used the heavy-HISQ approach,obtaining results at multiple values of the heavy quarkmass on a range of fine lattice spacings, fitting the heavy-quark mass dependence to obtain a physical result for B c → J/ψ in the continuum limit.We give each form factor in a parameterisation thatallows it to be fully reconstructed, including correlationsbetween the form factors (see Appendix A). This yieldsa total rate integrated over q for a light lepton in thefinal state, up to a CKM factor, of (repeating Eq. (40))1 | η EW V cb | Γ( B − c → J/ψµ − ν µ ) = 1 . × s − . (45)The 6% uncertainty is broken down in an error budget inTable XVII. Our result for the total rate gives a branch-ing fraction (repeating Eq. (43))Br( B − c → J/ψµ − ν µ ) = 0 . . (46)The first two uncertainties, from our lattice QCD cal-culation and from the current uncertainty in V cb are equalhere.We will discuss the implications of our results for thecase of a heavy ( τ ) lepton in the final-state elsewhere [22]. TABLE XVIII. Correlation matrix for z -expansion coeffi-cients of A a A a A a A a A a A a A -0.2828 1.0 -0.4169 -0.01778 a A -0.01528 -0.4169 1.0 -0.003817 a A z -expansion coefficientsof A A a A a A a A a A a A a A -0.06628 0.2101 -0.234 0.01736 a A a A -0.001229 -0.02275 0.0433 0.01512 We have demonstrated that the heavy-HISQ method-ology is a viable procedure for these calculations andthat the statistical challenge of analysing many latticeQCD correlation functions, while preserving covariancesimportant for constructing and fitting q and heavy massdependence, can be handled. This work forms an im-portant precursor for the calculation of B s → D ∗ s formfactors across the full q range using the same machin-ery, which is underway. This will allow the dominantuncertainty from external inputs to be reduced in thedetermination of V cb from experimental results [58]. Acknowledgements
We are grateful to the MILC collaboration for the useof their configurations and code. We thank C. Bouchard,B. Colquhoun, J. Koponen, P. Lepage, E. McLean andC. McNeile for useful discussions. Computing was doneon the Cambridge service for Data Driven Discovery(CSD3), part of which is operated by the University ofCambridge Research Computing on behalf of the DIRACHPC Facility of the Science and Technology FacilitiesCouncil (STFC). The DIRAC component of CSD3 wasfunded by BEIS capital funding via STFC capital grantsST/P002307/1 and ST/R002452/1 and STFC operationsgrant ST/R00689X/1. DiRAC is part of the nationale-infrastructure. We are grateful to the CSD3 supportstaff for assistance. Funding for this work came fromthe UK Science and Technology Facilities Council grantsST/L000466/1 and ST/P000746/1.
Appendix A: Reconstructing the Fit
Our parameterisation of the form factors for B c → J/ψ in the continuum limit is given in Eq. (38). It consists9
TABLE XX. Correlation matrix for z -expansion coefficientsof A A a A a A a A a A a A -0.2576 0.1908 0.02885 0.01033 a A -0.03401 -0.2558 0.03141 -0.01381 a A a A -0.01442 0.02789 -0.07449 -0.007872TABLE XXI. Correlation matrix for z -expansion coefficientsof A V .Corr a V a V a V a V a A a A -0.005599 0.01376 0.0007743 4.261e-05 a A -0.0007073 0.003171 0.00062 3.389e-05 a A z -expansion coefficientsof A
1. Corr a A a A a A a A a A a A -0.09714 1.0 -0.4303 0.04665 a A -0.0211 -0.4303 1.0 -0.1168 a A z -expansion coeffi-cients of A A a A a A a A a A a A a A a A -0.3849 0.5746 0.5712 0.04731 a A z -expansion coeffi-cients of A V .Corr a V a V a V a V a A a A -0.0007464 0.0106 0.0002107 7.768e-06 a A -0.001222 0.004827 0.0002132 1.101e-05 a A -0.002162 0.0002878 -3.916e-05 -2.191e-06 TABLE XXV. Correlation matrix for z -expansion coefficientsof A
2. Corr a A a A a A a A a A a A -0.5824 1.0 0.04233 0.04547 a A -0.1636 0.04233 1.0 -0.1036 a A -0.03372 0.04547 -0.1036 1.0TABLE XXVI. Correlation matrix for z -expansion coeffi-cients of A V .Corr a V a V a V a V a A a A a A a A z -expansion coeffi-cients of V .Corr a V a V a V a V a V a V -0.2571 1.0 -0.35 -0.01055 a V a V of a pole factor with no uncertainty and a polynomialin z for which the coefficients with their uncertaintiesare given in Table XVI. In this section we give the cor-relations between the z -expansion coefficients which arenecessary for reconstructing our results explicitly, as wellas instructions for using the included ancillary files toload the z -expansion parameters together with their cor-relations automatically into python [59].The correlation between two coefficients is defined inthe usual way asCorr( X, Y ) = (cid:104) ( ¯ X − X )( ¯ Y − Y ) (cid:105) (cid:112) σ ( X ) σ ( Y ) (A1)where σ ( X ) is the variance of X and ¯ X is the mean of X . The values are tabulated in Tables XVIII to XXVII.In this calculation and in the ancillary files we usethe gvar python package to track and propagate corre-lations. Included in the ancillary files are two text files; CORRELATIONS.txt contains a dictionary includ-ing the means and variances of the z -expansion param-eters on the first line and a dictionary detailing the cor-relations between these parameters on the second line, CHECKS.txt contains arrays of q values and formfactor mean and standard deviation values at the cor-responding values of q . This file is used by the pythonscript load fit.py as a simple check that the fit has beenloaded correctly. Running python load fit.py will loadthe parameters from CORRELATIONS.txt and com-pare values computed at hard coded intervals in q tothose in CHECKS.txt which were computed as part ofthis work. Running python load fit.py will also pro-duce some simple plots of the form factors across thefull q range. We have tested load fit.py using Python3.7.5 [59], gvar 9.2.1 [60], numpy 1.18.2 [61] and mat-plotlib 3.1.2 [62].0 [1] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 121801(2018), arXiv:1711.05623 [hep-ex].[2] C. Bernard et al. , Phys. Rev. D , 014506 (2009),arXiv:0808.2519 [hep-lat].[3] J. A. Bailey et al. (Fermilab Lattice, MILC), Phys. Rev.D , 114504 (2014), arXiv:1403.0635 [hep-lat].[4] J. Harrison, C. Davies, and M. Wingate (HPQCD),Phys. Rev. D97 , 054502 (2018), arXiv:1711.11013 [hep-lat].[5] D. Bigi, P. Gambino, and S. Schacht, Phys. Lett.
B769 ,441 (2017), arXiv:1703.06124 [hep-ph].[6] M. Bordone, M. Jung, and D. van Dyk, Eur. Phys. J. C , 74 (2020), arXiv:1908.09398 [hep-ph].[7] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P.Lepage, J. Shigemitsu, H. Trottier, and K. Wong(HPQCD, UKQCD), Phys. Rev. D75 , 054502 (2007),arXiv:hep-lat/0610092 [hep-lat].[8] D. Hatton, C. Davies, B. Galloway, J. Koponen, G. Lep-age, and A. Lytle (HPQCD), (2020), arXiv:2005.01845[hep-lat].[9] C. McNeile, C. T. H. Davies, E. Follana, K. Hornbostel,and G. P. Lepage (HPQCD), Phys. Rev.
D85 , 031503(2012), arXiv:1110.4510 [hep-lat].[10] C. McNeile, C. Davies, E. Follana, K. Hornbostel,and G. Lepage, Phys. Rev. D , 074503 (2012),arXiv:1207.0994 [hep-lat].[11] E. McLean, C. T. H. Davies, A. T. Lytle, and J. Kopo-nen, Phys. Rev. D99 , 114512 (2019), arXiv:1904.02046[hep-lat].[12] E. McLean, C. Davies, J. Koponen, and A. Lytle, Phys.Rev. D , 074513 (2020), arXiv:1906.00701 [hep-lat].[13] D. Hatton, C. Davies, G. Lepage, and A. Lytle, in (2019)arXiv:1908.10116 [hep-lat].[14] D. Hatton, C. Davies, G. Lepage, and A. Ly-tle (HPQCD), Phys. Rev. D , 114513 (2019),arXiv:1909.00756 [hep-lat].[15] L. J. Cooper, C. T. H. Davies, J. Harrison, J. Komijani,and M. Wingate, (2020), arXiv:2003.00914 [hep-lat].[16] J. P. Lees et al. (BaBar), Phys. Rev. Lett. , 101802(2012), arXiv:1205.5442 [hep-ex].[17] J. T. Wei et al. (Belle), Phys. Rev. Lett. , 171801(2009), arXiv:0904.0770 [hep-ex].[18] J. P. Lees et al. (BaBar), Phys. Rev.
D86 , 032012 (2012),arXiv:1204.3933 [hep-ex].[19] J. P. Lees et al. (BaBar), Phys. Rev.
D88 , 072012 (2013),arXiv:1303.0571 [hep-ex].[20] P. Gambino, M. Jung, and S. Schacht, Phys. Lett.
B795 ,386 (2019), arXiv:1905.08209 [hep-ph].[21] G. Caria et al. (Belle), Phys. Rev. Lett. , 161803(2020), arXiv:1910.05864 [hep-ex].[22] J. Harrison, C. Davies, and A. Lytle, in preparation(2020).[23] T. D. Cohen, H. Lamm, and R. F. Lebed, Phys. Rev.
D98 , 034022 (2018), arXiv:1807.00256 [hep-ph].[24] A. Sirlin, Nucl. Phys. B , 83 (1982).[25] J. D. Richman and P. R. Burchat, Rev. Mod. Phys. ,893 (1995), arXiv:hep-ph/9508250 [hep-ph].[26] G. C. Donald, C. T. H. Davies, J. Koponen, andG. P. Lepage (HPQCD), Phys. Rev. D90 , 074506 (2014),arXiv:1311.6669 [hep-lat]. [27] A. Bazavov et al. (MILC), Phys. Rev.
D87 , 054505(2013), arXiv:1212.4768 [hep-lat].[28] A. Bazavov et al. (MILC), Phys. Rev.
D82 , 074501(2010), arXiv:1004.0342 [hep-lat].[29] R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. Mc-Neile, Phys. Rev.
D88 , 074504 (2013), arXiv:1303.1670[hep-lat].[30] B. Chakraborty, C. T. H. Davies, P. G. de Oliviera, J. Ko-ponen, G. P. Lepage, and R. S. Van de Water, Phys. Rev.
D96 , 034516 (2017), arXiv:1601.03071 [hep-lat].[31] B. Chakraborty, C. T. H. Davies, B. Galloway, P. Knecht,J. Koponen, G. C. Donald, R. J. Dowdall, G. P. Lep-age, and C. McNeile, Phys. Rev.
D91 , 054508 (2015),arXiv:1408.4169 [hep-lat].[32] C. T. Sachrajda and G. Villadoro, Phys. Lett.
B609 , 73(2005), arXiv:hep-lat/0411033 [hep-lat].[33] D. Guadagnoli, F. Mescia, and S. Simula, Phys. Rev. D , 114504 (2006), arXiv:hep-lat/0512020.[34] C. Monahan, J. Shigemitsu, and R. Horgan, Phys. Rev. D87 , 034017 (2013), arXiv:1211.6966 [hep-lat].[35] G. P. Lepage, corrfitter Version 8.0.2(github.com/gplepage/corrfitter).[36] R. Dowdall, C. Davies, R. Horgan, G. Lepage, C. Mona-han, J. Shigemitsu, and M. Wingate, Phys. Rev. D ,094508 (2019), arXiv:1907.01025 [hep-lat].[37] M. Tanabashi et al. (Particle Data Group), Phys. Rev.
D98 , 030001 (2018).[38] R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. Mc-Neile, Phys. Rev.
D88 , 074504 (2013), arXiv:1303.1670[hep-lat].[39] R. Aaij et al. (LHCb), Phys. Rev. D , 032005 (2017),arXiv:1612.07421 [hep-ex].[40] A. M. Sirunyan et al. (CMS), Phys. Rev. Lett. ,132001 (2019), arXiv:1902.00571 [hep-ex].[41] R. Aaij et al. (LHCb), Phys. Rev. Lett. , 232001(2019), arXiv:1904.00081 [hep-ex].[42] R. J. Dowdall, C. T. H. Davies, T. C. Hammant,and R. R. Horgan, Phys. Rev. D86 , 094510 (2012),arXiv:1207.5149 [hep-lat].[43] E. J. Eichten and C. Quigg, Phys. Rev. D , 5845(1994), arXiv:hep-ph/9402210.[44] S. Godfrey, Phys. Rev. D , 054017 (2004), arXiv:hep-ph/0406228.[45] N. Devlani, V. Kher, and A. Rai, Eur. Phys. J. A ,154 (2014).[46] C. Boyd, B. Grinstein, and R. F. Lebed, Phys. Rev. D , 6895 (1997), arXiv:hep-ph/9705252.[47] I. Caprini, L. Lellouch, and M. Neubert, Nucl. Phys. B , 153 (1998), arXiv:hep-ph/9712417.[48] R. J. Hill, Proceedings, 4th Conference on Flavor Physicsand CP Violation (FPCP 2006): Vancouver, BritishColumbia, Canada, April 9-12, 2006 , eConf
C060409 ,027 (2006), arXiv:hep-ph/0606023 [hep-ph].[49] A. Bazavov et al. , Phys. Rev.
D98 , 074512 (2018),arXiv:1712.09262 [hep-lat].[50] G. Lepage, B. Clark, C. Davies, K. Hornbostel,P. Mackenzie, C. Morningstar, and H. Trottier, Nucl.Phys. B Proc. Suppl. , 12 (2002), arXiv:hep-lat/0110175.[51] G. P. Lepage, lsqfit Version 11.4(github.com/gplepage/lsqfit). [52] P. Zyla et al. (Particle Data Group), Prog. Theor. Exp.Phys. , 083C01 (2020).[53] D. Ebert, R. Faustov, and V. Galkin, Phys. Rev. D ,094020 (2003), arXiv:hep-ph/0306306.[54] R. Aaij et al. (LHCb), Phys. Lett. B , 29 (2015),arXiv:1411.6899 [hep-ex].[55] A. M. Sirunyan et al. (CMS), Eur. Phys. J. C ,457 (2018), [Erratum: Eur.Phys.J.C 78, 561 (2018)],arXiv:1710.08949 [hep-ex].[56] T. A. Aaltonen et al. (CDF), Phys. Rev. D , 052001(2016), arXiv:1601.03819 [hep-ex].[57] R. Aaij et al. (LHCb), (2018), arXiv:1808.08865 [hep-ex]. [58] R. Aaij et al. (LHCb), Phys. Rev. D , 072004 (2020),arXiv:2001.03225 [hep-ex].[59] G. Van Rossum and F. L. Drake, Python 3 ReferenceManual (CreateSpace, Scotts Valley, CA, 2009).[60] G. P. Lepage, gvar Version 9.2.1(github.com/gplepage/gvar).[61] S. van der Walt, S. C. Colbert, and G. Varoquaux, Com-puting in Science Engineering , 22 (2011).[62] J. D. Hunter, Computing in Science & Engineering9