$B_c \rightarrow J/\psi$ Form Factors for the full $q^2$ range from Lattice QCD

We present the first lattice QCD determination of the $B_c \rightarrow J/\psi$ vector and axial-vector form factors. These will enable experimental information on the rate for $B_c$ semileptonic decays to $J/\psi$ to be converted into a value for $V_{cb}$. Our calculation covers the full physical $q^2$ range of the decay and uses non-perturbatively renormalised lattice currents. We use the Highly Improved Staggered Quark (HISQ) action for all valence quarks on the second generation MILC ensembles of gluon field configurations including $u$, $d$, $s$ and $c$ HISQ sea quarks. Our HISQ heavy quarks have masses ranging upwards from that of $c$; we are able to reach that of the $b$ on our finest lattices. This enables us to map out the dependence on heavy quark mass and determine results in the continuum limit at the $b$. We use our form factors to construct the differential rates for $B_c^- \rightarrow J/\psi \mu^- \overline{\nu}_\mu$ and obtain a total rate with 6\% uncertainty: $\Gamma(B_c^-\rightarrow J/\psi \mu^-\overline{\nu}_{\mu})/|\eta_{\mathrm{EW}}V_{cb}|^2 = 1.74(12)\times 10^{13} ~\mathrm{s}^{-1}$. Including values for $V_{cb}$, $\eta_{\mathrm{EW}}$ and $\tau_{B_c}$ yields a branching fraction for this decay mode of 0.0151(10)(10)(3) ~with uncertainties from lattice QCD, $\eta_\mathrm{EW}V_{cb}$ and $\tau_{B_c}$ respectively.


I. INTRODUCTION
Accurate calculations of hadronic parameters are needed for the determination of Cabibbo-Kobayashi-Maskawa (CKM) matrix elements from the comparison of Standard Model theory and experimental results for exclusive flavor-changing decay rates. Leptonic decay rates require decay constants to be determined and semileptonic decay rates require form factors. Lattice QCD is now established as the method of choice for the calculation of these hadronic parameters, and efforts are ongoing both to improve the accuracy of the results and to expand the range of processes for which calculational results are available. Here we report on the first lattice QCD calculation of the form factors for B c to J=ψ semileptonic decay, a process under active study by the LHCb experiment [1]. Because the valence quarks involved in this process are all heavy, the calculation is under very good control in lattice QCD as we will show. This opens the prospect of a new exclusive process that can be used for the determination of jV cb j that has a reduced theoretical uncertainty.
Currently exclusive determinations of jV cb j are focused on B → D Ã and B → D semileptonic decays. The emphasis is on the former pseudoscalar to vector meson transition because of more favorable kinematic factors for the differential decay rate towards the zero recoil region. Lattice QCD calculations are generally more accurate in this region because the daughter meson has small spatial momentum in the lattice frame (where the parent meson is usually taken to be at rest). The emphasis 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][3][4]. Comparison was then made to results derived from experimental differential rates in this same limit. More recently (see, for example, [5,6]) it has become clear that extrapolating the experimental results to the zero recoil point comes with significant systematic uncertainties, associated with the underlying model dependence of such extrapolations, that were previously being underestimated. The way forward requires a much more direct comparison q 2 -bin by q 2 -bin of the experimental decay rate and that from theory, determined from lattice QCD form factors. For this we need lattice QCD form factors as a function of q 2 and preferably covering the full q 2 range of the decay. This allows the uncertainty in the extracted value of jV cb j to be optimized between the changing experimental and theoretical uncertainties as a function of q 2 . We show here that it is possible to calculate the form factors for B c → J=ψ semileptonic decay across the full q 2 range in lattice QCD.
The B c → J=ψ form factor calculation that we describe here acts as a prototype for upcoming calculations of form factors for B → D Ã and B s → D Ã s . B c → J=ψ is a more 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, from which appropriate correlation functions are constructed, are inexpensive to calculate. The correlation functions then have small statistical errors, even when the daughter has maximum spatial momentum. The absence of valence light quarks means that finite-volume effects are negligible, and sensitivity to the u=d quark mass in the sea should be small.
The main obstacle for the calculation of B c → J=ψ form factors is that of the discretization effects associated with the heavy quarks. The c quark is handled very accurately in lattice QCD as long as improved discretizations of the Dirac equation are used. A particularly accurate approach is that of HPQCD's highly improved staggered quark (HISQ) action [7], and it is the one that we will use here for all quarks. For c quarks a recent calculation using HISQ [8] obtained a 0.4% uncertainty in the J=ψ decay constant and showed good control of lattice discretization effects all the way to very coarse lattices with a spacing of 0.15 fm. Discretization effects are larger for the b quark. Indeed, since we expect discretization effects to grow as a power of the heavy quark mass in lattice units, am h , we need to work on fine lattices to approach the b quark mass, and this is what we will do here. The dominant discretization effects in the HISQ action behave as a 4 , since tree-level a 2 errors are removed, and those that include powers of α s are heavily suppressed. This means that discretization effects can be controlled for quarks with masses m h between that of the c and the b on lattices with lattice spacing below 0.1 fm. By working with a range of quark masses reaching up to that of the b on a range of lattice spacings, we can map out both the dependence on m h and the dependence on am h , both of which are smooth functions, and obtain a physical value in the continuum limit for the b quark. This "heavy-HISQ" approach was developed by HPQCD for B meson decay constants [9,10] and is now being extended to form factors [11]. Its efficacy has been demonstrated for B s → D s form factors for the full q 2 range in [12], and here we apply it to B c → J=ψ. A big advantage of this approach is that the lattice current operators that couple to the W boson can be normalized fully nonperturbatively, e.g., [12][13][14][15], avoiding the systematic errors associated with the perturbative normalization needed if a nonrelativistic approach is used for the b quark in lattice QCD [3,4].
A further motivation for studying B c → J=ψ semileptonic decay in detail is to calculate the ratio, RðJ=ψÞ, of the partial widths for the outgoing lepton to be a τ compared to that for it to be an e or μ. The analogous results for B decays, RðDÞ and RðD Ã Þ, e.g., [16][17][18][19], have been a source of tension between experiment and the Standard Model [20], implying lepton-universality violation. Recent results from Belle, on the contrary, show good agreement with the Standard Model [21]. This makes it very important to test lepton-universality violation in other processes, and we can obtain RðJ=ψÞ from our form factors for comparison with ongoing LHCb analyses. We will present those results and analyses of other lepton-universality violation tests separately [22]; here we focus on the form factors and differential rates for W decay to μ or e.
The subsequent sections are organized as follows: (i) In Sec. II we begin by outlining the relevant experimental observables and relate them to the invariant form factors coming from QCD matrix elements. (ii) Section III gives an overview of methods generic to the extraction of matrix elements from HISQ threepoint and two-point correlation functions in lattice QCD, and discusses our choices of operators, polarizations and momenta appropriate to extract the form factors specified in II. (iii) Section IV details the specifics of our lattice calculation including our nonperturbative current renormalization. (iv) In Sec. V we present the direct results of the lattice calculation and discuss the extraction of the physical continuum form factors as a function of q 2 . (v) Finally in Secs. VI and VII we use our form factors to compute the physical differential rates for B − c → J=ψμ −ν μ and discuss the significance of these results and implications for future work.

