Model-Independent Bounds on $R(J/\psi)$

We present a model-independent bound on $R(J/\psi) \! \equiv \! \mathcal{BR} (B_c^+ \rightarrow J/\psi \, \tau^+\nu_\tau)/ \mathcal{BR} (B_c^+ \rightarrow J/\psi \, \mu^+\nu_\mu)$. This bound is constructed by constraining the form factors through a combination of dispersive relations, heavy-quark relations at zero-recoil, and the limited existing determinations from lattice QCD. The resulting 95\% confidence-level bound, $0.20\leq R(J/\psi)\leq0.39$, agrees with the recent LHCb result at $1.3 \, \sigma$, and rules out some previously suggested model form factors.


Introduction
Within the Standard Model, lepton universality is broken only by the Higgs interaction, but the discovery of neutrino masses implies that at least one, and potentially several, relevant forms of beyond-Standard Model modification exist. The ratios of semileptonic heavymeson decay branching fractions to distinct lepton flavors represent a group of observables particularly sensitive to new physics, because the QCD dynamics of the heavy-meson decays decouples from the electroweak interaction at leading order: This expression implies that the ratios of semileptonic heavy-meson decay branching fractions can differ from unity at this level of precision only due to kinematic factors. Measurements from BaBar, Belle, and LHCb of the ratios R(D ( * ) ) for heavy-light meson decays B → D ( * ) ν with = τ to those with = µ or e (or their average) exhibit tension with theoretical predictions. The HFLAV averages [1] of the experimental results R(D * ) = 0.306(13)(7) [2][3][4][5][6][7][8][9][10] and R(D) = 0.407(39)(24) [2][3][4] represent a combined 3.8σ discrepancy [1] from the HFLAVsuggested Standard-Model value of R(D * ) = 0.258(5) [1] obtained by an averaging [11][12][13] that utilizes experimental data, lattice QCD results, and heavy-quark effective theory, and from R(D) = 0.300(8) [14], which is an average of lattice QCD results [15,16], as well as a value R(D) = 0.299(3) obtained by also including experimental data supplemented by heavy-quark effective theory [17]. In light of this tension, the LHCb Collaboration has measured the rates for the heavy-heavy semileptonic meson decays B + c → J/ψ + ν ( Fig. 1) in the = τ, µ channels, finding R(J/ψ) = 0.71 (17) (18) [18].
At present, only model-dependent calculations of R(J/ψ) exist (collected in Table 1) [19][20][21][22][23][24][25][26][27][28][29][30][31][32]. Although most models' central values cluster in LHCb's quoted theory range of 0.25-0.28, one notes a wide spread in their estimated uncertainty. We take as a reasonable estimate of the model range 0 < R(J/ψ) < 0.48, the union of the 95% confidence levels (CL) of the reported theoretical uncertainties, which in turn typically account only for parameter fitting. These results rely upon approximations such as nonrelativistic reduction, constituent quarks, or perturbative QCD to obtain transition form factors between the heavy-heavy B + c and J/ψ mesons. Without a clear understanding of the systematic uncertainties these assumptions introduce, the reliability of these predictions is suspect. In this paper, we present the first model-independent constraint, a 95% CL bound of 0.20 ≤ R(J/ψ) ≤ 0.39 within the Standard Model, in which uncertainties are all quantifiable. In order to obtain this result, we begin in Sec. 2 with a discussion of the V −A structure of the Standard Model and the form factors. In Sec. 3 we explain how heavy-quark spin symmetry can be applied at the zero-recoil point to relate the form factors, using the method of [20]. The initial lattice-QCD results of the HPQCD Collaboration [33,34] for two of the transition form factors are discussed in Sec. 4. The dispersive analysis framework utilized to constrain the form factors as functions of momentum transfer is presented in Sec. 5. The results of our analysis, as well as future projections for the bound, appear in Sec. 6, and we conclude in Sec. 7.
In the Standard Model, the factorization of Eq. (1.1) into a leptonic and a hadronic tensor reduces the problem of calculating R(J/ψ) to the computation of the hadronic matrix element J/ψ |(V −A) µ |B + c . Using this factorization, the hadronic matrix element can be written in terms of four transition form factors. These form factors enter the matrix element in combination with the meson masses, M ≡ M B + c and m ≡ M J/ψ , the corresponding meson momenta P µ and p µ , and the polarization µ of the J/ψ. The form factors themselves depend only upon q 2 = (P − p) 2 , the squared momentum transfer to the leptons. A number of form-factor decompositions exist in the literature; one common set [35] used in lattice-QCD [33] and model calculations is given by V (q 2 ), A i (q 2 ), i = 0, 1, 2, 3: where q µ ≡ (P −p) µ . While we have exhibited five form factors, only four are independent. In the physical set, A 0 (q 2 ) is defined as the unique form factor that couples to timelike virtual W polarizations (∝ q µ ), while A 3 (q 2 ) is simply a convenient shorthand for a combination appearing in intermediate stages of calculations, and in fact satisfies Furthermore, the finiteness of Eq. (2.1) as q 2 → 0 requires A 3 (0) = A 0 (0), which proves useful in constructing our bounds. In what follows, we also use the notation t ≡ q 2 , and define two important kinematic points, t ± = (M ± m) 2 . Using Eq. (2.1) or an equivalent basis, model predictions include uncontrolled approximations for the form factors. Some models construct wave functions for the two mesons, while others attempt to compute a perturbative distribution amplitude at q 2 → 0 and then extrapolate to larger values with some functional form. In addition, some models do not respect form-factor relations, such as the heavy-quark spin-symmetry relations discussed below. Due to these issues, the good agreement seen between the model predictions may more reflect the theoretical prejudice in modeling than a genuine estimate of the true Standard-Model value.
While this decomposition is useful for lattice QCD, it is not the best decomposition for the dispersive analysis. The second convention we use is the helicity basis, which exchanges the form factors V, A i for g, f , F 1 , and F 2 . 1 They are related by 1 Strictly speaking, F1,2 are helicity amplitudes (in conventional notation [36], proportional to H0,t, respectively), while f, g are two linear combinations of them: where k is defined in Eq. (2.4). Table 1. Model predictions of R(J/ψ) classified by method, which are abbreviated as: constituent quark model (CQM), relativistic quark model (RCQM), QCD sum rules (QCDSR), nonrelativistic quark model (NRQM), nonrelativistic QCD (NRQCD), and perturbative QCD calculations (pQCD).

