Form Factors and the Fragmentation Fraction Ratio
Abstract
We present a lattice quantum chromodynamics determination of the scalar and vector form factors for the decay over the full physical range of momentum transfer. In conjunction with future experimental data, our results will provide a new method to extract , which may elucidate the current tension between exclusive and inclusive determinations of this parameter. Combining the form factor results at nonzero recoil with recent HPQCD results for the form factors, we determine the ratios and . These results give the fragmentation fraction ratios and , respectively. The fragmentation fraction ratio is an important ingredient in experimental determinations of meson branching fractions at hadron colliders, in particular for the rare decay . In addition to the form factor results, we make the first prediction of the branching fraction ratio , where is an electron or muon. Current experimental measurements of the corresponding ratio for the semileptonic decays of mesons disagree with Standard Model expectations at the level of nearly four standard deviations. Future experimental measurements of may help understand this discrepancy.
HPQCD Collaboration
I Introduction
Studies of and meson decays at the Large Hadron Collider provide precision tests of the Standard Model of particle physics and are an important tool in the search for new physics. For example, the first observation of the rare decay , through a combined analysis by the LHCb and CMS collaborations Khachatryan et al. (2015); Aaij et al. (2017), tested the Standard Model prediction of the branching fraction. This decay is doublysuppressed in the Standard Model, but may have large contributions from physics beyond the Standard Model (see, for example, Ge et al. (2016)). Although the observed branching fraction is currently consistent with Standard Model expectations, there is still considerable room for new physics, given the experimental and theoretical uncertainties. Both LHCb and CMS are expected to reduce their errors significantly in Run II and tightening constraints on possible new physics requires a corresponding improvement in the theoretical determination of the Standard Model branching fraction.
Extraction of the meson branching fraction relies on the normalization channels and Adeva et al. (2009). The branching fraction can then be expressed as Khachatryan et al. (2015)
(1) 
where the are the fragmentation fractions, which give the probability that a quark hadronizes into a meson. The factors in this equation represent detector efficiencies and the factors denote the observed numbers of events.
The analysis of Khachatryan et al. (2015) used the value of , determined from LHCb experimental data Aaij et al. (2012, 2013); Storaci et al. (2013). The ratio depends on the kinematic range of the experiment, leading to the introduction of an additional systematic uncertainty in the value of to account for the extrapolation of the LHCb result to the CMS acceptance. Reducing sources of systematic uncertainties in the value of this ratio will improve the precision of the determination of the branching fraction. Indeed, an accurate value for the fragmentation fraction ratio is necessary for improved measurements of other meson decay branching fractions at the LHC Adeva et al. (2009).
The ratio of the fragmentation fractions, , can be expressed in terms of the ratios of form factors Fleischer et al. (2010, 2011),
(2) 
where is the scalar form factor of the semileptonic decay at . The first lattice calculations of the form factor ratios in Equation (2) using heavy clover bottom and charm quarks were published in Bailey et al. (2012). In addition, the form factors, and , for the semileptonic decay were determined with twisted mass fermions for the region near zero recoil in Atoui et al. (2014).
In this article we calculate the form factors, and , for the semileptonic decay . We present a determination of these form factors over the full physical range of momentum transfer, using the modified expansion for the chiralcontinuumkinematic extrapolation. We combine these form factor results with recent HPQCD results for the decay Na et al. (2015) to determine the ratios of and form factors relevant to the ratio of fragmentation fractions, .
We use the nonrelativistic (NRQCD) action for the bottom quarks and the Highly Improved Staggered Quark (HISQ) action for the charm quarks. Our form factors for have appeared already in Na et al. (2015). Here we first present form factor results and then proceed to the form factor ratios. We find
(3) 
This leads to
(4) 
and
(5) 
respectively. The uncertainties in these results are: the experimental statistical and systematic uncertainties; theoretical uncertainties (predominantly arising from a factor that captures deviations from naive factorization and, in Equation (5), an electroweak correction factor); and the uncertainties in our lattice input. In quoting these results, we have assumed that there are no correlations between the lattice results and the other sources of uncertainty.
In addition to determining the fragmentation fraction ratio relevant to the measurement of the branching fraction for the rare decay, , the semileptonic decay provides a new method to determine the CKM matrix element . There is a longstanding tension between determinations of from exclusive and inclusive measurements of the semileptonic meson decays (see, for example, Amhis et al. (2014, 2016a) and the review in Patrigiani et al. (2016)), although recent analyses suggest the tension has eased Gambino et al. (2016); Bigi and Gambino (2016). The decay has yet to be observed experimentally and consequently has received less theoretical attention than semileptonic decays of the meson. The studies that have been undertaken for the decay include calculations based on relativistic quark models Zhao et al. (2007); Bhol (2014), lightcone sum rules Li et al. (2009), perturbative factorization Fan et al. (2014) and estimates using the BetheSalpeter method Li et al. (2010); Chen et al. (2012). At present, there is one unquenched lattice calculation of the form factor at zero recoil Atoui et al. (2014). The FNAL/MILC collaboration has previously studied the ratio of the form factors of the and decays Bailey et al. (2012).
We determine the form factor for the semileptonic decay at zero momentum transfer to be and at zero recoil to be . Although experimental data is frequently presented in the form , the additional information provided by our calculation of the shape of the form factors throughout the kinematic range will, when combined with future experimental data, provide a new method to extract and may elucidate the puzzle of the tension between inclusive and exclusive determinations of this CKM matrix element.
In the next section we briefly outline the details of the calculation, including the gauge ensembles, bottomcharm currents and two and threepoint correlator construction. Our calculation closely parallels that presented in Na et al. (2015) for the semileptonic decay and we refer the reader to that work for further details. In Section III we discuss correlator fits to our lattice data and Section IV covers the chiralcontinuumkinematic extrapolations, which follows closely the methodology of Na et al. (2015). We explain how some of the correlations between the new data and the data are incorporated into the chiralcontinuumkinematic expansion. Section V presents our final results for the form factors, for and , and for and . We summarize in Section VI and in Appendix A we give the information necessary to reconstruct the form factors. The analogous details for form factors were summarized in Appendix A of Na et al. (2015).
Ii Ensembles, currents and correlators
Our determination of the form factors for the semileptonic decay closely parallels the analysis presented in Na et al. (2015). Here we simply sketch the key ingredients of the analysis and refer the reader to Sections II and III of Na et al. (2015) for more details of the lattice calculation.
We use five gauge ensembles, summarized in Table 1, generated by the MILC collaboration Bazavov et al. (2010). These ensembles include three “coarse” (with lattice spacing ) and two “fine” (with ) ensembles and incorporate flavors of AsqTad sea quarks.
Set  (sea)  