II. THEORETICAL BACKGROUND
Here we give the partial rates for where l is the final state lepton as differentials with respect to q 2 and angular variables defined in the standard way in Fig. 1. In this work we consider only the cases l ¼ μ and l ¼ e. The partial rates are obtained from the full differential decay rate assuming the J=ψ decay is purely electromagnetic and summing over μ þ μ − helicities (assuming that the μ þ μ − are massless and hence pure helicity eigenstates) [23]. This gives where N J=ψ and N B c are amplitudes proportional to the decay constant of the corresponding meson and ϵ is the J=ψ polarization vector. We will make use of these when we come to extract the form factors in Eq. (6) from our lattice correlation functions.

III. COMPUTATIONAL STRATEGY
We follow the strategy of previous heavy-HISQ calculations of form-factors [11,12,15]. We work with multiple heavy 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 are able to reach a mass very close to the physical b quark mass on our finest lattices. An important point is that, at every lattice spacing, we are able to cover the full q 2 range of the decay for the heavy quark masses that we can reach at that lattice spacing; i.e., the q 2 range expands on finer lattices in step with the heavy quark mass range [12]. We compute form factors by extracting combinations of the relevant matrix elements defined in Eq. (6) from correlation functions computed on the lattice. We then fit the form factors as a function of lattice spacing and heavy quark mass to determine their functional form in the continuum limit at the physical b mass.
By considering the insertion of complete sets of states we may express these correlation functions in terms of general fit forms (folding the two-point functions about the midpoint of the lattice so that t extends from 0 to N t =2), Here n, m correspond to on shell particle states with quantum numbers resulting in nonzero amplitudes and the o labels indicate energies and amplitudes corresponding to the time-doubled states typical of staggered quarks. The lowest energy, nonoscillating states are those corresponding to the J=ψ and H − c . We work with the H − c at rest, choosing to access the full q 2 range by giving momentum only to the J=ψ, and extract the matrix elements of these states from our lattice correlators. This gives and J 00 ðν;ΓÞ ¼ where ⃗ p 0 ν is the ν component of the J=ψ spatial momentum, with ν corresponding to the choice of polarization in Eq. (10), with currentcΓh. In the subsequent subsections we discuss the combinations of ν and Γ for which we must compute correlation functions in order to extract the full set of correlation functions defined in Eq. (6).

A. Extracting Vðq 2 Þ
The choice of operators used to extract Vðq 2 Þ is the most straightforward since the matrix element of the vector In this calculation we give the J=ψ spatial momentum ⃗ p 0 ¼ ðk; k; 0Þ. In order to isolate all the form factors we need one component of ⃗ p 0 to be zero. Keeping both of the others non-zero minmizes the discretization errors for a given magnitude of p 0 . Here we choose μ ¼ 3 and ν ¼ 1 and find where we have defined the relativistic normalization, with k the ν component of p 0 .

B. Extracting
In order to isolate A 0 ðq 2 Þ, following [26], we make use of the partially conserved axial-vector current (PCAC) Taking Γ μ ¼ γ 5 and ν ¼ 1 in Eq. (14) and multiplying by m c þ m b we then have C. Extracting A 1 ðq 2 Þ In order to isolate A 1 we use the axial-vector current Γ ¼cγ μ γ 5 h and J=ψ operatorcγ ν c and choose μ ¼ ν ¼ 3 along the spatial direction with zero J=ψ momentum. Using Eq. (6) this gives X which gives, in terms of J 00 D. Extracting A 2 ðq 2 Þ The extraction of A 2 is more complicated than the extraction of the other form factors since no trivial choice of directions in axial-vector and J=ψ operators isolates the contribution of A 2 relative to A 1 or A 0 . We use axial-vector current operator Γ ¼cγ 1 γ 5 h and J=ψ operatorcγ 1 c. This yields contributions from each form factor, ΦðkÞJ 00 We must then subtract the A 0 and A 1 contributions to obtain A 2 .