Heavy-Quark Spin Symmetry
Decays of heavy-light Qq systems possess enhanced symmetries in the heavy-quark limit because operators that distinguish between heavy quarks of different spin and flavor are suppressed by 1/m Q , and their matrix elements vanish when m Q → ∞. Consequently, all transition form factors Q q |O|Qq in this limit are proportional to a single, universal Isgur-Wise function ξ(w) [37,38], whose momentum-transfer argument is w, the dot product of the initial and final heavy-light hadron 4-velocities, v µ ≡ p µ M /M and v µ ≡ p µ m /m, respectively: At the zero-recoil point t = (M −m) 2 or w = 1, the daughter hadron m is at rest with respect to the parent M . Indeed, one notes that w equals the Lorentz factor γ m of m in the M rest frame. The maximum value of w corresponds to the minimum momentum transfer t through the virtual W to the lepton pair, which occurs when the leptons are created with minimal energy, t = m 2 .
In heavy-light systems, the heavy-quark approximation corresponds to a light quark bound in a nearly static spin-independent color field. In the weak decay Q → Q between two very heavy quark flavors, as the momentum transfer to the light quark t → 0, q no longer changes states, and therefore the wave function of this light spectator quark remains unaffected. One thus concludes that ξ(1) = 1 at the zero-recoil (Isgur-Wise) point, yielding a absolute normalization for the form factors. These results are accurate up to corrections of O(Λ QCD /m Q ).
In the decay B + c → J/ψ, the spectator light quark is replaced by another heavy quark, c. This substitution results in a pair of related effects on the enhanced symmetries of the heavy-quark limit [39]. First, the difference between the heavy-quark kinetic energy operators produces energies no longer negligible compared to those of the spectator c, and this effect spoils the flavor symmetry in heavy-heavy systems. Furthermore, the spectator c receives a momentum transfer from the decay ofb →c of the same order as the momentum imparted to thec, so one cannot justify a normalization of the form factors at the zero-recoil point based purely upon symmetry.
While the heavy-flavor symmetry is lost, the separate spin symmetries ofb andc quarks remain, with an additional spin symmetry from the heavy spectator c. Furthermore, the presence of the heavy c suggests a system that is closer to a nonrelativistic limit than heavy-light systems. In the B + c → J/ψ semileptonic decays, one further finds that suggesting that an expansion about the zero-recoil point may still be reasonable. Together, the spin symmetries imply that the four form factors are related to a single, universal function h (∆ in Ref. [39]), but only at the zero-recoil point, and no symmetry-based normalization for h can be derived [39].
Using the trace formalism of [40], in Ref. [39] it was shown how to compute the relative normalization between the fourQq →Q q form factors near the zero-recoil point [i.e., where the spatial momentum transfer to the spectator q is O(m q )]. Using these relations, h was derived for a color-Coulomb potential in Ref. [39]. This approximation was improved in Ref. [41], where a constituent quark-model calculation of BR(B + c → J/ψ + ν ) for = e, µ but not τ , was performed. The heavy-quark spin-symmetry relations of Ref. [39] were generalized in [20] to account for a momentum transfer to the spectator quark occurring at leading order in NRQCD, specifically, to the case v = v but w → 1. We reproduce here the relations of [20], where the form factors g(w = 1), F 1 (w = 1), F 2 (w = 1) are related to f (w = 1) by These relations reproduce the standard Isgur-Wise result [37,38,42] when σ = 0. The relation between F 1 (w = 1) and f (w = 1) follows directly from the definition of Eq. (2.3), independent of heavy-quark symmetries.
Terms that break these relations should be O(m c /m b , Λ QCD /m c ) ≈ 30%, and we allow conservatively for up to 50% violations. The heavy-quark spin symmetry further relates the zero-recoil form factors of B + c → J/ψ to those of B + c → η c , which will be useful in the future to obtain further constraints.
In analogy with the heavy-light systems, we can enforce a further constraint from heavy-quark symmetries. The universal form factor h represents the overlap element of the initial and final states, and therefore should be maximized at w = 1. This statement is an assumption, but a very mild one: In the heavy-light system, the slope of the Isgur-Wise function is rigorously negative at w = 1 [43][44][45], and it would indeed be very surprising if the same did not hold for the form factors of heavy-heavy mesons, which are even more similar to idealized quark-model states.