C1  2.647  0.005/0.050  2096  4  
C2  2.618  0.010/0.050  2256  2  
C3  2.644  0.020/0.050  1200  2  
F1  3.699  0.0062/0.031  1896  4  
F2  3.712  0.0124/0.031  1200  4 
In addition, we tabulate the light pseudoscalar masses on these ensembles, for both AsqTad and HISQ valence quarks, in Table 2. The difference in these masses captures discretization effects arising from partial quenching. We account for these effects in the chiralcontinuumkinematic expansion, which we discuss in more detail in Section IV.
Set  

C1  0.15971(20)  0.15990(20)  0.36530(29)  0.31217(20)  0.41111(12) 
C2  0.22447(17)  0.21110(20)  0.38331(24)  0.32851(48)  0.41445(17) 
C3  0.31125(16)  0.29310(20)  0.40984(21)  0.35720(22)  0.41180(23) 
F1  0.14789(18)  0.13460(10)  0.25318(19)  0.22855(17)  0.294109(93) 
F2  0.20635(18)  0.18730(10)  0.27217(21)  0.24596(14)  0.29315(12) 
In Table 3 we list the valence quark masses for the NRQCD bottom quarks and HISQ charm quarks Na et al. (2012); Bouchard et al. (2014). For completeness and ease of reference, we include both the treelevel wave function renormalization for the massive HISQ quarks Monahan et al. (2013) and the spinaveraged mass, corrected for electroweak effects, determined in Na et al. (2012).
Set  

