B s →Kℓν form factors from lattice QCD
C.M. Bouchard, G. Peter Lepage, Christopher Monahan, Heechang Na, Junko Shigemitsu
BB s → K(cid:96)ν form factors from lattice QCD
C.M. Bouchard, ∗ G. Peter Lepage, Christopher Monahan, Heechang Na, and Junko Shigemitsu (HPQCD Collaboration) Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA Laboratory of Elementary Particle Physics, Cornell University, Ithaca, New York 14853, USA Physics Department, College of William and Mary, Williamsburg, Virginia 23187, USA Department Physics and Astronomy, University of Utah, Salt Lake City, Utah 84112, USA (Dated: July 3, 2018)We report the first lattice QCD calculation of the form factors for the standard model tree-leveldecay B s → K(cid:96)ν . In combination with future measurement, this calculation will provide an alter-native exclusive semileptonic determination of | V ub | . We compare our results with previous modelcalculations, make predictions for differential decay rates and branching fractions, and predict theratio of differential branching fractions between B s → Kτ ν and B s → Kµν . We also present stan-dard model predictions for differential decay rate forward-backward asymmetries and polarizationfractions and calculate potentially useful ratios of B s → K form factors with those of the fictitious B s → η s decay. Our lattice simulations utilize nonrelativistic QCD b and highly improved staggeredlight quarks on a subset of the MILC Collaboration’s 2 + 1 asqtad gauge configurations, includingtwo lattice spacings and a range of light quark masses. PACS numbers: 12.38.Gc, 13.20.He, 14.40.Nd, 14.40.Df
I. INTRODUCTION
The decay B s → K(cid:96)ν occurs at tree-level in thestandard model via the flavor-changing charged-current b → u transition, making it an alternative to B → π(cid:96)ν in the determination of | V ub | from exclusive semileptonicdecays. The difference in these processes, a spectatorstrange quark in B s → K(cid:96)ν vs a spectator down quarkin B → π(cid:96)ν , is beneficial for lattice QCD simulations,because it improves the ratio of signal to noise. Thoughthis process has not yet been observed, its measurementis planned at LHCb and is possible during an Υ(5 S ) runat BelleII. This provides a prediction opportunity for lat-tice QCD.In addition to the calculation of form factors for B s → K , we also calculate their ratios with form fac-tors for the fictitious B s → η s decay. Such ratios areessentially free of our largest systematic error, perturba-tive matching. In combination with a future calculationof B s → η s using a highly improved staggered (HISQ) b quark, these ratios would yield a non-perturbative eval-uation of the matching factor for the b → u current withnonrelativistic QCD (NRQCD) b quark. This matchingfactor would be applicable to B s → K(cid:96)ν and B → π(cid:96)ν simulations using NRQCD b quarks.To include correlations among the data for both de-cays, correlation function fits must include vast amountsof correlated data. To make such fits feasible, we havedeveloped a new technique, called chaining, discussed inAppendix A. In addition, the use of marginalization tech-niques developed in Ref. [1] significantly reduces the timerequired for the fits. ∗ [email protected] The chiral, continuum, and kinematic extrapolationsare performed simultaneously using the modified z ex-pansion [2, 3] with the chiral logarithmic corrections fixedby the results of hard pion chiral perturbation theory(HPChPT) [4, 5]. The factorization of chiral correctionsand kinematics, as found at one-loop order by HPChPT,suggests the modified z expansion is a natural choice forcarrying out this simultaneous extrapolation. We referto the combination of HPChPT chiral logarithmic cor-rections and the modified z expansion as the HPChPT z expansion. II. FORM FACTORS AND MATRIX ELEMENTS
The vector hadronic matrix element is parametrizedby the scalar and vector form factors f , + (cid:104) K | V µ | B s (cid:105) = f + (cid:18) p µB s + p µK − M B s − M K q q µ (cid:19) + f M B s − M K q q µ , (1)where V µ = ¯ uγ µ b and q µ = p µB s − p µK . At intermediatestages of the calculation we recast f , + in terms of themore convenient form factors f (cid:107) , ⊥ (cid:104) K | V µ | B s (cid:105) = (cid:112) M B s (cid:18) p µB s M B s f (cid:107) + p µ ⊥ f ⊥ (cid:19) , (2)where p µ ⊥ = p µK − p µB s ( p K · p B s ) /M B s . In the B s mesonrest frame, the form factors f (cid:107) , ⊥ are simply related to thetemporal and spatial components of the hadronic vectormatrix elements, (cid:104) K | V | B s (cid:105) = (cid:112) M B s f (cid:107) , (3) (cid:104) K | V k | B s (cid:105) = (cid:112) M B s p kK f ⊥ . (4) a r X i v : . [ h e p - l a t ] O c t The scalar and vector form factors are related to f (cid:107) , ⊥ by f = (cid:112) M B s M B s − M K (cid:2) ( M B s − E K ) f (cid:107) + p K f ⊥ (cid:3) , (5) f + = 1 (cid:112) M B s (cid:2) f (cid:107) + ( M B s − E K ) f ⊥ (cid:3) , (6)where p K is the kaon three-momentum. This discussiongeneralizes in a straightforward way for the B s → η s matrix elements. III. SIMULATION
Ensemble averages are performed with the MILC Col-laboration’s 2 + 1 asqtad gauge configurations [6] listedin Table I. Valence quarks in our simulation are nonrel-ativistic QCD (NRQCD) [7] b quarks, tuned in Ref. [8],and highly improved staggered (HISQ) [9] light and s quarks, the propagators for which were generated inRefs. [2, 3]. Valence quark masses for each ensemble usedin the simulations are collected in Table I and correspondto pion masses ranging from, approximately, 260 MeV to500 MeV.Heavy-light B s meson bilinears Φ αB s are built fromNRQCD b and HISQ s quarks (for details see Ref. [8])and light-light kaon (and similarly for the η s ) bilinearsΦ K are built from HISQ light and s quarks (for detailssee Ref. [2]). From these bilinears we build two and threepoint correlation function data C αβB s ( t , t ) = 1 L (cid:88) x , y (cid:104) Φ βB s ( t, y ) Φ α † B s ( t , x ) (cid:105) , (7) C K, p ( t , t ) = 1 L (cid:88) x , y e i p · ( x − y ) (cid:104) Φ K ( t, y ) Φ † K ( t , x ) (cid:105) , (8) C αJ, p ( t , t, T ) = 1 L (cid:88) x , y , z e i p · ( z − x ) ×(cid:104) Φ K ( t + T, x ) J ( t, z ) Φ α † B s ( t , y ) (cid:105) , (9)where indices α, β specify b quark smearing. We gener-ate data for both a local and Gaussian smeared b quark,with smearing function φ introduced via the replacement (cid:80) y → (cid:80) y , y (cid:48) φ ( y (cid:48) − y ) in Eqs. (7) and (9). Three pointand daughter meson two point correlation function dataare generated at four daughter meson momenta, cor-responding to p L ∈ π { (000) , (100) , (110) , (111) } . Inthree point data, these momenta are inserted at x inFig. 1. The sum over x in Eqs. (8) and (9) is per-formed using random wall sources with U(1) phases ξ ,i.e. (cid:80) x → (cid:80) x , x (cid:48) ξ ( x ) ξ ( x (cid:48) ). In the three point correla-tor a B s meson source is inserted at timeslice t , selectedat random on each configuration to reduce autocorrela-tions. The current J is inserted at timeslices t such that t ≤ t ≤ t + T and the daughter meson is annihilated attimeslice t + T . Prior to performing the fits, all data are φ ( y ′ − y ) b ξ ( x ′ ) , ξ ( x ) ut, z t , y t , y ′ t + T, x ′ t + T, x Js FIG. 1: Setup for three point correlator data generation.shifted to a common t = 0. This three point correlatorsetup is depicted in Fig. 1. Additional details regardingthe two and three point correlation function generationcan be found in Ref. [10].The flavor-changing current J is an effec-tive lattice vector current V µ corrected through O ( α s , Λ QCD /m b , α s / ( am b )). The lattice currents thatcontribute through this order are V (0) µ = Ψ u γ µ Ψ b , (10) V (1) µ = − am b Ψ u γ µ γ · ∇ Ψ b . (11)Matrix elements of the continuum vector current (cid:104) V µ (cid:105) arematched to those of the lattice vector current accordingto (cid:104) V µ (cid:105) = (1 + α s ρ ( V µ )0 ) (cid:104)V (0) µ (cid:105) + (cid:104)V (1) , sub µ (cid:105) , (12)where (cid:104)V (1) , sub µ (cid:105) ≡ (cid:104)V (1) µ (cid:105) − α s ζ V µ (cid:104)V (0) µ (cid:105) . (13)The matching calculation is done to one loop using mass-less HISQ lattice perturbation theory [11]. In implement-ing the matching, we omit O ( α s Λ QCD /m b ) contribu-tions. Ref. [12], which used asqtad valence quarks, foundcontributions of this order to be negligible. In Ref. [10],which used HISQ valence quarks, these contributions tothe temporal component of the vector current were stud-ied and again were found to be negligible. We also omit O (Λ QCD /m b ) relativistic matching corrections. These,and higher order, omitted contributions to the matchingresult in our leading systematic error. An estimate of thiserror, and its incorporation in our fit results, is discussedin the following section. IV. CORRELATION FUNCTION FITS
Two and three point correlation function fit Ans¨atze,and the selection of priors, closely follows the methods ofTABLE I: Left to right: labels for the ensembles used in this analysis; lattice volume; inverse lattice spacing in r -units; light/strange sea-quark masses; tadpole improvement factor u = (cid:104) plaquette (cid:105) / ; number of configurations;number of time sources; valence u -quark mass; valence s -quark mass; b -quark mass; and the spin-averaged b ¯ b groundstate energies used to relate our B s meson simulation energies to their physical values. Ensemble L × N t r /a au m sea u N conf N tsrc am u am s am b aE sim b ¯ b C1 24 ×
64 2.647(3) 0.005/0.05 0.8678 1200 2 0.0070 0.0489 2.650 0.28356(15)C2 20 ×
64 2.618(3) 0.01/0.05 0.8677 1200 2 0.0123 0.0492 2.688 0.28323(18)C3 20 ×
64 2.644(3) 0.02/0.05 0.8688 600 2 0.0246 0.0491 2.650 0.27897(20)F1 28 ×
96 3.699(3) 0.0062/0.031 0.8782 1200 4 0.00674 0.0337 1.832 0.25653(14)F2 28 ×
96 3.712(4) 0.0124/0.031 0.8788 600 4 0.01350 0.0336 1.826 0.25558(28)
Ref. [10]. Two point B s data are fit to C αβB s ( t ) = N − (cid:88) n =0 b α ( n ) b β ( n ) † e − E sim( n ) Bs t + ˜ N − (cid:88) m =0 ˜ b α ( m ) ˜ b β ( m ) † ( − t e − ˜ E sim( m ) Bs t , (14)where tildes denote oscillating state contributions and E sim B s is the simulated B s energy. The physical groundstate B s mass is related to the simulation ground stateenergy by E (0) B s = E sim(0) B s + 12 ( M expt b ¯ b − E sim b ¯ b ) (15)where M expt b ¯ b = 9 . η b annihilation, andcharmed sea effects not present in our simulations, and E sim b ¯ b is the spin-averaged energy of b ¯ b states calculatedon the ensembles used in the simulation and listed in Ta-ble I. The b quark smearing is indicated by indices α, β .Kaon and η s two point correlator data are fit to an ex-pression of the form C p ( t ) = N − (cid:88) n =0 | d ( n ) p | (cid:0) e − E ( n ) t + e − E ( n ) ( N t − t ) (cid:1) + ˜ N − (cid:88) m =0 | ˜ d ( m ) p | ( − t (cid:0) e − ˜ E ( m ) t + e − ˜ E ( m ) ( N t − t ) (cid:1) . (16)Results of two point fits satisfy the dispersion relationand are stable with respect to variations in ( N, ˜ N ) andthe range of timeslices included in the fits, as demon-strated for kaon two point data in Ref. [10]. The zero momentum η s has no oscillating state contributions dueto mass degeneracy of its valence quarks. Three point correlation function data are described by C αJ, p ( t, T )= N − (cid:88) n,m =0 d ( n ) p A ( n,m ) J, p b α ( m ) † e − E ( n ) ( T − t ) e − E sim( m ) Bs t + N − (cid:88) n =0 ˜ N − (cid:88) m =0 d ( n ) p B ( n,m ) J, p ˜ b α ( m ) † ( − t e − E ( n ) ( T − t ) e − ˜ E sim( m ) Bs t + ˜ N − (cid:88) n =0 N − (cid:88) m =0 ˜ d ( n ) p C ( n,m ) J, p b α ( m ) † ( − T − t e − ˜ E ( n ) ( T − t ) e − E sim( m ) Bs t + ˜ N − (cid:88) n,m =0 ˜ d ( n ) p D ( n,m ) J, p ˜ b α ( m ) † ( − T e − ˜ E ( n ) ( T − t ) e − ˜ E sim( m ) Bs t , (17)where the three point amplitudes A , B , C , and D are pro-portional to the hadronic matrix elements. The groundstate hadronic matrix element is obtained from A (0 , √ A (0 , J, p = a (cid:104) K (0) p | J | B (0) s (cid:105) (cid:113) a E (0) K (cid:113) a E (0) B s , (18)where the factor of 4 / √ T betweenthe mother and daughter mesons. On the coarse ensem-bles we include data for T = 13 , ,
15 while for the fineensembles we include T = 23 ,
24 data.On each ensemble we perform a simultaneous fit to twoand three point correlation function data for the B s → K and B s → η s decays, at all simulated momenta, includ-ing both spatial and temporal currents, and for the tem-poral separations listed above. This ensures correlationsamong these data are accounted for in the analysis. How-ever, fits to such large data sets produce unwieldy datacovariance matrices and are typically not convergent, orrequire a prohibitively large number of iterations. Thiscan be partially addressed by thinning the data, e.g. bythe use of singular value decomposition (SVD) cuts, butthis reduces the accuracy of the fits. ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) ( , ) / ( , ) FIG. 2: (color online). Chained and marginalized fitresults for the ground state amplitude A (0 , V t , (1 , , of the B s → K decay on ensemble F2. Fit results are shown asa function of the number of ( top ) states explicitlyincluded in the fit and ( bottom ) total states accountedfor in the fit. Final results are taken from (6,1)/(8,8)fits, represented by gray bands.To address this problem we introduce a technique,which we refer to as chaining, to simplify fits to very largedata sets. Consider a data set consisting of N correlators,data = (correlator , correlator , . . . , correlator N ). Beforethe fit, all fit parameters are assigned priors. Chainingfirst fits correlator then uses the best fit mean values andcovariances to replace the corresponding priors in subse-quent fits. The updated set of priors is then used in thefit to correlator . In this and all subsequent fits, corre-lations are accounted for between the data being fit andthose priors which are best fit results from previous fits —this is an important step as it prevents “double counting”data. After this second fit, the priors are again updatedaccording to the best fit mean values and covariances.This process is repeated for all correlators. The collec-tion of best fit mean values and covariances following thefit to correlator N are the final fit results. Chaining isdescribed in greater detail in Appendix A.We combine the use of Bayesian [14], marginalized [1],and chained fitting techniques. Our final fit results use TABLE II: Fit results for the scalar and vector B s → K form factors on each ensemble and for each simulatedmomentum. Ensemble f B s K (000) f B s K (100) f B s K (110) f B s K (111)C1 0.8244(23) 0.7081(27) 0.6383(30) 0.5938(41)C2 0.8427(25) 0.6927(35) 0.6036(49) 0.536(12)C3 0.8313(29) 0.6953(33) 0.6309(30) 0.5844(46)F1 0.8322(25) 0.6844(35) 0.5994(43) 0.5551(56)F2 0.8316(27) 0.6915(38) 0.6119(43) 0.5563(61)Ensemble f B s K + (100) f B s K + (110) f B s K + (111)C1 2.087(16) 1.657(14) 1.378(13)C2 1.880(12) 1.412(16) 1.142(33)C3 1.773(11) 1.4212(84) 1.184(10)F1 1.878(13) 1.385(12) 1.158(13)F2 1.834(14) 1.396(10) 1.163(14) marginalization with a total of ( N, ˜ N ) = (8 ,
8) statesaccounted for, of which (6 ,
1) are explicitly fit. We re-fer to such fits with the shorthand notation, (6 , / (8 , B s → K form factor re-sults from the correlation function fits are tabulated inTable II and additional details are given in Appendix B.The form factors obtained from these fits preserve cor-relations resulting from shared gauge field configurationsand quark propagators used in data generation. Thepreservation of correlations is demonstrated in the toppanel of Fig. 3 where, e.g., significant correlations amongthe B s → K form factor fit results are seen at commonmomenta and nonzero correlations among form factorsfor the two decays is suggested. The bottom panel ofFig. 3 shows the distribution over all ensembles of corre-lations among form factors for the two decays. Account-ing for these correlations is useful in our determination ofthe ratio of form factors for the two decays. Fit resultsfor B s → η s , and the resulting form factor ratios, arepresented in Appendix D. V. CHIRAL, CONTINUUM, AND KINEMATICEXTRAPOLATION
The results of HPChPT [4, 5] suggest a factorization,to at least one-loop order, of the soft physics of logarith-mic chiral corrections and the physics associated withkinematics in the form factors describing semileptonicdecays of heavy mesons, f (cid:107) , ⊥ ( E ) = (1 + [logs]) K (cid:107) , ⊥ ( E ) . (19) f ( ) f ( ) f ( ) f ( ) f + ( ) f + ( ) f + ( ) f ( ) f ( ) f ( ) f ( ) f + ( ) f + ( ) f + ( ) f (000) f (100) f (110) f (111) f + (100) f + (110) f + (111) f (000) f (100) f (110) f (111) f + (100) f + (110) f + (111) B s → K B s → η s B s → K B s → η s cor ij FIG. 3: (color online). ( top ) Heat map of thecorrelation matrix for ensemble C1. ( bottom )Distribution of correlations among the form factors for B s → K and B s → η s for all ensembles.The logarithmic chiral corrections, calculated in Ref. [5]for several B ( s ) decays, are independent of E . An un-specified function K characterizes the kinematics.To obtain results over the full kinematic range onemust include lattice simulation data over a range of en-ergies. However, for any relevant physical scale Λ (e.g.Λ QCD , 1 /r , Λ ChPT , . . . ), E > ∼ Λ at nominal latticemomenta and there is no convergent expansion of theunknown function K ( E ) in powers of E/ Λ. This is aninherent limitation of characterizing the kinematics interms of energy. The energy of the daughter meson is apoor variable with which to describe the kinematics.In contrast, the z expansion [15–17] provides a con- vergent, model-independent characterization of the kine-matics over the entire kinematically accessible range.Combining a z expansion on each ensemble with theHPChPT inspired factorization of Eq. (19) allows a si-multaneous chiral, continuum, and kinematic extrapola-tion of lattice data at arbitrary energies. Because thechiral logs are the same for f (cid:107) and f ⊥ , linear combina-tions (i.e. f and f + ) factorize in the same way and havethe same chiral logs. Motivated by these observations,we construct a HPChPT-motivated modified z expan-sion, which we call the “HPChPT z expansion”, and fitthe lattice data of Tables II and IX, with accompanyingcovariance matrix, to fit functions of the form P , + ( q ) f , + ( q ) = (1 + [logs]) × K (cid:88) k =0 a (0 , +) k D (0 , +) k z ( q ) k , (20)where [logs] are the continuum HPChPT logs of Ref. [5],and generic analytic chiral and discretization effects areaccounted for by D k . Resonances above q but belowthe B s K production threshold, i.e. those in the range q < q < ( M B s + M K ) , are accounted for via theBlaschke factor, P = 1 − q /M . Though not observed,we allow for the possibility of a J P = 0 + state in P ,with choice of mass guided by Ref. [13]. Our fit resultsare insensitive to the presence of this state. The factor-ization suggested by HPChPT may not hold at higherorder [18] so we allow chiral analytic terms, which helpparametrize effects from omitted higher order chiral logs,to have energy dependence (i.e. to vary with k ).We note that Eq. (20) is the modified z expansion in-troduced in Refs. [2, 3], with the coefficients of the chirallogarithmic corrections fixed by the results of HPChPT.In the chiral and continuum limitslim m → m physical a → (1 + [logs]) a k D k = b k of Ref. [17] , (21)and Eq. (20) is equivalent to the Bourrely-Caprini-Lellouch parametrization [17] of the form factors.Following Ref. [17] we impose a constraint on a (+) K fromthe expected scaling behavior of f + ( q ) in the neighbor-hood of q . The resulting fit function for f + is P + ( q ) f + ( q , a )= (1+[logs]) K − (cid:88) k =0 a (+) k D (+) k ( a ) (cid:2) z ( q ) k − ( − k − K kK z ( q ) K (cid:3) . (22)We write f ( q , a ), z ( q ), and D k ( a ), explicitly exposingthe dependence on q and a . This is useful in explaining This assumes the general arguments on which the z expansion isbased hold for heavier than physical quark masses and at finitelattice spacing. the implementation of a second kinematic constraint weimpose on the form factors. At the kinematic endpoint q = 0, the continuum extrapolated form factors f and f + are equal, i.e. f (0 ,
0) = f + (0 , a (0)0 , a (0)0 D (0)0 (0) = − K (cid:88) k =1 a (0) k D (0) k (0) z (0) k + K − (cid:88) k =0 a (+) k D (+) k (0) (cid:2) z (0) k − ( − k − K kK z (0) K (cid:3) . (23)Imposing this constraint results in the fit function for f : P ( q ) f ( q , a ) = (1 + [logs]) × (cid:26) K (cid:88) k =1 a (0) k (cid:104) D (0) k ( a ) z ( q ) k − D (0)0 ( a ) D (0)0 (0) D (0) k (0) z (0) k (cid:105) + D (0)0 ( a ) D (0)0 (0) K − (cid:88) k =0 a (+) k D (+) k (0) (cid:2) z (0) k − ( − k − K kK z (0) K (cid:3)(cid:27) . (24)In the fit functions for f and f + , Eqs. (22) and (24), D k and [logs] are given by, D k = 1 + c ( k )1 x π + c ( k )2 (cid:0) δx π + δx K (cid:1) + c ( k )3 δx η s + d ( k )1 ( a/r ) + d ( k )2 ( a/r ) + e ( k )1 ( aE K ) + e ( k )2 ( aE K ) , (25)[logs] = − x π (log x π + δ F V ) − g x K log x K − g x η log x η , (26)with implicit indices in Eq. (25) specifying the scalaror vector form factor. We account for momentum-independent and momentum-dependent discretization ef-fects in D k . The values of aE K that enter the fit are thevalues from the simulation and are, of course, small. Fi-nite volume effects in the simulation are included via ashift δ F V in the pion log [19]. The infinite volume limitis taken by setting this shift to zero. Eq. (26) gives theHPChPT [5] result for the chiral logarithmic correctionto B s → K form factors. These expressions make use ofthe dimensionless quantities x π,K,η = M π,K,η (4 πf π ) , (27) δx π,K = ( M asqtad π,K ) − ( M HISQ π,K ) (4 πf π ) , (28) δx η s = ( M HISQ η s ) − ( M physical η s ) (4 πf π ) , (29)where M η = ( M π + 2 M η s ) /
3. We determine q and z on each ensemble using correlator fit results for meson masses and simulation momenta. Light and heavy quarkdiscretization effects are accommodated for by makingthe d ( k ) i mild functions of the masses, accomplished bythe replacements d ( k )1 → d ( k )1 (1 + l ( k )1 x π + l ( k )2 x π )(1 + h ( k )1 δx b + h ( k )2 δx b ) ,d ( k )2 → d ( k )2 (1 + l ( k )3 x π + l ( k )4 x π )(1 + h ( k )3 δx b + h ( k )4 δx b ) , (30)where δx b = am b − .
26 is chosen so that as am b variesover the coarse and fine ensembles − . < ∼ δx b < ∼ . O ( α s , Λ QCD /m b , α s / ( am b )) contributions to be ∼
4% ofthe total contribution to (cid:104) V (cid:105) . Of this 4% the major-ity, ∼ . O ( α s ) correctionand <
1% from the NRQCD matching via (cid:104) J (1) , sub0 (cid:105) . For (cid:104) V k (cid:105) we find contributions at this order to be ∼ ∼
1% coming from the O ( α s ) correction and <
1% fromthe NRQCD matching. The matching error results fromomitted higher order corrections, the size of which we es-timate from observed leading order effects, where we con-servatively use the larger 4%. Following the argumentsoutlined in Ref. [10] we estimate the matching error to bethe same size as the observed O ( α s , Λ QCD /m b , α s / ( am b ))contributions and take the matching error to be 4%. Thisis equivalent to taking the O ( α s ) matching coefficient tobe four times larger than the O ( α s ) matching coefficient ρ ( V )0 (13 times larger than ρ ( V k )0 ). This uncertainty isassociated with the hadronic matrix elements and there-fore, by Eqs. (3) and (4), with f (cid:107) and f ⊥ . To correctlyincorporate it in the results for f and f + we convert ourfit functions for f , + into f (cid:107) , ⊥ , multiply by (1 + m (cid:107) , ⊥ ),where m (cid:107) , ⊥ is a coefficient representing the matching er-ror with a prior central value of zero and width 0.04, thenconvert back to f , + before performing the fit. Schemat-ically, we modify the fit functions, defined in Eqs. (22)and (24), by f , f + → f (cid:107) , f ⊥ (31) f (cid:107) , f ⊥ → (1 + m (cid:107) ) f (cid:107) , (1 + m ⊥ ) f ⊥ (32)(1 + m (cid:107) ) f (cid:107) , (1 + m ⊥ ) f ⊥ → f corrected0 , f corrected+ , (33)then we use f corrected0 , + to fit the results of the correlationfunction fits of Sec. IV. Conversions between the formfactors f , + and f (cid:107) , ⊥ are performed using Eqs. (5) and(6).The results of a simultaneous fit to the data for f B s K , + and f B s η s , + , in which the maximum order of z [specified by K in Eqs. (22) and (24)] is 3 and χ / d . o . f . = 84 . /
70, areshown relative to the data in Fig. 4 for B s → K . Detailsof prior choices and fit results are given in Appendix C.We test the stability of this fit to the following modi-fications of the fit Ans¨atze:(1) Truncate the z expansion at O ( z ). zf + B s K f B s K C1C2C3 zf + B s K f B s K F1F2
FIG. 4: (color online). B s → K form factor results froma simultaneous chiral, continuum, and kinematicextrapolation via the HPChPT z expansion are shown( top ) relative to coarse ensemble data (C1, C2, and C3)and ( bottom ) relative to fine ensemble data (F1 and F2).(2) Truncate the z expansion at O ( z ).(3) Truncate the z expansion at O ( z ).(4) Drop O ( aE K ) momentum-dependent and O ( a )momentum-independent discretization terms inEq. (25).(5) Drop the am b -dependent discretization terms inEq. (30).(6) Drop the light-quark mass-dependent discretizationterms in Eq. (30).(7) Add the following next-to-next-to-leading-order(NNLO) chiral analytic terms to D k as defined in Eq. (25): c ( k )4 x π + c ( k )5 (cid:0) δx π + δx K (cid:1) + c ( k )6 δx η s + c ( k )7 x π (cid:0) δx π + δx K (cid:1) + c ( k )8 x π δx η s (34)+ c ( k )9 (cid:0) δx π + δx K (cid:1) δx η s + c ( k )10 x π ( a/r ) + c ( k )11 (cid:0) δx π + δx K (cid:1) ( a/r ) + c ( k )12 δx η s ( a/r ) . (8) Drop the sea- and valence-quark mass differenceterm (cid:0) δx π + δx K (cid:1) from Eq. (25).(9) Drop the strange quark mistuning term δx η s fromEq. (25).(10) Drop finite volume effects, i.e. set δ F V = 0 inEq. (26).The stability of the B s → K fit results to these modi-fications is shown in Fig. 5, where results are shown atthe extrapolated q = 0 point. This point is furthestfrom the data region where simulations are performedand therefore is particularly sensitive to changes in thefit function. In Fig. 5 our final fit result, as defined byEqs. (22) and (24) with K = 3 and by Eqs. (25)–(30), isindicated by the dashed line and gray band.Modifications 1, 2, and 3 vary the order of the trunca-tion in z and demonstrate that by O ( z ) fit results havestabilized and errors have saturated. We therefore con-clude that the error of the O ( z ) fit adequately accountsfor the systematic error due to truncating the z expan-sion.Momentum-dependent and momentum-independentdiscretization effects proportional to a are removed inmodification 4. This results in a modest increase in χ and a negligible shift in the fit result. This suggests our fi-nal fit, which includes the a effects, adequately accountsfor all discretization effects observed in the data.In modifications 5 and 6 we remove heavy- and light-quark mass-dependent discretization effects with essen-tially no impact on the fit. That our results are indepen-dent of light-quark mass dependent discretization effectssuggests that staggered taste violating effects are accom-modated for by a generic a dependence.Modification 7 tests the truncation of chiral analyticterms after next-to-leading-order (NLO) by adding theNNLO terms listed in Eq. (34). This results in a slightdecrease in χ but has no noticeable effect on the fitcentral value or error. From this we conclude that er-rors associated with omitted higher order chiral termsare negligible.Differences in sea and valence quark masses, due inpart to our use of HISQ valence- and asqtad sea-quarks,are neglected in modification 8. This results in a smallincrease in χ and negligible change in the fit results. Weaccount for these small mass differences in our final fit,though this test suggests they are unimportant in the fit.Effects due to strange quark mass mistuning on theensembles are omitted in modification 9, resulting in a
82 84 86 88 90 92 94 96 1 2 3 4 5 6 7 8 9 10 χ [ d . o . f . ] f , + B s K ( ) FIG. 5: (color online). The stability of the HPChPT z expansion is demonstrated by studying the fit resultsunder various modifications, discussed in the text. Thetop panel shows χ with 70 degrees of freedom (d.o.f.)for each test fit and the bottom panel shows formfactors extrapolated to q = 0.modest increase in χ and no change in the fit centralvalue and error. We include these effects in our final fit.Modification 10 results in nearly identical fit results,suggesting that finite volume effects are negligible in ourdata. We include these effects in our final fit results. VI. FORM FACTOR RESULTS
In this section we present final results, with a com-plete error budget, for the B s → K form factors. Weprovide the needed information to reconstruct the formfactors and compare our results with previous model cal-culations.Fig. 6 shows the results of the chiral, continuum, andkinematic extrapolation of Sec. V, plotted over the entirekinematic range of q . The form factors, extrapolated to q = 0, have the value f B s K , + (0) = 0 . A. Fit errors for the HPChPT z expansion The inputs in our chiral, continuum, and kinematicextrapolation fits are data (the correlator fit results for f and f + in Tables II and IX with the accompanyingcovariance matrix) and priors. The total hessian error ofthe fit can be described in terms of contributions fromthese inputs, as described in detail in Appendix A. Wegroup priors in a meaningful, though not unique, way and q [GeV ] f B s K f + B s K data 0 0.5 1 1.5 2 2.5 3 3.5 0 4 8 12 16 20 24 FIG. 6: (color online). B s → K form factor results froma simultaneous chiral, continuum, and kinematicextrapolation via the HPChPT z expansion. The q region for which lattice simulation data exist isindicated by the shaded region.discuss the error associated with the chiral, continuum,and kinematic extrapolation based on these groupings.As the priors are, by construction, uncorrelated with oneanother, we can group them together in any way we findmeaningful. The resulting error groupings are uncorre-lated and add in quadrature to the total error. In Fig. 7we plot the following relative error components as func-tions of q :(i) experiment: This is the error in the fit due to un-certainty of experimentally determined, and other,input parameters. It is the sum in quadrature ofthe errors due to priors for the “Group I” fit pa-rameters listed in Table VII. This error is inde-pendent of q and subdominant.(ii) kinematic: This error component is due to the pri-ors for the coefficients a (0 , +) k in Eqs. (22) and (24).A comparison of the fit results from modifications1, 2, and 3 in Fig. 5 shows that by O ( z ) the fitresults have stabilized and errors have saturated.The kinematic error therefore includes the errorassociated with truncating the z expansion. Theextrapolation to values of q for which we have nosimulation data is controlled by the z expansion.As a result, the growth in form factor errors awayfrom the simulation region is due almost entirelyto kinematic and statistical errors.(iii) chiral: This error component is the sum in quadra-ture of errors associated with priors for c ( k ) i inEq. (25). These terms are responsible for extrap-olating to the physical light quark mass and foraccommodating for the slight strange quark mis-tuning and the small mismatch in sea and valencequark masses due to the mixed action used in thesimulation. As shown in Fig. 7, these errors aresubdominant and do not vary significantly with q .(iv) discretization: We account for momentum-dependent discretization effects via the e ( k ) i , andmomentum-independent discretization effects viathe d ( k ) i , terms of Eq. (25). In addition we allow forheavy- and light-quark mass-dependent discretiza-tion effects via the h ( k ) i and l ( k ) i terms in Eq. (30).The discretization error component, which is es-sentially independent of q , is the sum in quadra-ture of the error due to the priors for these fitparameters.(v) statistical: The statistical component of the erroris due to uncertainty in the data, i.e. the errorsfrom form factor fit results of Table II. Simulationdata exist for q > ∼
17 GeV for f and over therange 17 GeV < ∼ q < ∼
22 GeV for f + . Extrap-olation beyond these regions leads to increasingerrors.(vi) matching: The matching error is due to the un-certainty associated with the priors for m (cid:107) , ⊥ intro-duced in Eq. (32) and discussed in the surroundingtext.In addition to the largest sources of error, which weaccount for directly in the fit, there are remaining sys-tematic uncertainties.We simulate with degenerate light quarks and neglectelectromagnetism. By adjusting the physical kaon mass( M K ± → M K ) used in the chiral, continuum and kine-matic extrapolation, we estimate the “kinematic” effectsof omitting electromagnetic and isospin symmetry break-ing in our simulation to be < ∼ . P , + , are shown in Fig. 8 where they arecompared with results from a model calculation usingperturbative QCD (pQCD) [21] and a relativistic quarkmodel (RQM) [22]. Our results provide significant clari-fication on the form factors at large q . B. Reconstructing B s → K(cid:96)ν
Form Factors
In the physical limit our form factor results areparametrized in a BCL [17] form with coefficients b (0 , +) k [see Eq. (21)]. Including the kinematic constraint and f r e l a ti v e e rr o r c o m pon e n t s [ % ] q [GeV ] experimentkinematicchiraldiscretizationstatisticalmatchingtotal 0 4 8 12 16 20 0 4 8 12 16 20 24 f + r e l a ti v e e rr o r c o m pon e n t s [ % ] q [GeV ] experimentkinematicchiraldiscretizationstatisticalmatchingtotal 0 4 8 12 16 20 0 4 8 12 16 20 24 FIG. 7: (color online). B s → K ( top ) f and ( bottom ) f + relative error components. The total error (solidline) is the sum in quadrature of the components.terms through order z , we have P ( q ) f ( q ) = (cid:88) k =1 b (0) k ( z k − z (0) k )+ (cid:88) k =0 b (+) k (cid:20) z (0) k − ( − k − k z (0) (cid:21) , (35) P + ( q ) f + ( q ) = (cid:88) k =0 b (+) k (cid:20) z k − ( − k − k z (cid:21) , (36)where z ( q ) = (cid:112) t + − q − √ t + − t (cid:112) t + − q + √ t + − t , (37) t + = ( M B s + M K ) , (38) t = ( M B s + M K )( (cid:112) M B s − (cid:112) M K ) , (39) P , + ( q ) = 1 − q /M , + , (40)and the resonance masses are M = 5 . M + = 5 . Top ) Physical extrapolated coefficients of the HPChPT z expansion for the B s → K form factors,defined in Eqs. (35) and (36) and ( bottom ) the associated covariance matrix. Coefficient Value b (0)1 b (0)2 b (0)3 b (+)0 b (+)1 -0.750(193) b (+)2 b (0)1 b (0)2 b (0)3 b (+)0 b (+)1 b (+)2 b (0)1 . × − . × − . × − . × − . × − . × − b (0)2 .
702 5 .
852 9 . × − . × − . b (0)3 . × . × − . × − . b (+)0 . × − . × − − . × − b (+)1 . × − . × − b (+)2 . b (0 , +) k , derived from the extrapolation fit results of Sec. V,and the associated covariance matrix, are given in Ta-ble III. Note that it is necessary to take into account thecorrelations among the coefficients to correctly reproducethe form factor errors. VII. PHENOMENOLOGY
With the benefit of ab initio form factors from latticeQCD, we explore the standard model implications of ourresults. In this section we make standard model predic-tions for several observables related to the B s → K(cid:96)ν decay for (cid:96) = µ and τ .The standard model B s → K(cid:96)ν differential decay rateis related to the form factors by d Γ dq = G F | V ub | π M B s (cid:16) − m (cid:96) q (cid:17) | p K | (cid:20)(cid:16) m (cid:96) q (cid:17) M B s p K | f + | + 3 m (cid:96) q ( M B s − M K ) | f | (cid:21) . (41)In Fig. 9 we plot predicted differential decay rates for B s → Kµν and B s → Kτ ν , divided by | V ub | , over thefull kinematic range of q . The ratio Γ / | V ub | can becombined with experimental results for the decay rates,typically differential decay rates integrated over q bins,to allow the determination of | V ub | . In Eqs. (42) and (43)we give numerical results for d Γ /dq , integrated over thekinematically accessible regions of q ,Γ( B s → Kµν ) / | V ub | = 7 . .
52) ps − , (42)Γ( B s → Kτ ν ) / | V ub | = 4 . .
60) ps − . (43) Combining our form factor results with the current inclusive and exclusive semileptonic determinations of | V ub | , exclusive | V ub | = 3 . × − , (44)inclusive | V ub | = 4 . × − , (45)we demonstrate in Fig. 10 the potential of this decay toshed light on this ∼ σ discrepancy. In this and subse-quent figures, dark interior bands represent the error inthe differential branching fractions omitting the error as-sociated with | V ub | . Experimental errors commensuratewith these predictions, especially for the B s → Kτ ν de-cay or at large q for the B s → Kµν decay, would allowdifferentiation between the current inclusive and exclu-sive values of | V ub | .Decays that couple to the τ have increased dependenceon the scalar form factor and to new physics models withscalar states (see, e.g., Refs. [25, 26] for a discussion ofnew physics in the closely related decay B → πτ ν ). Theratio of the B s → Kτ ν differential branching fraction tothat for B s → Kµν , R τµ ( q , q ) = (cid:82) q q dq d B /dq ( B s → Kτ ν ) (cid:82) q q dq d B /dq ( B s → Kµν ) , (46)is therefore a potentially sensitive probe of new physics.Integrating over the full kinematic range, we find R τµ ( m µ , q ) = 0 . , (47) For inclusive | V ub | we take the value from the Particle DataGroup [23]. For the exclusive determination we use the “globallattice + Belle” results reported by the FLAG-2 collabora-tion [24]. P f B s K q [GeV ] this workpQCDRQM 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 4 8 12 16 20 24 P + f + B s K q [GeV ]this workpQCDRQM 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 4 8 12 16 20 24 FIG. 8: (color online). Comparison of our B s → K ( top ) f and ( bottom ) f + form factors with those from apQCD model [21] and the RQM [22]. To easecomparison the vertical scale is reduced by multiplyingthe form factors by the Blaschke factor P , + .where q = ( M B s − M K ) . We plot the standard modelprediction for this ratio, as a function of q = ( q + q ) /
2, over the full kinematic range in Fig. 11.The angular dependence of the differential decay rate,neglecting final state electromagnetic interactions, isgiven by d Γ dq d cos θ (cid:96) = G F | V ub | π M B s (cid:16) − m (cid:96) q (cid:17) | p K |× (cid:104) M B s | p K | (cid:16) sin θ (cid:96) + m (cid:96) q cos θ (cid:96) (cid:17) f + 4 m (cid:96) q ( M B s − M K ) M B s | p K | cos θ (cid:96) f f + + m (cid:96) q ( M B s − M K ) f (cid:105) , (48)where θ (cid:96) is defined, in the q rest frame (i.e. where p (cid:96) + p ν is zero), as the angle between the final state lepton and | V ub | - d Γ / dq ( B s → K µ ν ) [ p s - G e V - ] q [GeV ]0.00.10.20.30.40.50.60.70.8 0 4 8 12 16 20 24 | V ub | - d Γ / dq ( B s → K τ ν ) [ p s - G e V - ] q [GeV ]0.000.050.100.150.200.250.300.350.40 0 4 8 12 16 20 24 FIG. 9: (color online). Predicted differential decayrates, divided by | V ub | , for ( top ) B s → Kµν and( bottom ) B s → Kτ ν .the B s meson. From this angular dependence we canextract a forward-backward asymmetry [27], A (cid:96) FB ( q ) = (cid:34) (cid:90) − (cid:90) − (cid:35) d cos θ (cid:96) d Γ dq d cos θ (cid:96) (49)= G F | V ub | π M B s (cid:16) − m (cid:96) q (cid:17) | p K | × m (cid:96) q ( M B s − M K ) f f + , (50)which is suppressed in the standard model by a factor of m (cid:96) /q . In Fig. 12 we show standard model predictionsfor the forward-backward asymmetry using the inclusiveand exclusive values for | V ub | . Integrating over the fullkinematic range of q gives (cid:90) q m µ dq A µ FB ( q ) / | V ub | = 0 . − , (51) (cid:90) q m τ dq A τ FB ( q ) / | V ub | = 1 . − . (52)2 d B / dq ( B s → K µ ν ) [ G e V - ] q [GeV ]inclusive | V ub | exclusive | V ub | d B / dq ( B s → K τ ν ) [ G e V - ] q [GeV ] inclusive | V ub | exclusive | V ub | FIG. 10: (color online). Predicted differential branchingfractions for the ( top ) B s → Kµν and ( bottom ) B s → Kτ ν decays using inclusive and exclusivesemileptonic determinations of | V ub | . In each band, thelight outer band includes all sources of error and thedark interior band neglects the uncertainty in | V ub | .Normalizing the forward-backward asymmetry by thedifferential decay rate removes | V ub | ambiguity and mosthadronic uncertainties,¯ A (cid:96) FB ( q , q ) = (cid:82) q q dq A (cid:96) FB ( q ) (cid:82) q q dq d Γ /dq , (53)and represents the probability the lepton will have a mo-mentum component, in this frame, in the direction of mo-tion of the parent B s meson. Integrating over q yields¯ A µ FB ( m µ , q ) = 0 . , (54)¯ A τ FB ( m τ , q ) = 0 . , (55)with central values equal to those obtained by takingthe ratio of results from Eqs. (51) and (52) with those q [GeV ] R τµ FIG. 11: (color online). Predicted differential branchingfraction ratio.from Eqs. (42) and (43). The errors, however, are ∼ × smaller when correlations are accounted for. The normal-ized standard model asymmetries are plotted in Fig. 13as a function of q .The production of right-handed final state leptons ishelicity-suppressed in the standard model, providing aprobe of new physics via helicity-violating interactions.The standard model differential decay rates for left-handed (LH) and right handed (RH) polarized final stateleptons in B s → K(cid:96)ν decays is [27] d Γ(LH) dq = G F | V ub | | p K | π (cid:16) − m (cid:96) q (cid:17) f ,d Γ(RH) dq = G F | V ub | | p K | π m (cid:96) q (cid:16) − m (cid:96) q (cid:17) (56) × (cid:34)
38 ( M B s − M K ) M B s f + 12 | p K | f (cid:35) , and the (cid:96) -polarization distribution is given by the differ-ence A (cid:96) pol ( q ) = d Γ(LH) dq − d Γ(RH) dq . (57)We plot the τ -polarization distribution, again using theinclusive and exclusive values of | V ub | from Eqs. (44)and (45), in Fig. 14. Because of their relatively smallmass, muons produced in the decay are predominantlyleft-handed and the plot of A µ pol is equivalent to the to-tal differential decay rate. Integrating the (cid:96) -polarizationdistributions over q gives (cid:90) q m µ dq A µ pol ( q ) / | V ub | = 7 . .
60) ps − , (58) (cid:90) q m τ dq A τ pol ( q ) / | V ub | = 0 . − . (59)3 A µ FB [ p s - G e V - ] q [GeV ]inclusive | V ub | exclusive | V ub | A τ FB [ p s - G e V - ] q [GeV ] inclusive | V ub | exclusive | V ub | FIG. 12: (color online). Differential decay rateforward-backward asymmetries for the ( top ) B s → Kµν and ( bottom ) B s → Kτ ν decays using inclusive andexclusive semileptonic determinations of | V ub | . Lightouter bands includes all sources of error and the darkinterior bands neglect uncertainty in | V ub | .As with the forward-backward asymmetry, we normal-ize the (cid:96) -polarization distribution by the differential de-cay rate to remove ambiguity associated with | V ub | andhadronic uncertainties. The resulting polarization frac-tion [27] is defined by¯ A (cid:96) pol ( q , q ) = (cid:82) q q dq A (cid:96) pol ( q ) (cid:82) q q dq d Γ /dq . (60)Integrating over q we find the standard model predictionfor the fraction of polarized leptons to be¯ A µ pol ( m µ , q ) = 0 . +18 − ) , (61)¯ A τ pol ( m τ , q ) = 0 . , (62)where the error associated with the numerical integra-tion of ¯ A µ pol ( ± . q [GeV ] − A µ FB q [GeV ] − A τ FB FIG. 13: (color online). Normalized differential decayrate forward-backward asymmetries for the ( top ) B s → Kµν and ( bottom ) B s → Kτ ν decays. A τ po l [ p s - G e V - ] q [GeV ]inclusive | V ub | exclusive | V ub | -1.0-0.50.00.51.01.52.0 0 4 8 12 16 20 24 FIG. 14: (color online). Standard model τ -polarizationdistribution for the differential decay rate of B s → Kτ ν .constraint that ¯ A (cid:96) pol <
1. The q dependence of the (cid:96) -4 q [GeV ] − A µ pol -1-0.8-0.6-0.4-0.2 0 0.2 0.4 0 4 8 12 16 20 24 q [GeV ] − A τ pol FIG. 15: (color online). Standard model (cid:96) -polarizationfraction for the differential decay rate of B s → K(cid:96)ν , for (cid:96) = µ, τ .polarization fraction is plotted in Fig. 15. VIII. SUMMARY AND OUTLOOK
Using NRQCD b and HISQ light and strange valencequarks with the MILC 2 + 1 dynamical asqtad configu-rations, we report on the first lattice QCD calculation ofthe form factors for the semileptonic decay B s → K(cid:96)ν .With the help of a new technique, called chaining, wefit the B s → K correlator data simultaneously with datafor the fictitious decay B s → η s . Fitting these data si-multaneously accounts for correlations — useful for con-structing ratios of form factors. We extrapolate our lat-tice form factor results to the continuum, to physicalquark mass, and over the full kinematic range of q usinga combination of the modified z expansion and HPChPTthat we refer to as the HPChPT z expansion.We then make standard model predictions for:(i) differential decay rates divided by | V ub | , an ob-servable that, when combined with experiment, will allow an alternative semileptonic exclusive de-termination of | V ub | ,(ii) differential branching fractions using both the in-clusive and exclusive semileptonic B → π(cid:96)ν deter-minations of | V ub | ,(iii) the ratio of differential branching fractions R τµ ( q ),(iv) the forward-backward asymmetry, using inclusiveand exclusive values of | V ub | ,(v) the normalized forward-backward asymmetry,(vi) the τ -polarization distribution in the differentialdecay rate for B s → Kτ ν , and(vii) the (cid:96) -polarization fraction in the differential decayrate for B s → K(cid:96)ν , for (cid:96) = µ, τ .In Appendix D we construct ratios of form factors for B s → K with those for B s → η s . In combination witha future calculation of B s → η s using HISQ b , theseratios can provide a nonperturbative determination of the b → u current matching factor. This would be relevantfor both B s → K(cid:96)ν and B → π(cid:96)ν simulations usingNRQCD b quarks.Our results, built on first principles lattice QCD formfactors, greatly clarify standard model expectations [27]based on model estimates of form factors [21, 22, 28],most notably at large q . Combining our form factors,which are most precise at large q , with model calcula-tions, typically more reliable at low q , would result in amore precise determination of f and f + . We are study-ing the possibility of further refining B s → K(cid:96)ν standardmodel predictions using such form factors.
Acknowledgements
This research was supported by the DOE and NSF. Wethank the MILC collaboration for making their asqtad N f = 2 + 1 gauge field configurations available. Compu-tations were carried out at the Ohio Supercomputer Cen-ter and on facilities of the USQCD collaboration fundedby the Office of Science of the U.S. DOE. Appendix A: Fitting Basics
Here we describe in more detail two aspects of ourstatistical analysis: 1) the definition of our error budgetsfor fit results; and 2) the technique for chained fits ofmultiple data sets. We also discuss a general procedurefor testing fit procedures. These are general techniquesapplicable to many types of fitting problems [29]. Finallywe illustrate these ideas with an example drawn from thispaper.
1. Fits and Error Budgets
The formal structure of a least-squares problem in-volves fitting input data y i with functions f i ( p ) by ad-5justing fit parameters p α to minmize χ ( p ) = (cid:88) ij ∆ y ( p ) i (cid:0) cov − y (cid:1) ij ∆ y ( p ) j , (A1)where cov ij is the covariance matrix for the input dataand ∆ y ( p ) i ≡ f i ( p ) − y i . (A2)There are generally two types of input data — actualdata, and prior data for each fit parameter — but welump these together here since they enter χ ( p ) in thesame way. So the sums here over i and j are over alldata and priors. Note that priors and data may be cor-related in some problems.The best-fit parameters p α are those that minimize χ : ∂ α χ ( p ) = 2 (cid:88) ij ∂ α f i ( p ) (cid:0) cov − y (cid:1) ij ∆ y ( p ) j = 0 (A3)where the derivative ∂ α ≡ ∂/∂p α . The inverse covariancematrix, ∂ α ∂ β χ ( p ) /
2, for the p α is then given by (cid:0) cov − p (cid:1) αβ = (cid:88) ij ∂ α f i ( p ) (cid:0) cov − y (cid:1) ij ∂ β f j ( p ) + O (∆ y ) , (A4)where we neglect terms proportional to ∆ y (which makessense for reasonable fits to accurate data). This is theconventional result.The uncertainties in the p α are due to the uncertain-ties in the input data y i , and, for very accurate data, de-pend linearly upon cov y . The relationship can be demon-strated by differentiating Eq. (A3) with respect to y j toobtain (cid:88) β (cid:0) cov − p (cid:1) αβ ∂p β ∂y j = (cid:88) i ∂ α f i ( p ) (cid:0) cov − y (cid:1) ij + O (∆ y ) , (A5)where again we neglect terms proportional to ∆ y . Solv-ing for ∂p β /∂y j gives: ∂p β ∂y j = (cid:88) αi (cov p ) βα ∂ α f i ( p ) (cid:0) cov − y (cid:1) ij (A6)In the high-statistics, small-error limit the covariancesin the p α are related to those in the y i by the standardformula (cov p ) αβ = (cid:88) ij ∂p α ∂y i (cov y ) ij ∂p β ∂y j , (A7)and, indeed, substituting Eq. (A6) into this equation re-produces Eq. (A4) for cov p .Eqs. (A6) and (A7) allow us to express the error σ g fora function g ( p ) of the best-fit parameter values in termsof the input errors: σ g ≡ (cid:88) αβ ∂ α g ( p ) (cov p ) αβ ∂ β g ( p ) = (cid:88) ij c ij (cov y ) ij (A8) where c ij ≡ (cid:88) αβ ∂ α g ( p ) ∂p α ∂y i ∂p β ∂y j ∂ β g ( p ) . (A9)and Eq. (A6) is used to evaluate ∂ ¯ p α /∂y i . We can thendecompose σ g into separate contributions coming fromthe different block-diagonal submatrices of cov y . Thesecontributions to σ g constitute the error budget for g ( p ).The c ij s in Eq. (A8) depend upon both the y i and theircovariance matrix, but that dependence can be neglectedto leading order in cov y . Consequently Eq. (A8) can beused to estimate the impact on σ g of possible modifica-tions to any element of cov y .Note that the data’s covariance matrix cov y can bequite singular if there are strong correlations in the data.This can make it numerically difficult to invert the matrixfor use in χ ( p ). This problem is typically dealt with byusing a singular value decomposition (SVD) to regulatethe most singular components of the covariance matrix.In our fits we rescale the covariance matrix by its diago-nal elements to obtain the correlation matrix, which wethen diagonalize. We introduce a minimum eigenvalueby setting any smaller eigenvalue equal to the minimum.We then reconstitute the correlation matrix, and rescaleit back into a (less singular) covariance matrix which weuse in the fit. This procedure, in effect, increases the er-ror in the data and so increases the uncertainties in thefinal fit results; it is a conservative move.It is common when using SVD to discard eigenmodescorresponding to the small eigenvalues. This is equivalentto setting the variance associated with these modes toinfinity in the fit. In our implementation, all eigenmodesare retained, but the small eigenvalues are replaced bya (larger) minimum eigenvalue. This is a more realisticestimate for the variances of these modes — that is, morerealistic than setting them to infinity — and gives moreaccurate fit results.
2. Chained Fits
Chained fits simplify fits of multiple data sets whose fitfunctions share fit parameters by allowing us to fit eachdata set separately. To illustrate, consider two sets ofdata, y i ( A ) and y j ( B ), that we fit with functions f i ( A, p )and f j ( B, p ), respectively — both functions of the samefit parameters p α (unlike the previous section, here wedo not lump the priors in with the y s). The fit proce-dure is straightforward in a Bayesian framework if y ( A )and y ( B ) are statistically uncorrelated. We first fit, say,data set y ( A ) to obtain best-fit estimates p ( A ) for theparameters and an estimate cov p ( A ) for the parameters’covariance matrix. We then fit data set y ( B ), but using p ( A ) and cov p ( A ) to form the prior for the fit parameters.This two-step fit merges the information containedin y ( A ) with that from y ( B ) by feeding the informationfrom the first fit into the second fit as prior information.6The order in which the data sets are fit doesn’t matter inthe high-statistics (Gaussian) limit; with larger errors, itis better to fit the more accurate data set first. The χ for the two-step fit is the sum of the χ s for each step.The situation is slightly more complicated if y ( A ) and y ( B ) are correlated. Then the best-fit parameters p ( A )from the first fit above are correlated with the seconddata set y ( B ). The p ( A )- y ( B ) covariance can be com-puted fromcov p ( A ) y ( B ) ≡ (cid:88) y ( A ) ∂p ( A ) ∂y ( A ) cov y ( A ) y ( B ) (A10)using Eq. (A6) in the previous section. This correlationmust be included in the second fit, to data set y ( B ).So the second fit uses the best-fit parameters p ( A ) fromthe first fit to construct the prior, together with cov p ( A ) for parameter-parameter covariances and cov p ( A ) y ( B ) forparameter-data covariances.We refer to a sequential fit of multiple data sets, wherethe best-fit parameters and covariance matrix from onefit are used as the prior for the next fit, as a chainedfit . It is essential in such fits to account for possiblecorrelations between the priors (from previous fits) andthe data being fit at each stage. The results of a chainedfit should agree with those of a simultaneous fit in thelimit of large (i.e., Gaussian) statistics.
3. Testing Fits
It is generally useful to have ways of testing particularfit strategies. One simple approach to testing is to cre-ate multiple fake data sets that are very similar to theactual data being fit, but where the exact values for thefit parameters are known ahead of time. Running severalsuch data sets through an analysis code tells you veryquickly whether, for example, your analysis code givesresults that are correct to within one sigma 68% of thetime, as is desired.It is easy to create fake data sets of this sort. Onesimple recipe is the following:1. Fit the actual data to obtain a set of parametervalues p ∗ α such that the fit function f i ( p ∗ ) closelymatches the mean values y i of the actual data. Cal-culate the difference between the actual means ofthe data and the fit values for p = p ∗ : δy i ≡ f i ( p ∗ ) − y i . (A11)2. Create a bootstrap copy y bs i of the original dataand replace its mean values by: y ∗ i = y bs i + δy i (A12)The fake data set then consists of the mean values y ∗ i and the covariance matrix cov y of the originaldata. The role of the bootstrap here is to generate fluctuations in the means with the same distribu-tion as the original data. These data sets will fluc-tuate around central values f i ( p ∗ ) rather than theoriginal means of the data.3. Repeat the second step to create any number ofadditional fake data sets.Each fake data set is fit using the same procedure thatwas used to analyze the original data. The results for thefit parameters are compared with the parameter values p ∗ used to define the correction δy i [Eq. (A11)], since, byconstruction, these are the correct values for the param-eters in the fake data.Typically only a handful of parameters from a fit areof interest. Their best-fit values from different fake datasets will differ, but they should all agree with the p ∗ val-ues to within the errors generated by the fake fit (thatis, to within one sigma 68% of the time, two sigma 95%of the time, and so on). Such tests can reveal, for exam-ple, potential problems coming from poor priors or inad-equate SVD cuts, or biases in particular combinations offit parameters.
4. Example
We compare chained and unchained fit results inFig. 16. Because unchained fits to very large data setsare unreliable, for purposes of comparison we divide thedata into the smallest subsets that allow the extractionof individual matrix elements. Such fits are uncorrelatedin that they neglect correlations among data at differentmomenta, for different currents, and among the two de-cays. The uncorrelated fits include only one decay mode( B s → K or B s → η s ), data for only one simulationmomentum (000, 100, 110, or 111), and only one cur-rent ( V t or V k ). These fits are still complicated, however,as they require the minimum amount of data needed toextract a single matrix element. This minimum num-ber of correlators consists of parent and daughter twopoint and three point data, i.e. B s → B s , η s → η s (000),and B s → V t → η s (000). Including correlations resultsin marked improvement in the accuracy of matrix ele-ments obtained from the noisiest data — that for V k atlarge momenta. This improvement can be traced to cor-relations of these data with the more precise data for V t (for the same decay and at a common momentum), asdemonstrated in Fig. 3.In addition to properly accounting for correlations inthe data, chaining reduces the time required to performthe fits. While the uncorrelated fits required a total of1 hour 14 minutes, the chained (8,8) fit required only24 minutes. The use of marginalization significantly re-duces the time required. The chained and marginalized(6,1)/(8,8) fit required only 57 seconds.7 -2-1 0 1 2 V t (000) V t (100) V t (110) V t (111) V k (100) V k (110) V k (111)uncorrelated (6,6)chained (8,8)chained, marginalized (6,1)/(8,8) FIG. 16: (color online). B s → η s three point ground state amplitudes, for varying currents and momenta, asobtained from different fitting strategies described in the text. Plotted central values indicate the number ofstandard deviations by which a fit result differs from an “uncorrelated” fit. The size of the error bars is the ratio ofthe plotted fit error to that from an uncorrelated fit.TABLE IV: B s priors and fit results for aE sim(0) B s . Ensemble Prior 2pt 2+3ptC1 0.537(53) 0.53780(72) 0.53801(31)C2 0.54(6) 0.54360(84) 0.54234(35)C3 0.54(8) 0.5362(15) 0.53575(36)F1 0.405(55) 0.4081(13) 0.40869(21)F2 0.407(60) 0.40770(64) 0.40710(23)
Appendix B: Correlator Fit Results
The method for selecting priors for correlator fits wasdescribed in detail in Appendix B of Ref. [10]. We use thesame method in this analysis. Tables IV, V, and VI tabu-late priors and fit results for ground state energies. Theycompares results obtained from fits to two point correla-tion function data to those from simultaneous fits to twoand three point correlation function data, as describedin Sec. IV. The combined fits show improved precisionfor the B s meson mass and the larger momenta daughtermeson energies, suggesting that the three point correla-tion function data provide additional information to thefit. Within errors, the two point and simultaneous twoand three point fit results are consistent. Appendix C: HPChPT z Expansion Fit Results
Group I parameters listed in Table VII insert error inthe fit based on uncertainty associated with input param-eters – quantities not determined by the data. Priors for r and M η phys s are taken from Ref. [30]. We base our priorchoice for the BB ∗ π coupling g BB ∗ π on the combinedworks in Ref. [31]. Resonance masses for the Blaschke factors P , + introduced in Eq. (20) are calculated rela-tive to the B s meson mass in our simulations, M B s K = M B s − ( M B s − M B ) + 400(1) MeV , (C1) M B s K + = M B s − ( M B s − M B ) + ∆ hyperfine B , (C2) M B s η s = M B s + 400(1) MeV , (C3) M B s η s + = M B s + ∆ hyperfine B s , (C4)and we refer to the shift relative to M B s as ∆ B s K,B s η s , + .The M B s − M B and hyperfine splittings are taken fromthe PDG [23]. We tested increasing the uncertainty inthe location of the scalar pole, which we have taken tobe 400(1) MeV above the J P = 0 − state. A splitting of400(50) MeV gives identical results for the form factors,in both the central value and error, but accommodatesfor part of the error in f via allowed uncertainty in M .To reconstruct the form factors in this case, correlationsbetween P and the coefficients of the z expansion mustbe accounted for. By effectively fixing M we arrive atthe same fit results and can neglect uncertainty in P andcorrelations with the coefficients. The 4% uncertainty as-sociated with the perturbative matching is accounted forby m (cid:107) and m ⊥ , where we use prior central values of zeroand width 0.04, as explained by Eq. (32) and surround-ing text. Matrix elements for B s → K and B s → η s usethe same matching factors so we use common m (cid:107) , ⊥ forboth data sets. We use values for r /a from Ref. [6] and M asqtad π,K from Ref. [32]. We use values for M HISQ π,K,η s and M B s from best fit results in this and an ongoing B → π analysis using HISQ valence quarks.The Group II parameters of Table VIII are quantitiesdetermined by the fit. We choose priors for a k to be0 ±
5, based roughly on the unitarity constraint, and ver-ified that fit results are insensitive to variations in theprior width from 1 to 10. Chiral analytic terms are writ-ten in terms of dimensionless parameters that are nat-urally O (1). For this reason we use priors of zero with8TABLE V: K priors and fit results. For each ensemble, the first row lists priors, the second row gives two pointcorrelator fit results, and the third row shows simultaneous two and three point correlator fit results. Ensemble aM (0) K aE (0) K (100) aE (0) K (110) aE (0) K (111) C1 0.312(17) 0.41(11) 0.48(23) 0.55(28)0.31211(15) 0.40657(58) 0.48461(76) 0.5511(16)0.31195(14) 0.40661(49) 0.48408(63) 0.5513(13)C2 0.329(24) 0.45(15) 0.55(15) 0.61(31)0.32863(18) 0.45406(85) 0.5511(16) 0.6261(75)0.32870(16) 0.45434(73) 0.5506(11) 0.6273(35)C3 0.356(25) 0.475(75) 0.58(20) 0.65(30)0.35717(22) 0.47521(85) 0.5723(11) 0.6524(30)0.35744(21) 0.47507(71) 0.57218(80) 0.6539(18)F1 0.229(60) 0.32(24) 0.39(34) 0.43(40)0.22865(11) 0.32024(66) 0.39229(86) 0.4515(25)0.22861(12) 0.32020(61) 0.39192(82) 0.4528(16)F2 0.246(36) 0.33(23) 0.40(30) 0.47(37)0.24577(13) 0.33322(52) 0.40214(73) 0.4623(14)0.24566(13) 0.33310(50) 0.40184(72) 0.4624(11)
TABLE VI: Like Table V but for the η s . Ensemble aM (0) η s aE (0) η s (100) aE (0) η s (110) aE (0) η s (111) C1 0.411(9) 0.487(12) 0.553(50) 0.61(11)0.41111(12) 0.48736(23) 0.55311(29) 0.61148(60)0.41107(11) 0.48726(23) 0.55294(29) 0.61135(52)C2 0.415(12) 0.52(5) 0.61(11) 0.68(23)0.41445(17) 0.51949(46) 0.6063(12) 0.6797(31)0.41446(15) 0.51934(44) 0.60647(67) 0.6794(18)C3 0.412(20) 0.518(40) 0.61(12) 0.69(35)0.41180(23) 0.51757(63) 0.60723(78) 0.6831(23)0.41175(20) 0.51742(57) 0.60720(67) 0.6843(14)F1 0.294(24) 0.37(10) 0.43(23) 0.48(34)0.294109(93) 0.36965(31) 0.43278(45) 0.4867(13)0.294066(88) 0.36988(26) 0.43301(38) 0.48729(88)F2 0.293(30) 0.369(89) 0.43(18) 0.49(30)0.29315(12) 0.36939(35) 0.43259(45) 0.48810(87)0.29310(12) 0.36927(35) 0.43197(48) 0.48729(97) width one for c and c . Based on previous analyses usingthe same ensembles we know that sea-quark effects aresmaller than those of the valence quarks, so we choosepriors for c to be 0 ± .
3. The leading order HISQ dis-cretization effects are O ( α s a ), so for the coefficients d and e which characterize the O ( a ) discretization ef-fects, we choose priors of 0 ± .
3. Coefficients d and e characterize O ( a ) effects and we use 0 ±
1. The coeffi-cients h and l characterize light- and heavy-quark mass-dependent discretization effects. These terms are writtenin terms of O (1) quantities and we take the coefficientsto have priors of 0 ± Appendix D: B s → η s Form Factors and Ratios
The results of B s → η s correlator fits are tabulated inTable IX and plotted as data points in the top two panelsof Fig. 17. From these plots one sees that simulationdata exhibit very small light sea quark mass and latticespacing dependence. These fit results are obtained froma single fit to both the B s → K and B s → η s datadescribed in Sec. IV. As a result, the B s → η s fit resultsof Table IX are correlated with the B s → K results ofTable II, as shown in Fig. 3.The B s → η s form factor data of Table IX is extrap-olated to the physical quark mass, the continuum limit,and over the entire kinematic range using the HPChPT z expansion described in Sec. V. This fit is also done si-9multaneously with the extrapolation of the B s → K data.The fit functions for the simultaneous chiral, continuum,and kinematic extrapolation of B s → η s are equivalentto those of Sec. V, with Eqs. (25) and (26) modified asfollows: D k = 1 + c ( k )1 x K + c ( k )2 δx K + c ( k )3 δx η s + d ( k )1 ( a/r ) + d ( k )2 ( a/r ) + e ( k )1 ( aE η s ) + e ( k )2 ( aE η s ) , (D1)[logs] = − g x K log x K − g x η log x η , (D2)with implicit indices in Eq. (D1) specifying scalar or vec-tor form factor. Results of this fit for the B s → η s formfactors are shown relative to data, and extrapolated overthe full kinematic range of q , in Fig. 17. The HPChPT z expansion stability analysis outlined in Sec. V involvedsimultaneous fits to both B s → K and B s → η s data.The B s → η s fit results for each of the modifications dis-cussed in that analysis are shown in Fig. 18. Becausethese results are from a simultaneous fit, the values of χ in Fig. 5 are applicable here as well and are reproducedfor convenience in Fig. 18. Note that the chiral analyticterms for B s → η s differ slightly from those for B s → K ,c.f. Eqs. (25) and (D1). As a result, the NNLO analyticterms added to the B s → η s fit function in modification7 differ from those listed in Eq. (34). Error breakdownplots for the B s → η s form factors are shown in Fig. 19.In the ratios of form factors, R (cid:107) ( q ) = f B s K (cid:107) ( q ) f B s η s (cid:107) ( q ) , (D3) R ⊥ ( q ) = f B s K ⊥ ( q ) f B s η s ⊥ ( q ) , (D4)the leading systematic error, that due to one-loop pertur-bative matching, largely cancels. Fig. 20 plots the ratiosas functions of q and shows that they are most preciselydetermined at q = ( M B s − M η s ) , where R (cid:107) (cid:0) ( M B s − M η s ) (cid:1) = 0 . , (D5) R ⊥ (cid:0) ( M B s − M η s ) (cid:1) = 0 . . (D6)The errors of the ratios are broken down into componentsin Fig. 21. Neglecting correlations among the B s → K and B s → η s decays yields ratios at this q with ∼ f B s η s (cid:107) , ⊥ using HISQ b quarks, these ratios will provide anonperturbative determination of the NRQCD b → u current matching factor, applicable to both B s → K and B → π . TABLE VII: Group I priors for the HPChPT z expansion for f B s K , + and f B s η s , + . Quantities listed infive consecutive rows have ensemble-dependent valuescorresponding to C1, C2, C3, F1, and F2. Group I Prior Fit r [fm] 0.3133(23) 0.3133(23) g B ∗ Bπ M η phys s [GeV] 0.6858(40) 0.6858(40)∆ B s K [GeV] 0.3127(10) 0.3126(10)∆ B s K + [GeV] -0.04157(42) -0.04157(42)∆ B s η s [GeV] 0.4000(10) 0.4000(10)∆ B s η s + [GeV] 0.0487(22) 0.0487(22) m (cid:107) m ⊥ r /a r /a r /a r /a r /a aM B s aM B s aM B s aM B s aM B s aM HISQ K aM HISQ K aM HISQ K aM HISQ K aM HISQ K aM asqtad K aM asqtad K aM asqtad K aM asqtad K aM asqtad K aM HISQ π aM HISQ π aM HISQ π aM HISQ π aM HISQ π aM asqtad π aM asqtad π aM asqtad π aM asqtad π aM asqtad π aM HISQ η s aM HISQ η s aM HISQ η s aM HISQ η s aM HISQ η s z expansion for f B s K , + and f B s η s , + . Fit result Fit resultGroup II Prior f B s K f B s K + f B s η s f B s η s + Group II Prior f B s K f B s K + f B s η s f B s η s + a h (0)1 a h (1)1 a h (2)1 c (0)1 h (0)2 c (1)1 h (1)2 c (2)1 h (2)2 c (0)2 h (0)3 c (1)2 h (1)3 c (2)2 h (2)3 c (0)3 h (0)4 c (1)3 h (1)4 c (2)3 h (2)4 d (0)1 l (0)1 d (1)1 l (1)1 d (2)1 l (2)1 d (0)2 l (0)2 d (1)2 l (1)2 d (2)2 l (2)2 e (0)1 l (0)3 e (1)1 l (1)3 e (2)1 l (2)3 e (0)2 l (0)4 e (1)2 l (1)4 e (2)2 l (2)4 B s → η s form factors on each ensemble and for eachsimulated momentum. Ensemble f B s η s (000) f B s η s (100) f B s η s (110) f B s η s (111)C1 0.8135(17) 0.7352(22) 0.6813(19) 0.6381(21)C2 0.8205(21) 0.7127(33) 0.6475(39) 0.5921(70)C3 0.8140(26) 0.7095(32) 0.6504(31) 0.6069(39)F1 0.8179(20) 0.7107(23) 0.6410(26) 0.5862(47)F2 0.8229(24) 0.7096(31) 0.6383(33) 0.5874(51)Ensemble f B s η s + (100) f B s η s + (110) f B s η s + (111)C1 1.843(10) 1.5476(62) 1.3400(63)C2 1.742(13) 1.3885(99) 1.150(17)C3 1.6802(95) 1.3855(84) 1.1771(85)F1 1.6928(71) 1.3497(55) 1.134(10)F2 1.7012(97) 1.3588(72) 1.155(11) zf + B s η s f B s η s C1C2C3 zf + B s η s f B s η s F1F2 q [GeV ] f B s η s f + B s η s data 0 0.5 1 1.5 2 2.5 0 4 8 12 16 20 FIG. 17: (color online). B s → η s form factor resultsfrom a simultaneous HPChPT z expansion are shown( top ) relative to coarse ensemble data (C1, C2, and C3),( middle ) relative to fine ensemble data (F1 and F2),and ( bottom ) in the continuum limit with physicalmasses, extrapolated over the full kinematic range.2
82 84 86 88 90 92 94 96 1 2 3 4 5 6 7 8 9 10 χ [ d . o . f . ] f , + B s η s ( ) FIG. 18: (color online). The stability of the HPChPT z expansion is demonstrated by studying the fit resultsunder various modifications, discussed in Sec. V of thetext. f r e l a ti v e e rr o r c o m pon e n t s [ % ] q [GeV ] experimentkinematicchiraldiscretizationstatisticalmatchingtotal 0 3 6 9 12 15 18 0 4 8 12 16 20 f + r e l a ti v e e rr o r c o m pon e n t s [ % ] q [GeV ] experimentkinematicchiraldiscretizationstatisticalmatchingtotal 0 3 6 9 12 15 18 0 4 8 12 16 20 FIG. 19: (color online). B s → η s ( top ) f and ( bottom ) f + relative error components. The total error (solidline) is the sum in quadrature of the components.3 q [GeV ] R || q [GeV ] R ⊥ FIG. 20: (color online). Ratios of B s → K to B s → η s form factors R (cid:107) , ⊥ as functions of q . R || r e l a ti v e e rr o r c o m pon e n t s [ % ] q [GeV ] experimentkinematicchiraldiscretizationstatisticalmatchingtotal 0 5 10 15 20 25 0 5 10 15 20 R ⊥ r e l a ti v e e rr o r c o m pon e n t s [ % ] q [GeV ] experimentkinematicchiraldiscretizationstatisticalmatchingtotal 0 5 10 15 20 25 30 35 40 0 5 10 15 20 FIG. 21: (color online). Relative error components for( top ) R (cid:107) and ( bottom ) R ⊥ as a function of q . The totalerror is the sum in quadrature of the components.4 [1] K. Hornbostel, G.P. Lepage, C.T.H. Davies, R.J. Dow-dall, H. Na, and J. Shigemitsu (HPQCD), Phys. Rev. D , 031504 (2012) [arXiv:1111.1363 [hep-lat]][2] H. Na, C. T. H. Davies, E. Follana, G. P. Lepage, andJ. Shigemitsu (HPQCD), Phys. Rev. D , 114506 (2010)[arXiv:1008.4562 [hep-lat]][3] H. Na, C. T. H. Davies, E. Follana, J. Koponen,G. P. Lepage, and J. Shigemitsu (HPQCD), Phys. Rev.D , 114505 (2011) [arXiv:1109.1501 [hep-lat]][4] J. Bijnens and I. Jemos, Nucl. Phys. B , 54 (2010);Erratum-ibid. B , 182 (2011) [arXiv:1006.1197 [hep-ph]][5] J. Bijnens and I. Jemos, Nucl. Phys. B , 145 (2011)[arXiv:1011.6531 [hep-ph]][6] A. Bazavov, C. Bernard, C. DeTar, S. Gottlieb,U. M. Heller, J. E. Hetrick, J. Laiho, L. Levkova,P. B. Mackenzie, M. B. Oktay, R. Sugar, D. Toussaint,and R. S. Van de Water (MILC), Rev. Mod. Phys. ,1349 (2010) [arXiv:0903.3598 [hep-lat]][7] G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea, andK. Hornbostel (HPQCD), Phys. Rev. D , 4052 (1992)[arXiv:hep-lat/9205007][8] H. Na, C. J. Monahan, C. T. H. Davies, R. Horgan,G. P. Lepage and J. Shigemitsu (HPQCD), Phys. Rev.D , 034506 (2012) [arXiv:1202.4914 [hep-lat]][9] E. Follana, Q. Mason, C. Davies, K. Hornbostel,G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong(HPQCD), Phys. Rev. D , 054502 (2007) [arXiv:hep-lat/0610092][10] C. M. Bouchard, G. P. Lepage, C. Monahan, H. Na,and J. Shigemitsu (HPQCD), Phys. Rev. D ,054509 (2013); Erratum-ibid. D , 079901 (2013)[arXiv:1306.2384 [hep-lat]][11] C. Monahan, J. Shigemitsu, and R. Horgan (HPQCD),Phys. Rev. D , 034017 (2013) [arXiv:1211.6966 [hep-lat]][12] E. Gulez, A. Gray, M. Wingate, C. T. H. Davies,G. P. Lepage, and J. Shigemitsu (HPQCD), Phys. Rev.D , 074502 (2006); Erratum-ibid D , 119906 (2007)[arXiv:hep-lat/0601021][13] E. B. Gregory, C. T. H. Davies, I. D. Kendall, J. Ko-ponen, K. Wong, E. Follana, E. G´amiz, G. P. Lepage,E. H. M¨uller, H. Na, and J. Shigemitsu (HPQCD), Phys.Rev. D , 014506 (2011) [arXiv:1010.3848 [hep-lat]][14] G. P. Lepage, B. Clark, C. T. H. Davies, K. Horn-bostel, P. B. Mackenzie, C. Morningstar and H. Trot-tier, Nucl. Phys. Proc. Suppl. , 12 (2002) [arXiv:hep-lat/0110175][15] C. G. Boyd, B. Grinstein, and R. F. Lebed, Phys. Rev.Lett. , 4603 (1995) [arXiv:hep-ph/9412324][16] C. M. Arnesen, B. Grinstein, I. Z. Rothstein, andI. W. Stewart, Phys. Rev. Lett. , 071802 (2005)[arXiv:hep-ph/0504209][17] C. Bourrely, I. Caprini, and L. Lellouch, Phys. Rev. D , 013008 (2009); Erratum-ibid. D , 099902 (2010)[arXiv:0807.2722 [hep-ph]][18] G. Colangelo, M. Procura, L. Rothen, R. Stucki, andJ. Tarrus JHEP (2012) 081 [arXiv:1208.0498 [hep-ph]][19] C. Bernard (MILC), Phys. Rev. D , 054031 (2002)[arXiv:hep-lat/0111051][20] C. T. H. Davies, C. McNeile, E. Follana, G. P. Lepage,H. Na, and J. Shigemitsu (HPQCD), Phys. Rev. D , 114025(2012) [arXiv:1207.0265 [hep-ph]][22] R. N. Faustov and V. O. Galkin, Phys. Rev. D , 094028(2013) [arXiv:1304.3255 [hep-ph]][23] J. Beringer et al. (Particle Data Group), Phys. Rev. D , 010001 (2012) and 2013 partial update for the 2014edition [ http://pdg.lbl.gov ][24] S. Aoki et al. (FLAG Working Group), arXiv:1310.8555[hep-lat] [ http://itpwiki.unibe.ch/flag ]; J. A. Bai-ley et al. (FNAL Lattice and MILC), Phys. Rev.D , 054507 (2009) [arXiv:0811.3640]; E. Gulez,A. Gray, M. Wingate, C. T. H. Davies, G. P. Lep-age, and J. Shigemitsu (HPQCD), Phys. Rev. D ,074502 (2006); Erratum-ibid. D , 119906 (2007) [hep-lat/0601021][25] Y. S. Tsai, Nucl. Phys. B, Proc. Suppl. 55, 293 (1997)[26] C.-H. Chen and C.-Q. Geng, JHEP (2006) 053[arXiv:hep-ph/0608166][27] Ulf-G. Meißner and W. Wang, JHEP (2014) 107[arXiv:1311.5420 [hep-ph]][28] R.C. Verma, J. Phys. G , 025005 (2012)[arXiv:1103.2973 [hep-ph]][29] The fitting software used in this paper is available on-line: see G. P. Lepage (2012). lsqfit v4.8.5.1. ZENODO.10.5281/zenodo.10236 for a general package for nonlinearleast-squares fitting, and G. P. Lepage (2012). corrfit-ter v3.7.1. ZENODO. 10.5281/zenodo.10237 for a gen-eral package for fitting 2-point and 3-point correlators.This software implements the strategies discussed in Ap-pendix A of this paper.[30] C. T. H. Davies, E. Follana, I. D. Kendall, G. P. Lep-age, and C. McNeile (HPQCD), Phys. Rev. D , 034506(2010) [arXiv:0910.1229 [hep-lat]][31] H. Ohki, H. Matsufuru, and T. Onogi, Phys. Rev. D , 094509 (2008) [arXiv:0802.1563 [hep-lat]]; W. Det-mold, C.-J. David Lin, and S. Meinel, Phys. Rev. D ,114508 (2012) [arXiv:1203.3378 [hep-lat]]; F. Bernardoni,J. Bulava, M. Donnellan, and R. Sommer (ALPHA),arXiv:1404.6951 [hep-lat];[32] C. Aubin, C. Bernard, C. DeTar, J. Osborn, Steven Got-tlieb, E. B. Gregory, D. Toussaint, U. M. Heller,J. E. Hetrick, and R. Sugar (MILC), Phys. Rev. D70