Lattice QCD Results
The state-of-the-art lattice QCD calculations for B + c → J/ψ are limited to preliminary results from the HPQCD Collaboration for V (q 2 ) at two q 2 values and A 1 (q 2 ) at three q 2 values [33,34]. These results were obtained using 2+1+1 HISQ ensembles, in which the smallest lattice spacing is a ≈ 0.09 fm, and the b quark is treated via NRQCD, and are reproduced in Fig. 2. For q 2 = t − , 0 A 1 (q 2 ) has also been computed on coarser lattices and for lighter dynamical b-quark ensembles, in order to check the accuracy and assess the uncertainty of the a ≈ 0.09 fm NRQCD results. At present, there are no lattice results . Below, we show that the most desirable piece of new information from the lattice is a computation of A 0 (0), which could cut our uncertainties in half.

Dispersive Relations
In this work we derive constraints on the form factors of B + c → J/ψ using analyticity and unitarity constraints on a particular two-point Green's function and a conformal parameterization in the manner implemented by Boyd, Grinstein, and Lebed (BGL) [46] for the decays of heavy-light hadrons to heavy-light or light-light hadrons. We utilize a slightly different set of free parameters to simplify the computation for our particular case of a heavy-heavy meson decaying to another heavy-heavy meson. Here we briefly sketch the necessary components, emphasizing where we differ from the literature.
To derive our constraints, one considers the two-point momentum-space Green's function Π µν J of a vectorlike quark current, J µ ≡QΓ µ Q . Π µν J can be decomposed in different ways [42,[47][48][49][50]; in this work we choose to separate it into spin-1 (Π T J ) and spin-0 (Π L J ) piecesà la [42] via From perturbative QCD (pQCD), the functions Π L,T J are known to contain first-and secondorder divergences, respectively, and must undergo subtractions in order to be rendered finite. The finite dispersion relations are: The freedom to chose a value of q 2 can be leveraged to compute χ(q 2 ) reliably in pQCD, far in q 2 from where the two-point function receives nonperturbative contributions from effects such as bound states and resonances. The formal condition on q 2 to be in the perturbative regime is which, for Q, Q = c, b, q 2 = 0 is clearly sufficient. Existing calculations of two-loop pQCD χ(q 2 = 0) modified by non-perturbative vacuum contributions [51][52][53][54][55] used in Ref. [42] can be applied here. An example of the state of the art in this regard (although slightly different from the approach used here) appears in Ref. [17]. The spectral functions Im Π J can be decomposed into a sum over the complete set of states X that can couple the current J µ to the vacuum: Each term in the sum is semipositive definite, thereby producing a strict inequality for each X in Eqs. (5.2). These inequalities can be made stronger by including multiple X at once, as discussed in Refs. [12,13,42]. For X we include only below-threshold B + c poles and a single two-body channel, B + c +J/ψ, implying that our results provide very conservative bounds.
In contrast to many prior dispersive analyses, B + c → J/ψ, like the Λ b → Λ c process studied in Ref. [50], does not give the lightest two-body threshold with the correct quantum numbers; these lighter thresholds must be taken into consideration. Depending upon the quantum numbers indicated by J, the first physically prominent two-body production threshold in t occurs at B ( * ) +D (see Table 2). In early literature such as [50], the branch cut starting at the threshold for the process of interest was the one used in the dispersive analysis, while the effect of the cut from the lower threshold up to this threshold was modeled and argued to amount to a slight loosening of the unitarity bound given below by Eq. (5.11). Here, however, we represent the analytic features more faithfully by using the lower threshold directly. With this fact in mind, we define a new variable t bd ≡ (M B ( * )+MD ) 2 that corresponds to the first branch point in a given two-point function, while the B + c +J/ψ branch point occurs at t + > t bd .
With these variables, one maps the complex t plane to the unit disk in a variable z (with the two sides of the branch cut forming the unit circle C) using the conformal variable transformation where t * is the branch point around which one deforms the contour, and t 0 is a free parameter used to improve the convergence of functions at small z. In this mapping, z is real for t ≤ t * and a pure phase for t ≥ t * . Prior work that computed the form factors between baryons whose threshold was above that of the lightest pair in that channel (i.e., Λ b → Λ c , Λ b → p) took t * = t + [42,49], which introduces into the region |z| < 1 a subthreshold branch cut, meaning that the form factors have complex nonanalyticities that cannot trivially be removed. To avoid this issue, we instead set t * = t bd , which is possible because we are only interested in the semileptonic decay region, m 2 ≤ t ≤ t − , which is always smaller than t bd . This choice ensures that the only nonanalytic features within the unit circle |z| = 1 are simple poles corresponding to single particles B ( * )+ c , which can be removed by Blaschke factors described below. The need to avoid branch cuts but not poles from |z| < 1 derives from the unique feature of the Blaschke factors, which can remove each pole given only its location (i.e., mass), independent of its residue. 2 In contrast, correctly accounting for a branch cut requires knowledge of both the location of the branch point and the function along the cut.
To remove these subthreshold poles, one multiplies by z(t; t s ) [using the definition of Eq. (5.5)], a Blaschke factor, which eliminates a simple pole t = t s . Using this formalism, the bound on each form factor F i (t) can be written as The function P i (t) in Eq. (5.6) is a product of Blaschke factors z(t; t p ) that remove dynamical singularities due to the presence of subthreshold resonant poles. Masses corresponding to the poles that must be removed in B + c → J/ψ are found in Table 2, organized by the channel to which each one contributes. These masses have either been measured by LHCb [58,59] or derived from model calculations [60], with uncertainties that are negligible for our purposes.
The weight function φ i (t; t 0 ) is called an outer function in complex analysis, and is given by where j = T, L (for which n j = 3, 2, respectively), the functionP i (t) is a product of factors z(t; t s ) or z(t; t s ) designed to remove kinematical singularities at points t = t s < t bc from the other factors in Eq. (5.6), and W i (t) is computable weight function depending upon the particular form factor F i . The outer function can be reexpressed in a general form for any particular F i as where n I is an isospin Clebsch-Gordan factor, which is 1 for B + c → J/ψ. The remaining factors are found in Table 3. Transforming the dispersion-relation inequality based upon Eq. (5.4) into z-space, Eq. (5.6) becomes which, upon dividing out the non-analytic terms, allows the expansion in z corresponding to an analytic function: (5.10) Inserting this form into Eq. (5.9), one finds that the bound can be compactly written as a constraint on the Taylor series coefficients: All possible functional dependences of the form factor F i (t) consistent with Eqs. (5.2) are now incorporated into the coefficients a in . It is useful to introduce a number of dimensionless parameters that are functions of the meson masses: and a parameter N related to t 0 in Eq. (5.5) by It is straightforward to compute the kinematical range for the semileptonic process given in terms of z: 14) The minimal (optimized) truncation error is achieved when z min = −z max , which occurs when N opt = λ κ . (5.15) Evaluating at N = N opt , one finds From these expressions, we find that the semileptonic decays have z max,τ ≈ 0.019 and z max,µ ≈ 0.027, where each has a 1.5% variation, depending upon whether the BD or B * D threshold is the lowest branch point, t bd .