C1  2.650  0.0489  0.6207  1.00495618  0.28356(15) 
C2  2.688  0.0492  0.6300  1.00524023  0.28323(18) 
C3  2.650  0.0491  0.6235  1.00504054  0.27897(20) 
F1  1.832  0.0337  0.4130  1.00103879  0.25653(14) 
F2  1.826  0.0336  0.4120  1.00102902  0.25558(28) 
To study semileptonic decays, we evaluate the matrix element of the bottomcharm vector current, , between and states. We express this matrix element in terms of the form factors and as
(6) 
where the momentum transfer is . In practice it is simpler to work with the form factors and , which are related to and via
(7)  
(8) 
Here is the energy of the daughter meson in the rest frame of the meson. In the following, we work in the rest frame of the meson and when we refer to the spatial momentum, , we mean the momentum of the meson.
NRQCD is an effective theory for heavy quarks and results determined using lattice NRQCD must be matched to full QCD to make contact with experimental data. We match the bottomcharm currents, , at one loop in perturbation theory through , where is the bare lattice mass Monahan et al. (2013). We rescale all currents by the nontrivial massive wave function renormalization for the HISQ charm quarks, tabulated in Table 3, Na et al. (2015).
We calculate and meson twopoint correlators and threepoint correlators of the bottomcharm currents, . We use smeared heavystrange bilinears to represent the meson and incorporate both deltafunction and Gaussian smearing, with a smearing radius of and on the coarse and fine ensembles, respectively. Threepoint correlators are computed with the setup illustrated in Figure 1. The meson is created at time and a current inserted at timeslice , between and . The daughter meson is then annihilated at timeslice . We use four values of : 12, 13, 14, and 15 on the coarse lattices; and 21, 22, 23, and 24 on the fine lattices. We implement spatial sums at the source through the random wall sources and Na et al. (2010). We generate data for four different values of the meson momenta, , , , and , where is the spatial lattice extent.
We fit meson twopoint functions to a sum of decaying exponentials in Euclidean time, ,
(9) 
Here the superscripts and indicate the smearing associated with the meson source (delta function or Gaussian); the and are amplitudes associated with the ordinary nonoscillatory states and the oscillatory states that arise in the staggered quark formalism; the meson energies are and for the nonoscillatory and oscillatory states, respectively; and is the number of exponentials included in the fit.
The ground state energy in NRQCD, , is related to the true energy in full QCD, , by
(10) 
because the quark rest mass has been integrated out in NRQCD. Here is the spinaveraged mass used to tune the quark mass and was determined in Na et al. (2012). We tabulate the values for in Table 3.
We fit the meson twopoint functions to the form
(11) 
For the threepoint correlator we use the fit ansatz
(12) 
The amplitudes for energy levels depend on the current , the daughter meson momentum , and the smearing of the meson source, .
The hadronic matrix element between and meson states is then given in terms of the ground state energies and amplitudes extracted from two and threepoint correlator fits by the relation
(13) 
For more details on this relation, see Section III of Na et al. (2015).
Iii Correlator fit and form factor results
We employ a Bayesian multiexponential fitting procedure, based on the
python
packages lsqfit
Lepage, G.P. (a) and corrfitter
Lepage, G.P. (b), that has been
used by the HPQCD collaboration for a wide range of lattice calculations.
Statistical correlations between data points, and correlations between data and
priors, are automatically captured with the gvar
class Lepage, G.P. (c),
which
facilitates the straightforward manipulation of Gaussiandistributed random
variables.
In this Bayesian multiexponential approach, one uses a number of indicators of fit stability, consistency, and goodnessoffit to check the fit results. For example, we check that, beyond a minimum number of exponentials, the fit results are independent of the number of exponentials included in the fit. Figure 2 illustrates the results of this test for the meson twopoint fits on ensemble set F1. The upper panel presents our results for four values of the spatial momentum, plotted as a function of the number of exponentials included in the plot. The lower panel shows the results obtained from three types of fits: a simultaneous fit to correlator data for all four spatial momenta, plotted with blue diamonds; a chained fit (discussed in detail in Appendix A of Bouchard et al. (2014)) to correlator data for all four spatial momenta simultaneously, shown with red squares; and an “individual” fit, plotted with purple circles. These individual fits include the correlator data for just a single daughter meson momentum in each fit.
We take the result for from the chained fit as our final result for each momentum. These results are tabulated in Table 4 and shown in Figure 2 as shaded bands in each plot. All three fit approaches give consistent results, as seen in the lower panel of Figure 2, but the simultaneous fits, with or without chaining, have the advantage that they capture the correlations between momenta, which is then reflected in the uncertainty quoted in the fit results. The chained fits give slightly better values of reduced . For example, for the ground state results plotted in the lower panel, the chained fits give for , while the simultaneous fits give . Both fits include 164 degrees of freedom. In addition, the chained fits are about ten percent faster than the simultaneous fits—14.6s to generate all the data in the lower plot for the chained fit compared to 16.4s for the simultaneous fit. This is not an important consideration for the twopoint fits, but becomes relevant for the larger threepoint fits, which can take many hours. Choosing to use chained fits for both two and threepoint fits ensures a consistent approach throughout the fitting procedure.
Set  

