Model-independent prediction of R(ηc)

We present a model-independent prediction for R(ηc)≡ℬℛ(Bc+ → ηcτ+ντ)/ℬℛ(Bc+ → ηcµ+νµ). This prediction is obtained from 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 prediction, R(ηc) = 0.29(5), agrees with previous model predictions, but without uncontrolled systematic uncertainties.


Introduction
The Higgs interaction is the only source of lepton universality violations within the standard model, but the observation of neutrino masses implies that at least one form of beyondstandard model modification exist. The ratios of semileptonic heavy-meson decays for distinct lepton flavors are particularly sensitive to new physics, because the QCD dynamics of the heavy-meson decays decouple from the electroweak interaction at leading order: This expression implies that the ratios of semileptonic heavy-meson decays can differ from unity at this level of precision only due to kinematic factors, although it is possible to further remove this dependence [1][2][3][4][5][6][7][8]. Measurements from BaBar, Belle, and LHCb of the ratios R(D ( * ) ) of heavy-light meson decays B → D ( * ) ν, with = τ to = µ, exhibit tension with theoretical predictions. The HFLAV averages [9] of the experimental results R(D * ) = 0.306(13)(7) [10][11][12][13][14][15][16][17][18] and R(D) = 0.407(39)(24) [10][11][12] represent a combined 3.8σ discrepancy [9] from the HFLAV-suggested Standard-Model value of R(D * ) = 0.258(5) [9] obtained by an averaging [7,19,20] that utilizes experimental form factors, lattice QCD results, and heavy-quark effective theory, and from R(D) = 0.300(8) [21], which is an average of lattice QCD results [22,23], as well as a value R(D) = 0.299(3) obtained by also including experimentally extracted form factors [24]. Recently, the LHCb collaboration has measured R(J/ψ) = 0.71 (17) (18) [25] which agrees with the Standard-Model bound of 0.20 ≤ R(J/ψ) ≤ 0.39 at 1.3σ [26]. In the future, it would be useful to consider thebc →cc analog of the B → D process, B + c → η c . Alas, measurements of R(η c ) are substantially harder than R(J/ψ) for a few reasons, foremost of which is there is no clean process like J/ψ → µ + µ − in which to reconstruct the η c , which will result in larger backgrounds. Additionally the transition to η c from excited states is poorly understood, and this further complicates extraction of signals [27]. Despite these present experimental difficulties, it would be valuable to have a theoretical prediction for R(η c ) from the Standard Model ready for it. The current state of affairs, though, is limited to model-dependent calculations (collected in table 1) [3,[28][29][30][31][32][33][34][35][36][37][38][39]. Although most models' central values cluster in the range 0.25-0.35, one notes a wide spread in their estimated uncertainty which typically account only for parameter fitting. These results rely upon some approximations to obtain the B + c → η c transition form factors. Without a clear understanding of the systematic uncertainties these assumptions introduce, the reliability of these predictions is suspect.
We fill a blank space in the literature by computing a model-independent prediction, R(η c ) = 0.29(5) from the Standard Model, in which all uncertainties are quantifiable. In order to obtain this result, we begin in section 2 with a discussion of the V − A structure of the Standard Model and the form factors. In section 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 [29]. The initial lattice QCD results of the HPQCD collaboration [40,41] for the transition form factors are discussed in section 4. The dispersive analysis framework utilized to constrain the form factors as functions of momentum transfer is presented in section 5. The results of our analysis, as well as future projections, appear in section 6, and we conclude in section 7.
After this calculation was completed, a similar calculation appeared [42] that is in good agreement with ours.
In the Standard Model, the factorization of eq. (1.1) into a leptonic and a hadronic tensor reduces the problem of calculating R(η c ) to the computation of the hadronic matrix where the axial component goes to zero due to parity. Using this factorization, the hadronic matrix element can be written in terms of two transition form factors. These form factors enter the matrix element in combination with the meson masses, M ≡ M B + c and m ≡ M ηc , and the corresponding meson momenta P µ and p µ . The form factors themselves depend only upon t ≡ q 2 = (P − p) 2 , the squared momentum transfer to the leptons. The hadronic matrix element in our convention is given by f + (t), f − (t): In this work, we will exchange f − for f 0 , which is given by We further introduce two important kinematic values t ± = (M ± m) 2 . This convention differs from that utilize by HPQCD for their lattice QCD results [40,41] by the mass dimension of f 0 . The conversion between the two is In the HPQCD convention, it can be seen that f HPQCD (0), which should be applied when fitting the functions. Using eq. (2.1) or an equivalent basis, form factors are computed from models with uncontrolled approximations. Some models construct wave JHEP12(2018)114 functions for the two mesons, while others compute a perturbative distribution amplitude at q 2 → 0 and then extrapolate to larger values. In addition, some models violate delicate form-factor relations, such as the heavy-quark spin-symmetry relations discussed below. Due to these issues, it is potentially treacherous to take the all too well agreement seen between the model predictions as a genuine estimate of the true standard model value instead of a theoretical prejudice in modeling.
The differential cross section for the semileptonic decay is where, in terms of the spatial momentum p of the η c in the B + c rest frame, Inspecting eq. (2.4), one can see that in the light leptonic channels ( = e, µ), the contribution from f 0 can be neglected, while in the τ channel it cannot.

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) [43,44], 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, the momentum transfer t to the light quark is insufficient to change its state, 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 heavy-quark spin symmetry re-