Results
Before presenting our bound on R(J/ψ), we summarize the constraints the form factors g, f, F 1 , F 2 are required to satisfy: • The coefficients a n of each form factor are constrained by n a 2 n ≤ 1 [Eq. (5.11)], in particular, for the cases n = 1, 2 investigated here.
• Using Eq. (3.3), the values g(t − ) and F 2 (t − ) are related to the value of f (t − ), which in turn is computed from lattice QCD, to within 50%.
• All form factors (which are defined to have the same sign convention as the Isgur-Wise function) are assumed maximal at the zero-recoil point t = t − since the universal form factor h represents an overlap matrix element between initial and final states. Although this condition is not required by the model-independent parametrization Eq. (5.10), it appears to be supported by all the models cited in Table 1 for which functional expressions of form factors are provided. We find this condition to be suitably implemented via the constraints F i (t − ) ≥ F i (0) and dF i dt t − ≥ 0, where F i represents any of the form factors.
Imposing these constraints, we perform our fit in two steps, reflecting the difference in information between the two form factors (V, A 1 ) for which lattice values have been computed, and the two (A 0 , A 2 ) without.
In the first step, random Gaussian-distributed points are sampled for the form factors g and f [equivalently, by Eq. (2.3), V and A 1 ] whose mean gives the HPQCD results. The combined uncertainties are given by the quadrature sum of the reported uncertainty δ lat of the form-factor points and an additional systematic uncertainty, f lat (expressed as a percentage of the form-factor point value) that we use to estimate the uncomputed lattice uncertainties (i.e., finite-volume corrections, quark-mass dependence, discretization errors). f lat is taken to be 1, 5, or 20% of the value of the form factor from the lattice. For our final result, we suggest using f lat = 20%, while the other two values are helpful for understanding future prospects with improved lattice data. Using these sample points, we compute lines of best fit, from which we produce the coefficients a n . The resulting bands of allowed form factors are shown for f lat = 20% in Fig. 2, alongside the HPQCD results.
In the second step, we compute F 1 and F 2 (which include A 0,2 ), for which no lattice information exists. One could adopt the tactic of randomly generating points with some prior distribution, which, once accounting for the constraints, could be used to suggest a mean value of R(J/ψ) with some prior-dependent uncertainty. We instead opt to remove this possible dependence by obtaining the numerical maximum and minimum R(J/ψ) values, subject to the computed f, g values and the constraints listed above. In this way, the only uncertainties included are those from the lattice-QCD results and the violations of the heavy-quark spin-symmetry relations. The resulting bands of form factors for F 1 and  Figure 2. B + c → J/ψ form factors V (q 2 ) (red circles) and A 1 (q 2 ) (blue triangles) from the HPQCD Collaboration [33,34]. The interior bars represent the statistical uncertainty quoted by HPQCD. The exterior bars represent the result of including our f lat = 20% systematic uncertainty. The colored bands DA (dispersive analysis) represent our one-standard-deviation (1σ) best-fit region. that produce the minimum and maximum values of R(J/ψ) subject to the constraints are plotted in Fig. 3.
Having computed all the form factors, we present the 95% CL ranges for R(J/ψ) as a function of the truncation power n = 1, 2 in the dispersive analysis coefficients of Eq. (5.10) and the 1, 5, 20% systematic uncertainty f lat associated with the lattice data. The full results are presented in Table 4. The bound on R(J/ψ) appears relatively insensitive to increasing the number of free parameters a n , and only mildly dependent upon f lat . One might be concerned that increasing n could dramatically change these results, but we note that the typical value of n a 2 n for n = 1 is O(10 −2 ), while for n = 2 we find a 2 ≈ 1. While the dispersive constraint is saturated in the n = 2 case, the bound on R(J/ψ) is only enlarged by 5%. However, the saturation of any particular a n is not necessary to find the effect of higher a n to be negligible. Since higher-order terms are suppressed by z max ≈ 0.03, in order for these terms to contribute strongly, one must have a n+1 z max a n . Such an a n+1 value would either have to violate saturation n a 2 n ≤ 1 once the lower-order terms a n are fixed, or else it would change the numerical results very little.
In Fig. 4 we plot the previous model-dependent values of R(J/ψ) alongside the LHCb result and our 95% CL bound of 0.20 ≤ R(J/ψ) ≤ 0.39, as a function of publication date.  One can see that, while many of the previous model results lie within our 95% CL band, some are either partially or entirely excluded. The anomalously low NRQCD result of Ref. [25] is in severe disagreement with our bound (the small R(J/ψ) of Ref. [25] can be attributed to a larger-than-typical muonic branching ratio), while all the other models that have included 1σ theory uncertainty estimates are seen to remain compatible with the dispersive bounds obtained from a fairly sparse set of lattice results. Our 95% CL band should be viewed as a model-independent upper limit (subject to the assumptions listed above) for the largest theory uncertainty any model can allow and remain consistent with analyticity, unitarity, and existing lattice "data." With better lattice results-and/or actual experimental measurements of the form factors at any values of q 2 -the allowed parameter space for any given model will become severely curtailed.   Table 1.
Considering that the fits to the lattice data already constrain the form factors g and f well, we can ask what new piece of information would most improve our bounds. Since F 2 essentially affects only the τ channel, it is obvious that this form factor is the one most important to reducing the range in R(J/ψ). By inspection of both F 1 and F 2 in Fig. 3, we find that the q 2 dependence of the form factors generating the maximum and minimum R(J/ψ) values are quite different in shape. The minimum-R(J/ψ) ones clearly prefer flat form factors, while the maximum ones are nearly zero at q 2 = 0 and rise to their maximum allowed values at q 2 = t − . This dependence suggests that obtaining a value of F 2 (q 2 = 0) would greatly improve the lower bound, while F 2 (t − ) would restrict the upper bound. With the LHCb result already lying slightly above our bound, it would be most incisive to reduce our upper bound. This zero-recoil form factor is directly related by F 2 (t − ) = 2A 0 (t − ) to a traditional lattice form factor, and therefore should be possible to compute.
To investigate the possible effect of this new information, we consider a synthetic point F 2 (t − ) = 2(1±f lat ). This particular value is chosen because it lies near the average of the minimum and maximum preferred values, and is similar to the values suggested in models.
Taking f lat = 20%, we find that the bound could tightened to [0. 20, 0.35]. If this point and the existing 5 lattice points reached f lat = 1%, one could anticipate a range of [0.21, 0.32]. So, an additional lattice point at F 2 (t − ) could improve the bound by the same amount as reducing the uncertainty f lat from 20% to 1% (as seen in Table 4), but with substantially less computing resources.

Discussion and Conclusion
In contrast to the model-dependent previous works, we have presented a model-independent bound on R(J/ψ), finding it constrained to lie in the range 0.20 ≤ R(J/ψ) ≤ 0.39 at the 95% CL. At this level, we find that the LHCb result is consistent with the Standard Model at 1.3 σ. The near-term outlook for a higher-statistics LHCb measurement, coupled with new lattice results, promises to reduce the uncertainty on the experimental and theoretical values dramatically.
Even without a lattice QCD calculation of the F 2 form factor, additional potential areas of improvement can be investigated. Experience in the heavy-light sector and the fact that the R(J/ψ) bounds require saturating a 2 n = 1 suggest that including multiple states appearing in the dispersion relation can provide complementary information to help constrain the form factors further, and in this case one can additionally include the lattice results for B → D ( * ) [15,16,[61][62][63][64] and Λ b → Λ c [65].