C1  1.18755(22)  1.21517(34)  1.24284(33)  1.27013(39) 
C2  1.20090(30)  1.24013(56)  1.27822(61)  1.31543(97) 
C3  1.19010(33)  1.23026(53)  1.26948(54)  1.30755(79) 
F1  0.84674(12)  0.87559(19)  0.90373(20)  0.93096(26) 
F2  0.84415(14)  0.87348(25)  0.90145(25)  0.92869(33) 
As a further test of the twopoint fits for the meson we determine the ratio on each ensemble. We plot the results in Figure 3. The shaded region corresponds to , where we set . In general, the data lie systematically above the relativistic value of unity, indicating that the statistical uncertainties of the fit results are sufficiently small that we can resolve discretization effects at . These discretization effects are less than in the dispersion relation.
Figure 4 shows the corresponding twopoint fit results for the ground state of the meson for ensemble sets C2 and F1. These ensemble sets have the same sea quark mass ratios, (see Table 1) and the difference between the results stems almost entirely from the lattice spacing. We take the values with as our final results, highlighted in the figure by the square data points and the shaded bands. We tabulate our final results in Table 5.
C1  C2  C3  F1  F2 
0.53714(60)  0.54332(65)  0.53657(86)  0.40873(53)  0.40819(44) 
For the threepoint correlator fits, we use a fitting procedure that diverges slightly from the approach taken in Na et al. (2015) and do not employ a “mixed” fitting strategy. Instead of combining “individual” and “master” fits (see Na et al. (2015) for full details), we use chained fits to correlators at all spatial momenta. This fitting approach ensures that we keep track of all statistical correlations between data at different momenta while maintaining fit stability, which was an issue for simultaneous fits attempted in Na et al. (2015).
To improve stability and goodnessoffit, we thin the threepoint correlator data on the fine ensembles by keeping every third timeslice. We illustrate the stability of these fits with the number of exponentials in the fit in Figure 5.
We test our choice by comparing fit results for the threepoint amplitudes with thinning (keeping both every third and every fifth timeslice) and without thinning and plot the results in Figure 6. We do not consider thinning by an even integer, which removes information about the oscillatory states generated by the staggered quark action.
In Figure 7 we present results for the threepoint fits when different combinations of sourcesink separations, , are used. For our final results we take the full set, on the coarse ensembles and on the fine ensembles.
We fit the threepoint correlator data after matching the bottomcharm currents to full QCD, as described briefly in Section II and in more detail in Na et al. (2015). In Na et al. (2015) this approach was compared with fitting the data first and then matching to full QCD and, as expected, the results are in good agreement within errors.
We summarize our final results for the form factors, and , for each ensemble and momentum in Tables 7 and 7. We represent the correlations between form factors at different momenta as a heat map in Figure 8 for ensemble set F2.
Set  

