Form Factors for the full range from Lattice QCD with non-perturbatively normalized currents
We present a lattice QCD determination of the scalar and vector form factors over the full physical range of momentum transfer. The result is derived from correlation functions computed using the Highly Improved Staggered Quark (HISQ) formalism, on the second generation MILC gluon ensembles accounting for up, down, strange and charm contributions from the sea. We calculate correlation functions for three lattice spacing values and an array of unphysically light -quark masses, and extrapolate to the physical value. Using the HISQ formalism for all quarks means that the lattice current coupling to the can be renormalized non-perturbatively, giving a result free from perturbative matching errors for the first time. Our results are in agreement with, and more accurate than, previous determinations of these form factors. From the form factors we also determine the ratio of branching fractions that is sensitive to violation of lepton universality: , where is an electron or a muon. We find , which is also more accurate than previous lattice QCD results. Combined with a future measurement of , this could supply a new test of the Standard Model. We also compare the dependence on heavy quark mass of our form factors to expectations from Heavy Quark Effective Theory.
The weak decay processes of mesons such as the and , containing quarks, are a key potential source of insights into physics beyond the Standard Model (SM). Flavour-changing decays have gained a lot of interest because of a number of related tensions between experimental measurements and SM predictions Wei et al. (2009); Lees et al. (2012a, b, 2013); Aaij et al. (2014a, b); Huschle et al. (2015); Aaij et al. (2016a, 2015a, 2015b, b); Wehle et al. (2017); Sato et al. (2016); Hirose et al. (2018); Aaij et al. (2018a, b, 2017); Sirunyan et al. (2018); Aaboud et al. (2018). These tensions drive the need for improved theoretical calculations in the SM using methods and studying processes where we have good control of the uncertainties.
Lattice QCD is the method of choice for providing the hadronic input known as form factors that determine, up to a normalisation factor, the differential branching fraction for exclusive decay processes such (we suppress all electric charge and particle-antiparticle labels here in referring to decay processes). The normalisation factor that can then be extracted by comparison of theory with experiment is the Cabibbo-Kobayashi-Maskawa matrix element, in this case Aubert et al. (2009, 2010); Bailey et al. (2015); Na et al. (2015); Glattauer et al. (2016). Determination of then feeds into constraints on new physics through, for example, tests of the unitarity triangle.
There has been a long-standing tension in determinations of between exclusive (from and decays), and inclusive (from , where is any charmed hadronic state) processes. The most accurate exclusive results came from studies of the decay at zero recoil. It now seems likely that the uncertainties there were being underestimated because of the use of a very constrained parameterisation in the extrapolation of the experimental data to the zero recoil limit Bernlochner et al. (2017); Bigi et al. (2017a); Grinstein and Kobach (2017); Tanabashi et al. (2018), but see also Lees et al. (2019); Gambino et al. (2019). This underlines the importance in future of comparing theory and experiment across the full range of squared 4-momentum transfer () (a point emphasised for in Koponen et al. (2013)). It also demonstrates the need for comparison of accurate results from multiple decay processes for a more complete picture. Improved methods for producing the theoretical input to , namely lattice QCD determinations of form factors, are clearly necessary.
Here we provide improved accuracy for the form factors for the decay using a new lattice QCD method that covers the full range of the decay for the first time. Preliminary results appeared in McLean et al. (2018). The form factors are more attractive than for a first calculation to test methodology. They are numerically faster to compute and have higher statistical accuracy and smaller finite-volume effects because no valence quarks are present. Chiral perturbation theory Laiho and Van de Water (2006) expects that the form factors should be relatively insensitive to the spectator quark mass and hence should be very similar between and . This is confirmed at the 5% level by lattice QCD calculations Bailey et al. (2012); Monahan et al. (2017). Hence improved calculations of form factors can also offer information on .
Given an experimental determination, the decay can supply a new method for precisely determining the Cabibbo-Kobayashi-Maskawa (CKM) element . It can also supply a new test of the SM through quantities sensitive to lepton universality violation. We give the SM result for , where or . An experimental value for comparison to this would help to clarify the tension found between the SM and experiment in the related ratios Amhis et al. (2017) (see also a preliminary new analysis by Belle Caria (2019)).
Three lattice QCD calculations of form factors have already been performed. The FNAL/MILC collaboration Bailey et al. (2012, 2015) used the Fermilab action for the and quarks and the asqtad action for the light quarks on MILC gluon field ensembles that include flavours of asqtad sea quarks. On the same gluon field ensembles the HPQCD collaboration has calculated the form factors using Nonrelativistic QCD (NRQCD) for the valence and the Highly Improved Staggered Quark (HISQ) action for the other valence quarks Na et al. (2015); Monahan et al. (2017). A further calculation has been done using maximally twisted Wilson quarks on gluon field ensembles Atoui et al. (2014). Preliminary results using domain-wall quarks are given in Flynn et al. (2018).
A considerable limitation in the FNAL/MILC and HPQCD/NRQCD studies is the requirement for normalisation of the lattice QCD current. The matching between this current and that of continuum QCD is done in lattice QCD perturbation theory through , giving a systematic error which can be sizeable. Systematic errors coming from the truncation of the nonrelativistic expansion of the current are also a problem. In the Fermilab formalism the missing terms become discretisation effects on fine enough lattices; in the NRQCD formalism they mix discretisation effects and (where is the quark mass) relativistic corrections. Here we dispense with both of these problems by using a relativistic formalism with absolutely normalised lattice QCD currents.
Another limitation present in each of the previous studies is that the lattice QCD results are limited to a region of high , close to zero recoil. The reason for this is mainly to avoid large statistical errors. The signal/noise degrades exponentially as the spatial momentum of the meson in the final state grows. For quark decays the maximum spatial momentum of the final state meson can be large (tending to for light mesons, where is the meson mass). Systematic errors from missing discretisation (and relativistic) corrections also grow away from zero recoil. This is particularly problematic if discretisation effects are as above and relatively coarse lattices are used (to reduce numerical cost). Working close to zero recoil means that the lattice results then have to be extrapolated from the high region into the rest of the physical range. Here we also overcome this problem by working with a highly improved quark action in which even errors have been eliminated at tree-level Follana et al. (2007). We cover a range of values of the lattice spacing that includes very fine lattices and include results from lighter than physical quarks and this enables us to cover the full range in our lattice calculation.
We perform our calculation on the second-generation MILC gluon ensembles Bazavov et al. (2013), including effects from 2+1+1 flavours in the sea using the HISQ action Follana et al. (2007). We also use the HISQ action for all valence quarks. Our calculation employs HPQCD’s heavy-HISQ approach. In this we obtain lattice results at a number of unphysically light masses for the (we refer to this generically as the heavy quark ), reaching the quark mass on the finest lattices. This allows us to perform a combined fit in and lattice spacing that we can evaluate in the continuum limit at (and as a function of to compare, for example, to expectations from Heavy Quark Effective Theory (HQET)). By using only HISQ quarks, we can normalise all the lattice currents fully non-perturbatively and avoid systematic errors from current matching.
This calculation adds to a growing number of successful demonstrations of the heavy-HISQ approach. The method was developed for determination of the quark mass, meson masses and decay constants McNeile et al. (2010, 2012a, 2012b) and is now also being used by other groups for these calculations Bazavov et al. (2018); Petreczky and Weber (2019). A proof-of-principle application of heavy-HISQ to form factors was given for and in Lytle et al. (2016); Colquhoun et al. (2016), covering the full range for these decays and this work builds on those results. The axial form factor at zero recoil was calculated using heavy-HISQ in McLean et al. (2019).
This article is structured in the following way: Section II lays out our lattice QCD approach to calculating the form factors and then Section III presents our results, along with several consistency checks and our determination of . We also give curves showing the heavy-quark mass dependence of some features of the form factors that can be compared to HQET. For those simply hoping to use our calculated form factors, Appendix A gives the parameters and covariance matrix required to reconstruct them.
Ii Calculation Details
ii.1 Form Factors
In this section we specify our notation for the form factors and matrix elements. The differential decay rate for decays is given in the SM by
where is the mass of the lepton, is the electroweak correction, is the momentum transfer and , are the scalar and vector form factors that parameterize the fact that the decay process involves hadrons. We use superscript ‘’ to denote the strange spectator valence quark. The allowed range of values if the final states are on-shell is
The form factors are determined from matrix elements of the electroweak current between and states, where is the vector component and is the axial vector component. In a pseudoscalar-to-pseudoscalar amplitude, only contributes, since does not satisfy the parity invariance of QCD. In terms of form factors, the vector current matrix element is given by
Analyticity of this matrix element demands that
Via the partially conserved vector current relation (PCVC), the form factor is also directly related to the matrix element of the scalar current ;
In our calculation we determine the form factors by computing matrix elements of the temporal vector current and the scalar current . The form factors can be extracted from this combination using expressions derived from Eqs. (3) and (5) (once the currents have the correct continuum normalisation - see Section II.4):
Our goal is to compute and throughout the range of values . We extend the range to in order to take advantage of the constraint from Equation (4).
|1||0.0376||0.45||0.5||0, 0.056||14, 17, 20|
|0.65||0, 0.142, 0.201|
|0.8||0, 0.227, 0.323|
|0.036||0.433||0.5||0, 0.0279||14, 17, 20|
|0.0234||0.274||0.427||0, 0.113, 0.161||22, 25, 28|
|0.525||0, 0.161, 0.244|
|0.65||0, 0.244, 0.338|
|0.8||0, 0.338, 0.438|
|0.0165||0.194||0.5||0, 0.202, 0.281||31, 36, 41|
|0.65||0, 0.202, 0.281, 0.382|
|0.8||0, 0.281, 0.382, 0.473|
ii.2 Lattice Calculation
This calculation closely follows the approach employed in our calculation of the axial form factor at zero recoil McLean et al. (2019). Here, however, we must give spatial momentum to the charm quark in the final state so that we can cover the full range of the decay.
The gluon field configurations used in this calculation were generated by the MILC collaboration Bazavov et al. (2010, 2013). The relevant parameters for the specific ensembles we use are given in Table 1. The gluon fields are generated using a Symanzik-improved gluon action with coefficients matched to continuum QCD through Hart et al. (2009). The gluon fields include the effect of 2+1+1 flavours of quarks in the sea (, , , , where ) using the HISQ action Follana et al. (2007). In three of the four ensembles (sets 1, 3 and 4), the bare light quark mass is set to . The fact that the value is unphysically high is expected to have only a small effect on the form factors here, since we have no valence light quarks. We quantify this small effect by including a fourth ensemble (set 2) with roughly physical .
We use a number of different masses for the valence heavy quark . This allows us to resolve the dependence of the form factors on the heavy quark mass, so that a fit in can be performed and the results of the fit evaluated at . With a heavy quark mass varying both on a given ensemble and between ensembles, we can resolve both the discretisation effects that grow with large () masses and the physical dependence of the continuum form factors on . Using unphysically light -quarks also reduces the range, meaning that we can obtain lattice results across the full range while the statistical noise remains under control.
Staggered quarks have no spin degrees of freedom. Spin-parity quantum numbers are accounted for by construction of appropriate fermion bilinears and including an appropriate space-time dependent phase with each operator in the path integral. We categorize these phases according to the standard spin-taste notation, , where is the spin structure of the operator in the continuum limit, and is the ‘taste’ structure which accounts for the multiple possible copies of the operator constructed from staggered quark fields.
We have designed this calculation to use only local operators (combining fields at the same space-time point and having spin-taste) for the calculation of the current matrix elements that we require. This is an advantage since point-split operators can lead to noisier correlation functions. The spin-taste operators we use are: scalar , pseudoscalar , vector , and temporal axial-vector .
We compute a number of correlation functions on the ensembles detailed in Table 1. Valence quark masses, momenta and other inputs to the calculation are given in Table 2. We use random wall sources to generate all staggered propagators from the source since this gives improved statistical errors Davies et al. (2010a). First we compute two-point correlation functions between meson eigenstates of momentum ,
where represents a functional integral over all fields, are valence quark fields of the flavours the meson is charged under, is the spin-taste structure of and the division by the number of tastes is required to normalise closed loops made from staggered quarks Follana et al. (2007). We compute these correlation functions for all values, i.e. .
We compute correlation functions for a heavy-strange pseudoscalar, , with spin-taste structure , at rest. In terms of staggered quark propagators this takes the form
where is a staggered propagator for flavour , and the trace is over color. Here and , and the sum is over spatial sites labelled , . We also compute correlators for a charm-strange pseudoscalar meson , with structure . For these correlators we need both zero and non-zero spatial momentum. Non-zero spatial momentum is given to the by imposing twisted boundary conditions on the gluon fields when computing the charm quark propagators Guadagnoli et al. (2006). Then
where denotes a propagator with momentum twist . We compute these correlation functions using several different twists to produce the range of momenta given in Table 2. We design the propagators to have momentum , by imposing a twist in each spatial direction.
Necessary for extracting the vector current matrix element, we also compute correlation functions for a non-goldstone pseudoscalar heavy-strange mesons at rest, denoted . This has spin-taste structure . correlators are computed using
where we use the notation .
We also compute correlators for mesons, heavy-charmed pseudoscalars, using the same form as those for , Equation (8). These are used to find decay constants, which are useful in some of our continuum and fits. In our fits to heavy-quark mass dependence we will use the mass of the heavy-heavy pseudoscalar meson, as a physical proxy for the quark mass. To quantify mistuning of the charm and strange quark masses, we also require masses for and mesons, identical to with replaced and quarks respectively. We compute correlators for each of these at rest, using a spin-taste structure , taking the same form as those of the , Equation (8). Note that all of the mesons discussed here are artificially forbidden to annihilate in our lattice QCD calculation. We expect this to have negligible effect, for the purposes of this calculation, on the masses of the and the Chakraborty et al. (2015); the is an unphysical meson that can be defined in this limit in a lattice QCD calculation and is convenient for tuning the quark mass Davies et al. (2010b); Dowdall et al. (2013).
Three-point correlation functions are needed to allow determination of the current matrix elements for decay. We require two sets of such correlation functions, one with a scalar and one with a temporal vector current insertion. The first takes the form
In terms of the staggered quark formalism, both the source and sink are given structure , and the current insertion . We combine staggered propagators to construct these correlation functions as:
where we fix , and , and once again the charm propagator is given the appropriate twist . We compute these correlation functions for all values within , using 3 values to make sure that excited state effects are accounted for. The values vary with lattice spacing to give approximately the same physical range and always include both even and odd values. The values are given in Table 2.
The three-point correlation function with temporal vector current insertion is given by
This is generated using spin-taste at the source, at the sink, and at the current insertion. To achieve this we compute
The non-goldstone is required here to ensure that taste cancels in the correlation function. The difference between the and the , for example in their masses, is generated by taste-exchange discretisation effects. In practice it is very small for heavy mesons Follana et al. (2007), being suppressed by the heavy meson mass.
ii.3 Analysis of Correlation Functions
We now describe our simultaneous multi-exponential fits to the correlation functions using a standard Bayesian apporach Lepage et al. (2002); cor (2018). The parameters that we wish to determine are ground-state energies, two-point amplitudes and ground-state to ground-state matrix elements. Our correlation functions, however, are contaminated by contributions from excited states. These excited states must be included in our fits so that the systematic error on the ground-state parameters from the presence of the excited states is fully taken into account. Multi-exponential fits are then mandatory, guided by Bayesian priors for the parameters, discussed below. To reduce the number of exponentials needed by the fits, we drop values of the correlation functions when they are within of the end-points (where excited states contribute most). We use values of varying from 2 to 10 throughout the correlator fits. We take results from fits using 5 exponentials ( in the fit forms below), where good values are obtained and the ground-state parameters and their uncertainties have stabilised.
Two-point correlation functions are fit to the form
and , are fit parameters. The second term in Equation (II.3) accounts for the opposite-parity states that arise from the staggered quark time doublers and are known as oscillating states (see Appendix G of Follana et al. (2007)). These oscillating states do not appear when is a Goldstone-taste pseudoscalar with a quark and antiquark of the same mass, so in the and cases the second term is not required.
Figure 1 shows the quality of our results. We plot effective energies and amplitudes for the and mesons on the fine ensemble. The effective energy is defined as
Here is a time-smeared mean correlator to reduce the impact of the oscillating states:
This again converges to the ground-state amplitude from the fit.
For the three-point correlation functions we use the fit form
This includes fit parameters common to the fits of (when ), (when ) and two-point correlators, along with new fit parameters that are related to the current matrix elements. We perform a single simultaneous fit containing each correlator computed () at every and every , for each ensemble.
These simultaneous fits are very large and this causes problems for the covariance matrix which must be inverted to determine . We take two steps toward mitigating this. The first is to impose an svd (singular value decomposition) cut . This replaces any eigenvalue of the covariance matrix smaller than with , where is the largest eigenvalue of the matrix. The small eigenvalues are driven to zero if the statistics available are not high enough. The application of the svd cut makes the matrix less singular, and can be considered a conservative move since the only possible effect on the error of the final results is to inflate them. An appropriate value for is found by comparing estimates of covariance matrix eigenvalues between different bootstrap samples of the data using the Corrfitter package cor (2018). The resulting varies between ensembles since it depends on the statistical quality of the dataset, but we find them to be of order .
The other step we take towards a stable fit is employing a chained-fitting approach. We first perform an array of smaller fits, each fitting the correlators relevant only to one and one value. In the case of set 4, for example, this results in 11 separate fits. Then, a full simultaneous fit of all of the correlators is carried out, using as priors the results of the smaller fits. This both speeds up the full fit and improves stability of the results.
The priors for the fits were set up as follows. We set gaussian priors for the parameters , and log-normal priors for amplitudes , ground-state energies , and excited-state energy differences . Using log-normal distributions forbids ground-state energies, excited state energy differences and amplitudes moving too close to zero or becoming negative, improving stability of the fit.
Priors for ground state energies and amplitudes are set according to an empirical-Bayes approach, plots of the effective amplitude of the correlation functions are inspected to deduce reasonable priors. The ground-state oscillating parameters , , are given the same priors as the non-oscillating states, with uncertainties inflated by 50%. The resulting priors always have a standard deviation at least 10 times that of the final result. The logs of the excited-state energy differences are given prior values where was taken as 0.5 GeV. The log of oscillating and non-oscillating excited state amplitudes are given priors of . The ground-state non-oscillating to non-oscillating three-point parameter, is given a prior of , and the rest of the three-point parameters are given .
The physical quantities that we need here are extracted from the ground-state fit parameters and given in Tables 3 and 4. are the ground-state meson energies in lattice units. For mesons at rest, this corresponds to the mass of the meson, i.e. . The annihilation amplitude for an -meson at rest is given in lattice units by
If is a pseudoscalar operator , the decay constant can be found from this via
where , are the quark flavours that is charged under. We use this to determine the meson decay constant in Table 3. The current matrix elements that we are focussed on here can be extracted from the fit parameters via
These can be converted into values for the form factors once the currents have been normalised (Section II.4).
Figure 2 shows the results of a number of tests we performed on the fits to correlators on the fine ensemble. Each test modifies one of the features of the fits and we then plot the resultant value of the key output parameter . The robustness of the fits can be gauged by the effect of these changes, which are all small.
Figure 3 shows a comparison of the meson dispersion relation on the fine (set 1) and superfine (set 3) lattices. The dispersion relation is sensitive to discretisation effects in our quark action. The figure shows them to be small (see Donald et al. (2012) for more discussion of discretisation effects in dispersion relations for mesons using HISQ quarks).
ii.4 Current Normalization
In the HISQ formalism, the local scalar current (multiplied by the mass difference of flavours it is charged under) is conserved, and hence requires no renormalization. This is not the case for the local temporal vector current . We use this instead of the conserved vector current because it is much simpler, but we then require a renormalisation factor to match to the continuum current. This is simple to obtain fully non-perturbatively within this calculation Koponen et al. (2013); Donald et al. (2014), at no additional cost.
When both meson states in the matrix elements are at rest (the zero recoil point), the scalar and local vector matrix elements are related via the PCVC relation:
can be extracted from this relation using the matrix elements we have computed. The values found on each ensemble and for each are given in Table 5.
We also remove tree-level mass-dependent discretisation effects from the current using a normalization constant, derived in Monahan et al. (2013) and discussed in detail in McLean et al. (2019). values are also tabulated in Table 5; they have only a very small effect.
Combining these normalizations with the lattice current from the simultaneous correlation function fits, we find values for the form factors at a given heavy mass, lattice spacing, and :
where is defined in Equation (II.1) and we have made the dependence of the matrix elements on explicit.
ii.5 Obtaining a result at the Physical Point
We now discuss how we fit our results for and as a function of valence heavy quark mass, sea light quark mass and lattice spacing. Evaluating these fits at the mass of the , with physical and masses and zero lattice spacing will then give us the physical form factor curves from which to determine the differential decay rate, using Equation (II.1).
Following McLean et al. (2019) we use two methods; one a direct approach to fitting the form factors and the other in which we fit the ratio
in which discretisation effects are somewhat reduced. We will take our final result from the direct approach and we describe that here. We use the ratio approach as a test of uncertainties and we describe that in more detail in Appendix B.
We use identical fit functions for both approaches. We feed into the fit our results from Tables 3 and 4, retaining the correlations (not shown in the Tables) between values for different heavy quark masses and values on a given gluon field ensemble that we are able to capture in our simultaneous fits (Section II.3). We also include, where needed, correlated lattice spacing uncertainties.
ii.5.1 Kinematic Behaviour
Our fit form is a modified version of the Bourrely-Caprini-Lellouch (BCL) parameterisation for pseudoscalar-to-pseudoscalar form factors Bourrely et al. (2009):
The functon maps to a small region inside the unit circle on the complex plane, defined by
Here and we choose to be . This choice means that maps to and the fit functions simplify to . For the physical range of for to decay, the range covered by is , resulting in a rapidly converging series in powers of . We truncate at ; adding further powers of does not effect the results of the fit.
The factors in front of the sums in the BCL parameterisation account for lowest mass pole expected in the full plane for each form factor coming from the production of on-shell and states in the crossed channel of the semileptonic decay. Note that these poles, even though they are below the cut (for production) that begins at , are at much higher values than those covered by the semileptonic decay here (with maximum given by ).
We must estimate , the scalar heavy-charm meson mass, at each of the heavy masses we use. For this we use the fact that the splitting is an orbital excitation and therefore largely independent of the heavy quark mass. The splitting has been calculated in Dowdall et al. (2012) to be GeV at the quark mass. Combined with an mass from our lattice results, we construct the mass as . We do not include the uncertainty on in the fit, since any shift in the precise position of the pole will be absorbed into the other fit parameters.
To estimate , the vector heavy-charm meson mass, we use the fact that the hyperfine splitting should vanish in the infinite limit. then takes the approximate form . To reproduce this behaviour we use the ansatz , and fix at the quark mass using the value of from Dowdall et al. (2012). This gives .
ii.5.2 Heavy Quark Mass and Discretisation Effects
To account for dependence on the heavy quark mass and discretisation effects in a general way, we use the following form for each of the coefficients:
To understand this form, focus first on the terms inside the sum. Powers of allow for variation of the coefficients as the heavy quark mass changes, using an HQET-inspired form since this is a heavy-light to heavy-light meson transition. is proportional to at leading order in HQET, so is a suitable physical proxy for the heavy quark mass. We take here to be 0.5GeV. The other two terms in the sum allow for discretisation effects. These can be set by two scales. One is the variable heavy quark mass