IV. LATTICE CALCULATION
In this section we give details of the gauge field configurations used in this calculation, as well as the values of masses, momenta and staggered spin-taste operators used in the calculation. We use the second generation MILC gluon ensembles including light, strange and charm sea quarks [27,28] and compute HISQ charm and heavy quark propagators. The details of these gauge configurations are given in Table I. As stated earlier, we work with the H c at rest on the lattice. The arrangement of operators in the three point correlation functions is shown in Fig. 2. We refer to c 1 as the "spectator" charm quark, c 2 as the "active" quark and h as the "extended" heavy quark. We use a single random wall source at time T for both the spectator and active quark, with a phase patterning the source for the active quark to achieve an appropriate staggered spin-taste for the meson. We use the spectator propagator at time 0 as the source for the extended heavy quark propagator with a spin-taste operator corresponding to the H c quantum numbers and then tie the heavy quark propagator together TABLE I. Details of the gauge field configurations used in our calculation [27,28]. Set 1 is referred to as "fine," set 2 as "superfine," set 3 as "ultrafine" and set 4 as "physical fine." The physical value of ω 0 was determined in [29] to be 0.1715(9) fm and the values of ω 0 =a were taken from [12,30,31]. with the active charm quark propagator at the current at time t.
Following the notation of [7] for staggered spin-taste operators, we choose lattice meson and current operators which have the same quantum numbers as those detailed in Sec. III. We always take the currents to be the local ones, γ i ⊗ γ i and γ i γ 5 ⊗ γ i γ 5 . This means that the current insertions only need to be renormalized once. In order to make the desired matrix element an overall taste singlet, this requires that we use different tastes of J=ψ and H c . The specific combinations we use are given in Table II. This means that we have three different tastes of J=ψ and two different tastes of H c . Different tastes of a meson differ in mass by discretization effects that are small for HISQ quarks [7,27] and particularly small for heavy mesons. When fitting two-and three-point correlators the energies corresponding to different spin-taste operators use separate priors (i.e., we assume that they are different) except for γ 3 ⊗ γ 3 and γ 1 ⊗ γ 1 whose energies we fix to be equal. We will show later how close the different taste meson masses are.
We also calculate pseudoscalar heavyonium two-point correlation functions with spin-taste γ 5 ⊗ γ 5 . This allows us to determine the mass of the ground-state heavyonium meson that we denote η h . The mass of the η h is a useful physical proxy for the heavy quark mass when we map out the dependence on heavy quark mass of our results.
We choose values of am h for the heavy quark masses at each value of the lattice spacing following [11]. They range from above the c mass up to a value of 0.8. This corresponds to a growing range in m h as the lattice spacing becomes finer. The valence masses used to compute propagators are given in Table III. We also choose the values of T for our 3-point correlation functions following [11].
Since we must compute heavy quark propagators at multiple masses, from multiple T, the choice to insert momentum via the active charm quark minimises the number of inversions. We put momentum on the c quark propagator via a momentum twist [32,33]. The twists were chosen to span evenly the physical q 2 range for the largest value of am h , approximated from aM H c values given in [11] and the physical J=ψ mass. The values of twists and T are given in Table IV. Note that because we choose twists to evenly span the physical q 2 range for the maximum value of am val h on a given lattice, the greatest values of twist, θ, will correspond to negative values of q 2 for the smaller values of am val h used. We include these negative q 2 points in our analysis, since the form factors are analytic for −∞ < q 2 < M min pole where M min pole is the lowest mass subthreshold resonance.

A. Nonperturbative current renormalization
The renormalization factors relating our lattice current operators to the partially conserved lattice currents were computed in [11] for the axial-vector and in [12] for the vector current using the PCAC and PCVC relations respectively. These factors were not computed for am h ¼ 0.6 on set 1 and amh ¼ 0.65 on set 4, and for these we interpolate between the other values given in [11,12], noting that the differences between results on set 1 and set 4 are very small. In order to account for correlations FIG. 2. Arrangement of propagators in the three point function, we refer to c 1 and c 2 as the "spectator" and "active" charm quarks respectively and h as the "extended" heavy quark. J represents the insertion of a vector or axial-vector current.
TABLE II. Spin-taste operators used to isolate form factors. The first column is the operator used for the H c , the second for the J=ψ and the third column is the operator used at the current. between Z factors, which are neglected in our interpolation, we take the conservative approach of setting the uncertainty on our interpolated values equal to the largest uncertainty of the computed values on the corresponding set. We tabulate the renormalization factors in Table V, where we also include the ma-dependent discretization correction terms, Z disc , for the HISQ-quark tree level wave function renormalization computed in [34]. Z disc starts at ðam h Þ 4 because of the high level of improvement in the HISQ action, so is very close to 1. The total renormalization factor is given by Z V;A Z disc . Note that the pseudoscalar current used to determine A 0 is absolutely normalized and needs no Z factor. The statistical uncertainty of each lattice matrix element is typically at least an order of magnitude greater than that of the corresponding current renormalization factor. We therefore ignore correlations between different Z factors and between Z factors and our correlator data.