C1  0.8885(11)  0.8754(14)  0.8645(13)  0.8568(13) 
C2  0.8822(13)  0.8663(15)  0.8524(16)  0.8418(18) 
C3  0.8883(13)  0.8723(16)  0.8603(16)  0.8484(21) 
F1  0.90632(98)  0.8848(13)  0.8674(13)  0.8506(17) 
F2  0.9047(12)  0.8855(16)  0.8667(15)  0.8487(19) 
Set  

C1  1.1384(35)  1.1081(20)  1.0827(21) 
C2  1.1137(29)  1.0795(22)  1.0470(21) 
C3  1.1260(34)  1.0912(24)  1.0552(28) 
F1  1.1453(29)  1.0955(24)  1.0549(24) 
F2  1.1347(42)  1.0905(26)  1.0457(33) 
Iv Chiral, continuum and kinematic extrapolations
The form factor results presented in the previous section are determined at finite lattice spacing, with sea quark masses that are heavier than their physical values. These form factors are therefore functions of the momentum transfer, the lattice spacing, and the sea quark masses. The form factors determined from experimental data are functions of a single kinematic variable only. Typically this variable is the momentum transfer, , or the daughter meson energy, , but the form factors can also be expressed in terms of the variable, defined by
(14) 
where or the variable,
(15) 
Here and is a free parameter, which we take to be to ensure consistency with the analysis of Na et al. (2015). In Figure 9 we compare our results for the form factors, and , with the corresponding form factors for the decay, taken from Na et al. (2015), as a function of the variable. From the plot, we see that there is little dependence on the light spectator quark species in the form factor results.
To relate the form factor results determined at finite lattice spacing and unphysical sea quark masses to experimental data, we must therefore perform continuum and chiral extrapolations, along with a kinematic extrapolation in terms of one of the choices of kinematic variable. We combine these extrapolations through the modified expansion, introduced in Na et al. (2010, 2011), and applied to heavylight decays in Bouchard et al. (2013a, b); Bouchard et al. (2014). Our analysis of the chiralcontinuumkinematic extrapolation for decay closely parallels that for the decay in Na et al. (2015), so we only briefly outline the key components and refer the reader to Na et al. (2015) for details.
We express the dependence of the form factors on the variable through a modification of the BCL parameterization Bourrely et al. (2009)
(16)  
(17) 
Here the are Blaschke factors that take into account the effects of expected poles above the physical region,
(18) 
where we take Gregory et al. (2010), and . We find little dependence on the value of , in line with the results of Na et al. (2015). The expansion coefficients include lattice spacing and light quark mass dependence and can be written as
(19) 
where the include all lattice artifacts and chiral logarithms. These coefficients are given by
(20) 
where
(21)  
(22)  
(23) 
and the , , , and are fit parameters, along with the . We use the fit function form of Na et al. (2015), with a new fit parameter, , to account for the tuning of the valence strange quark mass on each ensemble. We tabulate the meson masses required to calculate in Table 2.
We further modify the expansion parameterization of the form factors to accommodate the systematic uncertainty associated with the truncation of the matching procedure at . We introduce fit parameters and , with central value zero and width and rescale the form factors, and according to
(24) 
We take the systematic uncertainties in these fit parameters as 3% and refer the reader to the detailed discussion of this approach in Na et al. (2015).
In Figure 10 we plot our fit results for , as a function of the variable. We obtain a reduced of with 36 degrees of freedom (dof), with a quality factor of . The value (or value) corresponds to the probability that the from the fit could have been larger, by chance, assuming the data are all Gaussian and consistent with each other. We plot the lattice data and the results of the chiralcontinuumkinematic extrapolation for as the upper, red shaded band and for as the lower, purple shaded band. We use the fit ansatz outlined above, including terms up to in the modified expansion, and refer to these results as the “standard extrapolation”. We tabulate our choice of priors and the fit results in Appendix A, and provide the corresponding expansion coefficients and their correlations in Table 11. Following Na et al. (2015) and the earlier work of Na et al. (2010, 2011), we group the priors into Group I and Group II variables, and add a third group. Broadly speaking, Group I priors are the typical fit parameters, Group II includes the input lattice scales and masses, and Group III priors are physical input masses. See the appendix of Na et al. (2015) for more details.
To test the convergence of our fit ansatz, we follow a procedure similar to that outlined in Na et al. (2015). This can be summarized as modifying the fit ansatz in the following ways:

include terms up to in the expansion;

include terms up to in the expansion;

add lightquark mass dependence to the fit parameters ;

add strangequark mass dependence to the fit parameters ;

add bottomquark mass dependence to the fit parameters ;

include discretization terms up to ;

include discretization terms up to ;

include discretization terms up to ;

include discretization terms up to ;

omit the term;

incorporate a 2% uncertainty for higherorder matching contributions;

incorporate a 4% uncertainty for higherorder matching contributions;

incorporate 4% and 2% uncertainties on coarse and fine ensembles, respectively, for higherorder matching contributions.
We show the results of these modifications in Figure 11. This plot demonstrates that the fit has converged with respect to a variety of modifications of the chiralcontinuumkinematic extrapolation ansatz. As part of this process, we also tested the significance of the Blaschke factor in the fit results. In line with the results of Na et al. (2015), we found that, while the results agreed within uncertainties, removing the Blaschke lowered the central value and increased the uncertainty of the result. This test is not strictly a test of convergence and is therefore not included in Figure 11.
To determine the ratio of form factors, we simultaneously fit the lattice form factor data for the and decays in a single script. We take the form factor results from Table III of Na et al. (2015) for the decay. Fitting the results simultaneously ensures that statistical correlations between the two data sets, such as those stemming from the lattice spacing determination on each ensemble set, are included in the final result for the ratio at zero momentum transfer. We do not reanalyze the to account for statistical correlations between the correlators themselves, which have negligible effect on the final result, given the current precision. This analysis would require fitting both and two and threepoint correlators simultaneously. To ensure that these statistical correlations are not important, we tested the correlations between the threepoint correlators on different ensemble sets. We show an example of the corresponding correlations as a heat map in Figure 12, from which one can see that statistical correlations are less than . We have found that correlations of this size have negligible impact at our current level of precision.
We fit the form factor data using the standard extrapolation ansätze for both the and data. For the decay, we choose the priors for the coefficients in the modified expansion to be equal to those for the corresponding expression for the expansion. These priors reflect the close agreement between the values for the and decays, illustrated in Figure 9. We list our choice of priors and the fit results for the ratio of form factors in Appendix A, and provide the corresponding expansion coefficients and their correlations in Table 12.
V Results
v.1 Form factors
We plot our final results for the form factors, and , as a function of the momentum transfer, , in Figure 13.
Our final result for the form factor at zero momentum transfer is
(25) 
We provide an estimate of the error budget for this result in Table 8. For the ratio of form factors, we find
(26) 
and
(27) 
with corresponding error budgets in Table 9. We show the extrapolation bands as a function of momentum transfer for both (purple hatched band) and (plain turquoise band) semileptonic decays in Figure 14.
We find agreement, within errors, with the results of Bailey et al. (2012), which are
(28)  
(29) 
Here we have combined the uncertainties quoted in Bailey et al. (2012), which are statistical and systematic, in quadrature.
For the form factor at zero recoil, , which is often quoted as
(30) 
where , we find
(31) 
This result is in good agreement with the value of determined in Atoui et al. (2014), with a slightly smaller uncertainty. The corresponding values for the form factors are Na et al. (2015) and Bailey et al. (2012) (where the quoted uncertainty includes only statistical uncertainties).