JHEP12(2018)114
In the decay B + c → η c , the spectator light quark is replaced by another heavy quark, c and some of these things will change. The presence of the heavy quark breaks some of the enhanced symmetries of the heavy-quark limit [45]. First, the difference between the heavyquark kinetic energy operators produces energies no longer negligible compared to those of the spectator c, spoiling 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 → η c 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 two form factors are related to a single, universal function h (∆ in ref. [45]), but only at the zero-recoil point, and no symmetry-based normalization for h can be derived [45].
Using the trace formalism of [46], in ref. [45] 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. [45]. This approximation was improved in ref. [47], where a constituent quark-model calculation of BR(B + c → η c + ν ) for = e, µ but not τ , was performed. The heavy-quark spin-symmetry relations were generalized in [29] to account for a momentum transfer to the spectator quark occurring at leading order in NRQCD. We reproduce here the relation of [29], where the form factors f + (w = 1) and f 0 (w = 1) are related by where r ≡ m/M , ρ ≡ m Q /m Q , and σ ≡ m q /m Q . These relations reproduce the standard Isgur-Wise result [43,44,48] when σ = 0. 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 → η c to those of B + c → J/ψ, which will be useful in the future to obtain further constraints on all six form factors.

Lattice QCD results
The state-of-the-art lattice QCD calculations for B + c → η c are limited to preliminary results from the HPQCD Collaboration for f + (q 2 ) at 4 q 2 values and f 0 (q 2 ) at 5 q 2 values [40,41]. f +,DA (q 2 ) f 0,DA (q 2 ) f +,lat (q 2 ) f 0,lat (q 2 ) Figure 2. B + c → η c form factors f + (q 2 ) (red circles) and f 0 (q 2 ) (blue triangles) from the HPQCD collaboration. 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.
These results were obtained using 2+1+1 HISQ ensembles, in which the smallest lattice spacing is a ≈ 0.09 fm (0.46 GeV −1 ), and the b quark is treated via NRQCD, are reproduced in figure 2. For q 2 = t − , 0 f 0 (q 2 ) has also been computed on coarser lattices and for lighter dynamical b-quark ensembles, which are used to check the accuracy and assess the uncertainty of the a ≈ 0.09 fm NRQCD results. In contrast to the situation for R(J/ψ), for R(η c ) both form factors have some lattice calculations, so the complications in treating unknown form factors is not required. Instead, the dispersive relations are sufficiently constraining that a rigorous error budget smaller than our naive 20% is the easiest way to reduce the error in R(η c ).

Dispersive relations
In this work we fit the form factors of B + c → η c using analyticity and unitarity constraints on two-point Green's functions and a conformal parameterization in the manner implemented by Boyd, Grinstein, and Lebed (BGL) [49] for the decays of heavy-light hadrons. This parameterization was extended to heavy-heavy systems in [26] with slightly different set of free parameters to simplify the computation, which we will utilize. Here we briefly sketch the necessary components.
The freedom to chose a value of q 2 allows us to compute χ(q 2 ) reliably in pQCD, far from where the two-point function receives nonperturbative contributions. 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 [54][55][56][57][58] used in ref. [48] 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. [24]. 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. [7,20,48]. For X we include only below-threshold B + c poles and a single two-body channel, B + c + η c , implying that our results provide very conservative bounds.
For B + c + η c , there are lighter two-body threshold with the correct quantum numbers that must be taken into consideration. The first physically prominent two-body production threshold in t occurs at B +D (see table 2). With this fact in mind, we define a new variable t bd ≡ (M B + M D ) 2 that corresponds to the first branch point in a given two-point function, while the B + c + η c 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

JHEP12(2018)114
Type Vector 1 − BD 6.337, 6.899, 7.012 Scalar 0 + BD 6.700, 7.108 Table 2. Lowest B + c states needed for Blaschke factors with t < t bc (whose relevant two-body threshold is indicated by "Lowest pair") for the J P channels of interest.
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 + [48,52], 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. 1 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 are from model calculations [61], 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 The analytic significance of Blaschke factors for heavy-hadron form factors was first noted in refs. [59,60].
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 → η c . The remaining factors are found in table 3. Transforming the dispersion-relation inequality, eq. (5.6), into z-space: which, upon dividing out the non-analytic terms, allows the expansion in z of 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:

JHEP12(2018)114
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 = λ κ . Evaluating at N = N opt , one finds From these expressions, we find that the semileptonic decays have z max,τ ≈ 0.022 and z max,µ ≈ 0.030, where each has a 1.3% variation, depending upon whether the BD or B * D threshold is the lowest branch point, t bd .

Results
Before presenting our prediction for R(η c ) we summarize the constraints the form factors f 0 and f + are required to satisfy: • The coefficients a n of each form factor are constrained by n a 2 n ≤ 1 from eq. (5.11), in particular, for the cases n = 1, 2, 3 investigated here. • Using eq. (3.3), the value of f + (t − ) is required to agree with f 0 (t − ), which is calculated from lattice QCD, within 50%.
Imposing these constraints, we perform our fit. Our third assumption relating the form factors through heavy quark spin symmetry is unimposed in [42], allowing us to reduce the uncertainty for f + (t − ). Gaussian-distributed points are sampled for the form factors f 0 and f + whose means are given by 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. This is a more conservative method than the χ 2 procedure [42]. 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 .   Table 5. Coefficients of f + and f 0 in the expansion from eq. (5.10) with n = 3 and f lat = 20%.
The resulting bands of allowed form factors are shown for f lat = 20% in figure 2, alongside the HPQCD results. One can see from table 4 that the overall error in our predictions is roughly equal to f lat , indicating that the underlying uncertainties in the lattice data are the dominant error in our results and that a 1% precise prediction is possible if the associated lattice data existed, at which point electroweak corrections become relevant. By varying the uncertainty in the heavy quark spin symmetry relations from 50% to 10%, we found a negligible contribution to the overall error, and therefore improving them are likely to have little effect on the results in the foreseeable future.
Having computed the form factors, we present predicted values for R(η c ) as a function of the truncation power n = 1, 2, 3 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, but as a conservative value, we suggest using the n = 3, f lat = 20% value of R(η c ) = 0.29 (5). In contrast to the case of R(J/ψ), we have more than three data points, and can therefore investigate the convergence more carefully. For f + , the series appears to rapidly converge such that neither a 2 nor a 3 can be distinguished from zero. The value we obtain of n a 2 f + ,n = 0.0016(2) could be used to slightly strengthen bounds in future dispersive analyses in the vector channel. For f 0 , the typical value of n a 2 n for n = 1 is O(10 −2 ), but for n = 2, 3 we find that a 2 2 ≈ 1 despite a 2 = 0.0(7) reflecting that while on average a 2 should be negligible, large fluctuations are permitted with the present uncertainties. Although the dispersive constraint is saturated in the n ≥ 2 case, the predictions for R(η c ) aren't observed to change outside of the uncertainties for increasing n. This confirms that while neglected higher-order terms can potentially have a 2 n ≈ 1, the suppression even for z max ≈ 0.03 is sufficient that the rapid convergence is still secured.
All model-dependent values for R(η c ) presented in table 1 comply with our result of R(η c ) = 0.29(5), albeit some, e.g. the anomalously large value of R(η c ) = 0.6(3) of [36], have seen their parameter space reduced. This general agreement gives us confidence in our result.

JHEP12(2018)114
The B + c → η c process has sufficient q 2 data, with the notable exception being f + (t − ), to compute R(η c ). Following [26], we reanalyze our dispersive fits with a synthetic data point f + (t − ) = 1 ± f lat to investigate it's potential constraining power. The resulting fits are found are indistinguishable from our current results within uncertainty. This implies that reducing the O(m c /m b , Λ QCD /m c ) theoretical uncertainty on the heavy quark spin symmetry constraints (eq. (3.3)) will have a negligible effect, since they are only applicable at t − . and increasing their accuracy is equivalent to a lattice QCD calculation of f + (t − ). Therefore, the best direction for improve would be obtained by future lattice results that can fully account for the systematics we have tried to estimate.

Discussion and conclusion
In this work we have presented a model-independent prediction of R(η c ) = 0.29 (5). While the near-term outlook for an experimental measurement of R(η c ) from LHCb measurement is poor, near-term lattice results promises to reduce the theoretical uncertainty to the 1% level which is sufficient to require consideration of electroweak corrections.
Even without improved lattice QCD calculations, potential areas of improvement are possible. Experience in the heavy-light sector and the fact that the R(J/ψ) bounds saturates the dispersive relations suggest that including multiple states that appear in the dispersion relation can provides complementary information to help constrain the form factors further, additionally one could include the lattice results for B → D ( * ) [22,23,[62][63][64][65] and Λ b → Λ c [66]. This would allow for a global, coupled set of predictions for the semileptonic ratios.