V. RESULTS
The correlator fits discussed in Sec. III were done simultaneously for all correlators on each set using the corrfitter PYTHON package [35]. This allows us to determine the ground-state fit parameters that we need along with a complete covariance matrix for them that can be fed into our continuum extrapolation fits. The priors for our correlator fits are listed in Table VI and take the same form on all sets. Our correlator fits include an SVD cut on the correlation matrix following the tools included in [35] (see the Appendix D of [36] for a discussion of this). We also truncate the time-length of the three point and two point correlators that we fit to remove excited-state dominated data points at early times and increase fit stability. The specifics of these fit choices are given in Table VII, along with variations used to test the stability of our analysis in Sec. V C.
As discussed in Sec. III we isolate form factors from the matrix elements extracted from correlation functions and combine these with current renormalization factors to obtain A 0 , A 1 , A 2 and V across a range of q 2 for several different heavy quark masses. The form factor results are given in Tables VIII-XI and plotted in Fig. 3. We see that the least well-determined form factor is A 2 , especially close to zero recoil. This can be traced back to the division by k 2 TABLE V. Z factors from [11,12] for the axial-vector and vector operators used in this work, together with the discretization corrections, Z disc . Z A and Z V values for am h ¼ 0. 6 where M max H c is the value of M H c corresponding to the largest am h , taken from [11]. While Ω J=ψ was chosen to follow the relativistic dispersion relation, Ω H c was chosen heuristically to give prior values approximately following the observed H c masses on each set while remaining suitably loose so as not to constrain the fit results.
Prior needed to determine it after the subtraction of the Oð1ÞA 1 contribution in Eq. (22). In Fig. 4 we plot J=ψ , which is the combination of A 2 and j ⃗ p 0 j 2 entering the helicity amplitudes in Eq. (5). We see that this has much clearer behaviour as a function of q 2 . TABLE VIII. Lattice form factor results for set 1. ak here is the value of the x and y components of the lattice momentum for the J=ψ. ak is calculated from the corresponding twist in Table IV. (12) (56) 0.8974(81) 1.17(74) 1.77(12) 0.141943 0.867 (43) 0.8875(81) 1.16(46) 1.743(91) 0.177429 0.849 (35) 0.8752(82) 1.15(34) 1.708(75) TABLE IX. Lattice form factor results for set 2. ak here is the value of the x and y components of the lattice momentum for the J=ψ. ak is calculated from the corresponding twist in Table IV.  X. Lattice form factor results for set 3. ak here is the value of the x and y components of the lattice momentum for the J=ψ. ak is calculated from the corresponding twist in Table IV. (25) (31) 0.52(12) 0.960(62) TABLE XI. Lattice form factor results for set 4. ak here is the value of the x and y components of the lattice momentum for the J=ψ. ak is calculated from the corresponding twist in Table IV.  (10) We also give values for the meson masses in Tables XII  for the J=ψ and Table XIII for the H c 1 and η h . From Table II we have three different J=ψ tastes, using one local and two different 1-link point-split operators.
The difference in mass of these different taste J=ψ mesons is tiny, barely visible above statistical uncertainties (at a level of 1 MeV) even on our coarsest lattices. For the H c we use two spin-taste operators, both local. As expected for pseudoscalars the taste-differences are slightly larger than for the vector J=ψ, and they are visible above the smaller statistical errors in this case. The mass differences are still only a few MeV, however, for a meson of a mass of several GeV, and they fall on the finer lattices since they are a discretization effect.
FIG. 3. The points show our lattice QCD results for each form factor as given in Tables VIII-XI as a function of squared 4-momentum transfer, q 2 . The legend gives the mapping between symbol colour and shape and the set of gluon field configurations used, as given by the lattice spacing, and the heavy quark in lattice units. The blue curve with error band is the result of our fit in lattice spacing and heavy quark mass, evaluating the result in the continuum limit and for the b quark mass, to give the physical form factor for B c → J=ψ. 1 Note that there is a systematic discrepancy between our values for aM H c ðγ 5 ⊗ γ 5 Þ on set 1, in Table XIII, and those in [11]. We found that increasing the statistical uncertainties on our masses on set 1 by a factor of 10 to cover this discrepancy did not change our final results significantly.
We calculate the value of q 2 for each form-factor based on the meson masses of the corresponding meson tastes, along with the lattice momentum of the J=ψ determined from the twists in Table IV.

A. Extrapolation to the physical point
In order to extract physical continuum form factors we must fit our lattice results to a fit function including discretization effects arising from the use of large values of am h as well as the physical dependence upon M H c and q 2 . In order to fit the q 2 dependence we make use of the z-expansion (see, for example, [37][38][39]) which has become standard for the description of form factors. We map the physical q 2 range to within the unit circle via the change of variables, The physical q 2 range is from 0 to t − where In Eq. (23) t þ is the pair production threshold for a bc current, and the start of a cut in the complex q 2 plane. This point is mapped to z ¼ −1. We take where H is an hū pseudoscalar meson. Note there is a sizeable gap between t 1=2 þ and t 1=2 − , and it is independent of m h up to binding energy effects. In this calculation we do not compute M H entering t þ . Instead we use an estimate of it derived from the M H c values that we have: This ensures that when we evaluate our fit function at the physical heavy quark mass then the correct physical continuum value of t þ is recovered. We take the D Ã mass from experiment since our valence c masses are tuned to the physical value throughout and the mistuning of the light quark masses in 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 with no uncertainty to simplify the necessary covariances required to reconstruct our fit function. Their uncertainties are negligibly small in any case.
The choice of t 0 in Eq. (23) determines the value of q 2 at which z ¼ 0. Here we choose t 0 ¼ t − . This choice means that z is simply related to the variable w which is the dot product of the 4-velocities of the B c and J=ψ and known to The combination of A 2 and j ⃗ p 0 j 2 entering the helicity amplitudes in Eq. (5). This combination does not exhibit as much noise as A 2 alone, seen in Fig. 3. TABLE XII. J=ψ masses for the local spin-taste operator γ 1 ⊗ γ 1 and 1-link operators γ 1 ⊗ 1 and γ 1 ⊗ γ 1 γ 2 used in our calculation, see Table II.   (17) be a good variable for a description of b → c form factors in the heavy quark effective theory approach [37,38] to B → D Ã decays. We expect that still to be a useful consideration here because the B c meson has similar properties to that of a heavy-light meson, interpolating between heavy-light and heavyonium [10]. Heavyonium mesons also display similar mass-suppressed differences in amplitudes between ground-state vector and pseudoscalar mesons, so that we expect a B c → J=ψ decay to behave rather similarly to B → D Ã . Physical particles withbc quark content, masses between t 1=2 þ and t 1=2 − (i.e., above the physical q 2 range for the decay) and the appropriate quantum numbers to couple to the current operator result in the appearance of a subthreshold pole in the corresponding form factor. Following [37,38] we include these poles in our fit using the form (Blaschke factor), Pðq 2 Þ is zero when q 2 ¼ M 2 pole and has modulus 1 along the cut in the q 2 plane for real q 2 > t þ where multiparticle production can occur. For the pole masses corresponding to B c states we use the values listed in [4], which we reproduce here in Table XV. We use four masses in each of the vector and axial channels and three in the pseudoscalar channel. These values are derived from experiment 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 setting up our fit in z-space the values are used without uncertainties. As we will show in Sec. V C our results are not sensitive to the exact form of the pole expression above [Eq. (26)]. The pole masses are shifted from those of the B c to account for our use of unphysical heavy-quark masses. We use again ensuring that when we extrapolate to the physical heavy quark mass the physical continuum values of Table XV are recovered. We fit each form factor, Fðq 2 Þ, to the fit function, where Pðq 2 Þ is the appropriate pole form for that form factor [constructed using 1 − states for Vðq 2 Þ, 1 þ states for A 1 ðq 2 Þ and A 2 ðq 2 Þ or 0 − states for A 0 ðq 2 Þ] as in Eq. (26). The remainder of the fit function is a polynomial in z with separate coefficients, a n , for each form factor that take the form, a n ¼ The Δ ðjÞ h allow for the dependence on the heavy quark mass using the η h mass as a physical proxy for this. We have Δ We take Λ ¼ 0.5 GeV. The polynomials in am c =π and am h =π in Eq. (28) allow for discretization effects that are either constant with the heavy quark mass (for example coming from the c-quark action) or that come from the heavy-quark action and so depend on am h . Because of the form of the HISQ action only even powers of a can appear here.
The remainder of Eq. (27), N n , takes into account the effect of mistuning the valence and sea quark masses for each form factor,   [40] except for the unphysical η s which we take from lattice calculations of the mass of the pion and kaon [41]. The η s mass is used to set mass mistuning terms in our fit and so include an uncertainty. The other masses are used as kinematic parameters in setting up our fit in z-space and used without uncertainties. The J=ψ mass is also used to tune the c quark mass, and there we do include its (negligible) uncertainty.  XV. Expected B c pseudoscalar, vector and axial vector masses below BD Ã threshold that we use in our pole factor, Eq. (26). Pseudoscalar values for the ground-state and first radial excitation are taken from experiment [40,[42][43][44]; the other values are taken from [4] and are derived from lattice QCD calculations [45] and model estimates [46][47][48]. We use the J=ψ meson mass to tune the c quark mass for the reasons discussed in [8]. We take We use the η s mass [41] to determine the mistuning of the s quark mass in the sea, taking We use the η s masses in lattice units from [11], together with the corresponding values of am val s , to do this. The values of the physical J=ψ and η s masses that we use are given, with their uncertainties, in Table XIV. To determine the mistuning of the u=d ¼ l quark mass in the sea we take with Δ½m s =m l ¼ 27.18ð10Þ from [49]. We take priors of 0(1) for each b n for each form factor, multiplying terms of order Oða 2 Þ by 0.5 because a 2 errors are removed in the HISQ action at tree-level [7]. We use priors of 0.0(0.1) for B n for each form-factor, following the results of the analysis of m sea c effects on w 0 in [31], and priors of 0.0(0.5) for C n and D n for each form-factor since sensitivity to u=d and s sea quark masses enter only at 1-loop. All remaining priors are taken as 0(1). We have checked that the prior width is conservative using the empirical Bayes criterion [50].
In doing our fit to Eq. (27) we impose the kinematical constraint, This condition holds in the continuum limit with physical sea quark masses for each heavy quark mass. We impose it using our lattice meson masses at each value of am h and allowing for discretization and quark mass mistuning effects effects in this implementation. We do this by imposing Δ kin here is a nuisance term made up of leading order discretization and mistuning effects to account for the use of lattice masses rather than values in the physical continuum limit. We take where α and β are priors taken as 0(1). We find that the fit returns values for α and β well within their prior widths. The fit was done simultaneously across all form factors in order to preserve correlations important for constructing helicity amplitudes. We use the nonlinear least squares fitting routine from the PYTHON package lsqfit [51]. It returns a χ 2 =d:o:f: of 0.18. The fit curves evaluated at the b quark mass in the continuum limit are plotted along with the raw lattice results in Fig. 3.
In Fig. 5 we show the lattice results and the fit function in z-space in a way that makes clearer how the fit is operating. The figure shows that the form factor multiplied by the pole function takes a very simple, linear, form in z-space. The constant term in the z-space polynomial varies with heavy quark mass but the slope changes very little.
The extrapolated continuum values of a n ¼ b 000 n are given in Table XVI. In the continuum limit at the physical b mass Eq. (27) becomes, for each form factor,  [37] includes additional outer functions along with the pole factor and allows unitarity constraints, which we do not make use of in our fit, to be applied. In Appendix B we convert our physical continuum results to the BGL scheme and check these unitarity constraints. We see that the bounds are far from saturation and so implementing them would have no effect on our results.

B. Heavy mass dependence
Our fits allow us to reconstruct the physical dependence of the form factors on the heavy quark mass. This allows us to check that the dependence is relatively benign and that we have adequate coverage of the heavy quark mass range. It could also be used to test theoretical understanding from model calculations.
This test of the physical heavy mass dependence is complicated slightly by the fact that our continuum fit results depend on m h both through the M η h dependence appearing in Δ ðiÞ h , and through M H c appearing in z, Pðq 2 Þ and our expression for the resonance masses. To evaluate our continuum form factors away from the physical b mass we must first determine the M η h dependence of M H c . To do this we fit our lattice data for the γ 5 ⊗ γ 5 H c and η h masses, and FIG. 6. The points show our lattice QCD results for each form factor as given in Tables VIII-XI as a function of the η h mass M η h , with data points corresponding to the same J=ψ spatial momentum (given in Tables VIII-XI) connected. We also use Eq. (42) to plot our continuum result (solid colored curves) at multiple, evenly spaced, fixed values of J=ψ momentum within the semileptonic region 0 ≤ q 2 ≤ q 2 max . The legend gives the mapping between symbol color and shape and the set of gluon field configurations used, as given by the lattice spacing, and the heavy quark in lattice units. Note that for the form factor A 2 we exclude from the plot the inaccurate lattice data for am h ¼ 0.5 on set 4. This form ensures the correct value of M H c as m h → m c . We take M phys η c ¼ 2.9839=GeV from [52] and neglect its very small uncertainty. We take prior widths of 0(1) for A 0 , B 0 , C 0 , D 0 , X i , Y i and Z i giving a good fit with χ 2 =d:o:f: ¼ 0.99 and Q ¼ 0.46. Our fitted parameters Z i may then be used to generate the continuum value of M H c at a given M η h by Note that this parametrization of the H c mass is only used here to demonstrate the heavy mass dependence of the form factors and will not have any impact on subsequent results. In Fig. 6 the form factors at fixed values of the J=ψ momentum are plotted against M η h . Here we choose values of the J=ψ momentum which evenly span the semileptonic range at the physical b quark mass and only plot the mass region for which the resulting value of q 2 is between 0 and q 2 max . We also plot our lattice data where we have connected points corresponding to a fixed J=ψ spatial momentum. Figure 6 illustrates that the continuum form factors have only mild heavy mass dependence across the range of masses used in this study. This demonstrates that we have good coverage of the heavy quark mass range and that our extrapolation to the b mass using these points is reliable. This is consistent with what is seen for other b → c form factors, e.g., [11,12].

C. Tests of the stability of the analysis
Our results for the form factors are dependent to some extent upon choices made in the correlator fits as well as choices made in the continuum/heavy-quark mass fit. Here we test the impact of those choices on both the continuum fit coefficients and the final outcome of our calculation to make sure that it is robust. It is convenient in that test to have a single quantity which we can compare, and we take that to be the total rate of B − c → J=ψl −ν l decay, i.e., ΓðB − c → J=ψμ −ν μ Þ=jη EW V cb j 2 . This is obtained by determining the helicity amplitudes from our form factors and then integrating in q 2 over the differential rate they give [see Eqs. (5) and (1)]. The results for the differential rates and total rate will be discussed in more detail in Sec. VI; here we focus on the stability of the final result under variations of fit choice.
We first look at the choices of the correlator fit parameters: ΔT 3pt , ΔT J=ψ 2pt , ΔT H c 2pt , the value of SVD cut and the number of exponentials used in the fit. In order to verify that our results are independent of such choices we repeat the full analysis using all combinations of the variations listed in Table VII. The total rate computed using each of these fit variations is plotted in Fig. 7, where we see that our final result is insensitive to such variations. Figure 8 shows a similar stability plot for a subset of the z-expansion coefficients for separate form factors (Eq. (38). We again see that the z-expansion coefficients are stable, within their uncertainty band, under changes in the correlator fits. Plots for other form factors not shown here are very similar.
We now turn to the effect of choices in the extrapolation to the physical point. The pole form of Eq. (26) removes q 2 dependence from poles in the q 2 plane above the physical region for the decay. For b → c decays these poles are substantially above the physical region. For example, here q 2 max is ð3.18 GeVÞ 2 whereas the lowest pole, for the pseudoscalar form factor is at the B c mass of 6.275 GeV. Hence we do not expect the exact positions or number of poles to have a big effect on the fits. The magnitude of Pðq 2 Þ does depend on how many poles are included [4], however, and this affects the normalization of the quantity form factor times Pðq 2 Þ that is expanded in z [Eq. (27)]. This in turn affects the prior width that must be allowed on the z-expansion coefficients to achieve a good fit.
We have investigated the effect of including fewer poles in Eq. (26) by repeating our analysis including only the first N poles resonances listed in Table XV. We then take a prior width on the z-expansion coefficients of 5.0 − N poles . We are able to obtain a good fit, with χ 2 =d:o:f: ≈ 0.2 in all cases. Since there are only three poles for A 0 expected below t þ , we include only three poles for that form factor even in the N poles case. Figure 9 shows these results, plotting against the lefthand y-axis the magnitude of the coefficient corresponding to the order z term, a 1 , coming from the fits as a function of the number of poles included. Results are given for each form factor. We see that as we include fewer poles, increasingly large z-expansion coefficients are needed FIG. 7. Plot showing the stability of the total rate for B − c → J=ψμ −ν μ under variations of the correlator fits. The x axis value corresponds to N ¼ δ 3 þ 4δ 2 þ 16δ 1 þ 64δ 4 where δ n is the value of δ corresponding to the fit given in Table VII for set n. The black horizontal line and red error band correspond to our final result, and the blue points and blue error band correspond to the combination of fit variations associated to N. Our result for the total rate is very stable to these variations. partly in order to account for the removal of physical q 2 dependence from missing poles but also because of the normalization change. Figure 9 also gives, marked as black crosses and using the right-hand y axis, the total width for B c → J=ψ obtained from that fit. We see that the total width is very stable as a function of N poles as the change in Pðq 2 Þ is compensated by the different coefficients obtained for the z-expansion.
We have also investigated the effects of including fewer z terms in Eq. (27) as well as fewer am c , am h and 2Λ QCD =M η h terms in Eq. (28). Figure 10 gives the total width obtained using these variations, where we see that our result is insensitive to the removal of the highest order terms. We also investigate the effect of increasing or decreasing the prior widths of the parameters b ijk n in Eq. (28) by a factor of 2. These results are also shown in Fig. 10 where we see only a very small effect on the central value of the total width. Note that the maximum number of poles included for A 0 is 3. The black crosses and error bars give the total width for the l ¼ μ case, Γ=jη EW V cb j 2 , determined from that fit, using the right-hand y axis. The grey band corresponds to our final result for the total width using N poles ¼ 4, and prior values for b ijk n of 0(1). This shows how the different coefficients as a function of N poles give a very stable result for the total width. Eq. (28). Oðn 1 ; n 2 ; n 3 ; n 4 Þ corresponds to the result including terms of highest order Oðð2Λ=M η h Þ n 1 ; ðam c Þ 2n 2 ; ðam h Þ 2n 3 ; z n 4 Þ. The vertical black line is our final result, corresponding to Oð3; 3; 3; 3Þ, and the grey band is its uncertainty. We also include variations in which we multiply our prior widths either by a factor of 2 or 0.5, labelled as 0ðσÞ → 0ðσ × 2Þ and 0ðσÞ → 0ðσ=2Þ respectively. Our result for the total rate is very stable to these variations.

VI. DISCUSSION
Using the form factors from Sec. VA we construct the helicity amplitudes using Eq. (5). The helicity amplitudes are plotted in Fig. 11, where we see that H 0 and H t are singular at q 2 ¼ 0 as we would expect from the 1= ffiffiffiffiffi q 2 p factors appearing in their definitions. This singular behavior is canceled by the ðq 2 − m 2 l Þ 2 factor appearing in the differential decay rate [Eq. (1)]. Together with the factor of j ⃗ p 0 j this results in physical rates which vanish at the extremes of the physical range, q 2 ¼ m 2 l and q 2 ¼ q 2 max . Using these helicity amplitudes we construct the differential decay rates using Eqs. (1), (2), (3) and (4). We take m μ ¼ 105.66 MeV and m e ¼ 0.51100 MeV [52] where we do not include their negligible uncertainties. The differential rates with respect to the angular variables defined in Fig. 1 are plotted in Fig. 12 and the differential rate dΓ=dq 2 is plotted in Fig. 13. In each case we plot the rate for l ¼ μ normalizing each differential rate by the total integrated decay rate. This means that any constant factors, such as jV cb j 2 or BðJ=ψ → μ þ μ − Þ cancel out. Where an integration over q 2 is necessary we use a simple trapezoidal interpolation in order to ensure covariances are carried through correctly, taking sufficiently many points that the results are insensitive to using additional points.
The angular differential rates show the expected qualitative picture for a pseudoscalar to vector meson decay [25]. The dependence on θ J=ψ is that for the case where the vector meson is seen in final states that carry spin (i.e., here in J=ψ → μ þ μ − ). This differs from, for example, 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 Þ 2 distribution in the W − rest frame. The picture switches for a c decay producing a virtual W þ [25].
We also compute the forward-backward asymmetry, A FB , the lepton polarization asymmetry (for the lepton l produced from the virtual W), A λ l , and the fraction of J=ψ produced with longitudinal polarization, F J=ψ L . These are defined as where we have chosen the forward direction for the purpose of A FB as being in the direction of the J=ψ momentum in the B c rest frame. These quantities are plotted in Fig. 14 for the cases l ¼ e and l ¼ μ. We see in the top plot that A FB is negative except near q 2 ¼ 0 for l ¼ μ. This may be be understood from Fig. 11 where we see that H 2 þ − H 2 − , which is the dominant contribution to A FB for m 2 l ≪ q 2 , is less than or equal to zero across the full physical q 2 range. The behavior of A FB near q 2 ¼ 0 in the l ¼ μ case originates from the −2m 2 l =q 2 H t H 0 cosðθ W Þ term in Eq. (2). When m 2 l =q 2 ≈ Oð1Þ it is apparent from Fig. 11 that this term will dominate over the H 2 þ − H 2 − contribution. In the middle plot of Fig. 14 we see that A λ e is equal to unity across the full q 2 range, in line with the expectation that in the massless limit the lepton is produced in a purely left-handed helicity eigenstate. In the bottom plot of F J=ψ L , approaches unity near q 2 ¼ 0 where H 0 and H t dominate the total rate, and goes to 1=3 at q 2 max where H 0 ¼ H þ ¼ H − and H t ¼ 0.
We also compute the total decay rate for the cases l ¼ e and l ¼ μ. We find and with the ratio Γ e =Γ μ ¼ 1.004018ð95Þ. We see that effects from including m μ amount to 0.4% of the rate. We can compare our results for the total rate to those from earlier, nonlattice QCD approaches such as potential models and QCD sum rules. In these other approaches it is much harder to quantify sources of uncertainty and derive an error budget. One way to obtain a picture of possible systematic errors is simply to compare results from different model calculations. Reference [53] uses a relativistic quark model to obtain a value 10.5 × 10 −12 GeV for Γ=jV cb j 2 and gives a table of comparison to other model approaches. Results vary over a range up to double this value, which implies a worrying lack of control of systematic effects in at least some of these calculations. In future it may be possible for such calculations to be tuned against our lattice QCD result for B c → J=ψ and then yield useful information on B c decay to excited charmonium which is much harder to determine accurately in lattice QCD.
From our value for Γ=jη EW V cb j 2 in Eq. (44) we can derive a result for the total width for B − c → J=ψμ −ν μ using a value for η EW (from Sec. II) and for jV cb j. We take jV cb j ¼ 41.0ð1.4Þ × 10 −3 from an average of inclusive and exclusive determinations and with the error scaled by 2.4 to allow for their inconsistency [52]. This gives Here the first uncertainty is from our lattice QCD calculation and the second from the uncertainty in jη EW V cb j. These uncertainties are approximately equal. Using the experimental average value for the B c lifetime of 0.510ð9Þ × 10 −12 s [52,54,55] then gives the branching fraction, The uncertainties are from lattice QCD, jη EW V cb j and the lifetime respectively. The accuracy of our result means that it can be used to normalize other B c decay modes. In the Particle Data Tables [52] B c branching fractions are given in the modified form Bðb → B þ c ÞΓ i =Γ. Here the first factor is the probability for a produced b quark to hadronise as a B c or one of its excitations. Using the value of the modified branching fraction for B − c → J=ψμ −ν μ obtained by CDF [56] of 8.7ð1.0Þ × 10 −5 [52] we infer The first uncertainty is the combined uncertainty from Eq. (47), and the second uncertainty is from the experimental value of the modified branching fraction. This value can be compared to the value Bðb → B þ Þ ¼ 0.407ð8Þ [52] to show how much less likely it is that a produced b quark will form a B c . We have focused here on the rates for production of a μ or e from the virtual W in the final state. We will discuss results for B − c → J=ψτ −ν τ in a follow-up paper [22].

A. Error budget
In order to estimate the contribution of the different sources of systematic uncertainty we use the inbuilt error budget routine in lsqfit [51], which computes the partial variance of our result with respect to the priors and data entering the fit (see e.g., [57] Appendix A). The error budget for the total rate for the l ¼ μ case is given in Table XVII. The size of the systematic uncertainties reflect the extent to which the corresponding fit parameters are determined by the prior values together with how much their values impact the value of Γ μ . Not surprisingly we see that large contributions come from the discretization effects from the heavy quark mass and the statistical uncertainties on the finest (ultrafine, set 3) lattices. The relatively large uncertainty coming from sea charm quark mass mistuining terms reflects the fact that we do not have access to configurations with multiple sea charm masses at a given lattice spacing, and hence this term is not well constrained. Excluding the sea charm mistuning term from our fits has only a small effect at of order 0.1σ.
Note that several of the largest uncertainties, those from statistics on set 3 and am h discretization effects and M η h dependence, may be straightforwardly and systematically reduced in the future by increasing statistics on ultrafine lattices (set 3) and obtaining results on exafine (a ¼ 0.03 fm) lattices respectively. Subleading errors may also be improved, with more precise determinations of w 0 or the inclusion of additional lattices with physical u=d quarks.

VII. CONCLUSIONS
We have carried out the first lattice QCD calculation of the form factors within the Standard Model for B − c → J=ψl −ν l decay. This is a process under active study by LHCb [1,58], and more accurate form factors are needed (and crucially with quantified errors) to reduce uncertainties in the experimental determination of the leptonuniversality violating ratio RðJ=ψÞ.
We give all four form factors across the full q 2 range of the decay. The calculation is done using the HISQ action for all quarks, enabling us to normalise the lattice current operators that couple to the W fully nonperturbatively. We have used the heavy-HISQ approach, obtaining results at multiple values of the heavy quark mass 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 parametrization that allows it to be fully reconstructed, including correlations between the form factors (see Appendix A). This yields a total rate integrated over q 2 for a light lepton in the final state, up to a CKM factor, of [repeating Eq. (44)] The 7% uncertainty is broken down in an error budget in The first two uncertainties, from our lattice QCD calculation and from the current uncertainty in V cb are equal here.
We will discuss the implications of our results for the case of a heavy (τ) lepton in the final-state elsewhere [22].
We have demonstrated that the heavy-HISQ methodology is a viable procedure for these calculations and that the statistical challenge of analyzing many lattice QCD correlation functions, while preserving covariances important for constructing and fitting q 2 and heavy mass dependence, can be handled. This work forms an important precursor for the calculation of B s → D Ã s form factors across the full q 2 range using the same machinery, which is underway. This will allow the dominant uncertainty from external inputs to be reduced in the determination of V cb from experimental results [59].

APPENDIX A: RECONSTRUCTING THE FIT
Our parametrization of the form factors for B c → J=ψ in the continuum limit is given in Eq. (38). It consists of a pole factor with no uncertainty and a polynomial in z for which the coefficients with their uncertainties are given in Table XVI. In this section we give the correlations between the z-expansion coefficients which are necessary for reconstructing our results explicitly, as well as instructions for using the included in Supplemental Material [60] to load the z-expansion parameters together with their correlations automatically into PYTHON [61].
The correlation between two coefficients is defined in the usual way as 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: Extrapolation of data in M η h to the physical value of M η b , errors proportional to powers of am c , am h , mass mistuning effects δ m , and uncertainties in determination of the lattice spacing. The second section gives a breakdown of the statistical uncertainties corresponding to each of our data sets. Finally "other priors" includes all of the remaining sources of uncertainty such as Δ kin and the current renormalization factors Z. where σ 2 ðXÞ is the variance of X andX is the mean of X.

Source
The values are tabulated in Tables XVIII to XXVII.
In this calculation and in the Supplemental Material [60] we use the GVAR PYTHON package to track and propagate correlations. Included in the Supplemental Material [60] are two text files; CORRELATIONS.txt contains a dictionary including the means and variances of the z-expansion       parameters on the first line and a dictionary detailing the correlations between these parameters on the second line, CHECKS.txt contains arrays of q 2 values and form factor mean and standard deviation values at the corresponding values of q 2 . This file is used by the PYTHON script load_fit.py as a simple check that the fit has been loaded correctly. Running python load_fit.py will load the parameters from CORRELATIONS.txt and compare values computed at hard coded intervals in q 2 to those in CHECKS.txt which were computed as part of this work. Running python load_fit.py will also produce some simple plots of the form factors across the full q 2 range. We have tested load_fit.py using PYTHON [64].

APPENDIX B: CONVERTING TO THE BGL PARAMETRIZATION AND TESTS OF UNITARITY BOUNDS
The now standard BGL parametrization [37] makes use of the form factors in the helicity basis. These are given, in terms of the form factors used in this work, by [65], where ⃗ p 0 is the J=ψ spatial momentum in the B c rest frame. The BGL scheme then parametrizes these form factors using the expansion in z space, where the pole function P i is the same as we have defined in Eq. (26) and the outer functions, ϕ, are defined in [65]. Note that here we use the subthreshold resonance masses given in Table 2 of [65] in the pole functions P i . In order to compute the outer functions we evaluate χ LðTÞ ðAEuÞ using the expressions in [37], using the pole mass m b ¼ 4.78 GeV from [52] and α s ¼ 0.22 as in [37]. Their numerical values are given in Table XXVIII. In order to convert our results into this parametrization scheme we use our fitted continuum parametrization to output form factor values at a large number of q 2 values in the semileptonic region q 2 < t − and then fit these using Eq. (B5) truncated at n ¼ 3. We tabulate these values, together with P 3 n¼0 ða BGL n Þ 2 , in Table XXIX where we see that the unitarity bound P 3 n¼0 ða BGL n Þ 2 ≤ 1 is satisfied comfortably for each form factor.