Bottomonium precision tests from full lattice QCD: Hyperfine splitting, Υ leptonic width, and b quark contribution to e+e-→hadrons

We calculate the mass difference between the Υ and ηb and the Υ leptonic width from lattice QCD using the highly improved staggered quark formalism for the b quark and including u, d, s and c quarks in the sea. We have results for lattices with lattice spacing as low as 0.03 fm and multiple heavy quark masses, enabling us to map out the heavy quark mass dependence and determine values at the b quark mass. Our results are MΥ −Mηb 1⁄4 57.5ð2.3Þð1.0Þ MeV (where the second uncertainty comes from neglect of quark-line disconnected correlation functions) and decay constants, fηb 1⁄4 724ð12Þ MeV and fΥ 1⁄4 677.2ð9.7Þ MeV, giving ΓðΥ → eþe−Þ 1⁄4 1.292ð37Þð3Þ keV. The hyperfine splitting and leptonic width are both in good agreement with experiment, and provide the most accurate lattice QCD results to date for these quantities by some margin. At the same time results for the time moments of the vector-vector correlation function can be compared to values for the b quark contribution to σðeþe− → hadronsÞ determined from experiment. Moments 4–10 provide a 2% test of QCD and yield a b quark contribution to the anomalous magnetic moment of the muon of 0.300ð15Þ × 10−10. Our results, covering a range of heavy quark masses, may also be useful to constrain QCD-like composite theories for beyond the Standard Model physics.


I. INTRODUCTION
Weak decay matrix elements calculated in lattice QCD are critical to the flavor physics program of overdetermining the Cabibbo-Kobayashi-Maskawa (CKM) matrix to find signs of new physics. For that program the weak decays of b quarks are particularly important since they give access to the least well-known CKM elements, V ub and V cb . These CKM matrix elements can be determined using either exclusive decay channels and lattice QCD form factors or inclusive decay channels and measured spectral shape functions. There continues to be some tension between the exclusive and inclusive determinations [1] that needs further improvements to both approaches to resolve. On the lattice QCD side this means developing improved approaches to B meson weak decay matrix elements, such as [2,3], but also providing more stringent tests of lattice QCD results in b physics to make sure that sources of systematic error are under full control. Here we provide three such tests using bottomonium correlation functions: the ground-state hyperfine splitting (the mass difference between the ϒ and η b ), the ϒ leptonic width and the b quark contribution to Rðe þ e − → hadronsÞ. The last two, being electromagnetic processes, can be compared with experiment free from CKM uncertainties. We obtain the most accurate results to date for these quantities and are able to include the effect of the b quark's electric charge in the calculation for the first time.
We use the highly improved staggered quark (HISQ) discretization of the quark action for these calculations. The HISQ action was developed in [4] to have small discretization errors with the leading errors, quadratic in the lattice spacing, removed. This makes the action particularly good for heavier quarks when discretization errors appear as powers of the quark mass in lattice units, which can be relatively large. This action enabled the first accurate lattice calculations in charm physics [5][6][7][8]. More recently it has been used to achieve sub-1% accuracy in the charmonium hyperfine splitting, J=ψ leptonic width, m c and c quark vacuum polarization contribution to the anomalous magnetic moment of the muon [9]. The calculation used a range of lattice spacing values from 0.15 to 0.03 fm with u, d, s and c quarks in the sea and included the effect of the c quark's electric charge [9].
Here we will extend the latter calculation to bottomonium. Because the b quark mass is much larger than that of c, we need fine lattices to reach the b with a quark mass in lattice units, am b < 1 and controlled discretization errors. Our strategy, known as the heavy-HISQ approach [10], is to perform calculations for a range of masses between c and b on lattices with a range of fine lattice spacings. We can then map out the dependence on the heavy quark mass of both the quantity being calculated and its discretization errors. This enables us to determine a physical result at the b quark mass. This approach has been very successful for decay constants and spectroscopy for heavy-light (B, B s and B c ) mesons [10][11][12] and is now being used for the form factors for B meson weak decays [2,3]. Here we will apply this approach to the ϒ for the first time.
There are alternative nonrelativistic approaches that can be used at the b quark mass on coarser lattices; see [13] for the determination of the ϒ and ϒ 0 leptonic widths using lattice nonrelativistic QCD. The nonrelativistic expansion of the Hamiltonian and the currents that appear in matrix elements gives systematic uncertainties from missing higher-order relativistic corrections and from renormalization of the lattice current to match the continuum current. These uncertainties hinder tests at high precision.
In contrast the HISQ action is relativistic and the HISQ vector current can be matched accurately and fully nonperturbatively to that in continuum QCD [14]. As we will show below, this enables us to improve the lattice QCD accuracy on the bottomonium hyperfine splitting to better than 5% and to achieve percent-level precision on the ϒ and η b decay constants and on moments that parametrize the b quark contribution to Rðe þ e − → hadronsÞ.
Our results cover the range of heavy quark masses from c to b and we will give results for decay constant to mass ratios over this range. These could be useful both for tuning of phenomenological models of QCD and as constraints on QCD-like composite theories for beyond the Standard Model physics.
The paper is organized as follows. In the next section we give details of the lattice calculation we perform. This includes a general description of the fits that we use to determine the heavy mass dependence of the quantities calculated. We then present results for the bottomonium hyperfine splitting in Sec. III, the decay constants for both the η b and ϒ in Sec. IV and the time moments of the vector current-current correlators in Sec. V. Each section includes a description of the calculation and then a discussion subsection with comparison to experiment and previous lattice QCD results. Section IV on decay constants also includes plots of decay constant to mass ratios and the ratio of vector to pseudoscalar decay constants over the quark mass range from c to b. We then give our conclusions in Sec. VI.

II. LATTICE CALCULATION
We use ensembles of lattice gluon field configurations provided by the MILC Collaboration [15] at values of the lattice spacing, a ≈ 0.09, 0.06, 0.045 and 0.03 fm. The configurations are generated with an α s a 2 -improved discretization of the gluon action [16] and include the effect of u, d, s and c quarks in the sea with the HISQ formalism [4]. The u and d masses are taken to be the same and we denote this mass m l . For most of the ensembles we have unphysically heavy u=d quarks with m l =m s ≈ 0.2 but we employ two ensembles with physical values of m l and lattice spacing values ∼0.09 and ∼0.06 fm. We expect sea quark mass effects to be small for the ϒ because it has no valence light quarks. However, an analysis of such effects is needed for accurate results. Table I lists the parameters of the ensembles that we use. The lattice spacing is determined in terms of the Wilson flow parameter w 0 [17]. On these ensembles we calculate quark propagators from random wall sources using the HISQ action and with a variety of masses, m h , from that of the c quark upward. The valence heavy quark masses that we use on each ensemble are listed in Table II. The value of ϵ Naik used in the coefficient of the Naik term in the HISQ action [4] is taken as the tree-level function of the quark mass given in [6].
We combine the quark propagators into (connected) meson correlation functions for both pseudoscalar (η h ) and vector (ϕ h ) mesons, using the local γ 5 and γ i operators converted to appropriate form for staggered quarks [8,9]. Note that we do not include quark-line disconnected correlation functions that take account of the heavy quark/antiquark annihilation to gluons. We expect the effect of the disconnected correlation functions to be very small in the heavyonium system. In [9] our result for the mass difference between J=ψ and η c mesons was accurate enough, for the first time, to see a difference with experiment of 7.3(1.2) MeV. We concluded that this was the effect of the missing disconnected correlation function on the η c mass. Here we will test for a similar effect on the η b .
On the coarsest two ensembles we use 16 time sources on each gluon field configuration for high statistics; we take eight time sources on the other ensembles. We use at least 100 configurations on each ensemble for a good statistical sample. In generating the very fine lattice (set 6 in Table I) a slow evolution in Monte Carlo time of the topological charge was observed [15], so that the ensemble does not explore many topological sectors. However, it has been shown that the impact of this on calculations for heavy mesons is negligible [20].  [15] used here with HISQ sea quark masses, m sea l (l ¼ u=d), m sea s and m sea c given in lattice units. The lattice spacing is given in units of w 0 [17]; the physical value of w 0 was determined to be 0.1715(9) fm from f π [18]. Sets 1 and 2 are "fine" (a ≈ 0.09 fm), sets 3 and 4 are "superfine" (a ≈ 0.06 fm), set 5 "ultrafine" (a ≈ 0.045 fm) and set 6 "exafine" (a ≈ 0.03 fm). The final two columns give the extent of the lattice in each spatial direction (L s ) and time (L t  We fit the correlation functions from each ensemble using a multiexponential constrained fit [21] and following the method in [9]. The fit form used for the pseudoscalar correlators as a function of t, the time separation between source and sink, is and the vector fit form is Here The term that oscillates in time in the vector case results from the use of staggered quarks. E 0 is the mass of the lowest lying state (either pseudoscalar or vector) and A 0 is related to the meson decay constants. The ground-state pseudoscalar meson we will denote as η h and the vector as ϕ h . We fit the correlation functions for all masses on a given ensemble simultaneously (with two exceptions, see Table II). This means that the correlations between results for different masses are carried through the rest of the calculation. The correlations between the ϕ h and η h correlators are safely neglected as the uncertainty in the ϕ h results dominates that for the η h . Results for the groundstate mesons are listed in Table II. We also have a limited amount of data which include the effects of quenched QED (electrically charged valence quarks, but not sea quarks). This allows us to assess the impact of QED and appropriately account for it in our error budgets. As in [9] we use photon fields in Feynman gauge in the QED L formalism [22]. Our quenched QED calculations [9,19] used a valence quark electric charge of 2=3e (i.e., the charge on a c quark), where e is the magnitude of the charge on an electron. We can use these results to determine the electromagnetic correction for the b electric charge of −ð1=3Þe. Given the smallness of α QED we take QED corrections to be linear in the quark charge squared, Q 2 , and simply rescale the effect of QED by a factor of 1=4 from that for Q ¼ ð2=3Þe. Results are given in Table III in the form of the ratio, R 0 , of results in QCD þ QED to those in pure QCD at a fixed value of the valence quark mass in lattice units. We see there that the impact of QED is tiny but visible. We showed in [9] that QED finite-volume effects were negligible for the electrically neutral charmonium mesons. This will continue to be true for the heavier mesons that we study here and so we ignore such effects.
As we showed in [9], fixing the lattice spacing from w 0 and f π , as we have done, means that QED corrections to the lattice spacing should be at the sub-0.1% level (coming from QED effects in the quark sea). We can therefore compare QCD plus quenched QED to pure QCD using the same value of the lattice spacing (i.e., that from Table I). QED affects the tuning of the lattice quark masses, however. We use the simple and natural scheme of tuning the b quark mass in both the QCD þ QED and pure QCD cases so that the ϒ mass has the physical value. We can then use our results to determine the impact of QED on the quantities we study, taking the renormalization of the quark mass into account. We will give that information, after fitting, as the renormalized ratio R, which is the ratio of the QCD þ QED result to that in pure QCD when both theories have a b quark mass separately tuned so that the ϒ mass takes the experimental value in both cases.
For each quantity that we study we must fit our results, in physical units, as a function of heavy quark mass and lattice spacing to determine results in the continuum limit at the physical b quark mass. We will use the ϕ h mass as a physical proxy for the heavy quark mass and then the physical point is defined by the ϕ h mass becoming equal to that of the ϒ.
In previous studies of heavy meson masses and decay constants using the heavy-HISQ method the HPQCD Collaboration have used fit forms that capture the heavy mass dependence as a polynomial in the inverse H s or η h mass [10,11]. In the case of heavy-light mesons this form is justified by the heavy quark effective theory (HQET) expansion. In the case of heavyonium the HQET expansion is not valid but the same form may still be expected to work as a Taylor expansion over a finite region in m h . Here, however, we choose to use a fit form that is more agnostic with regard to the dependence on the heavy quark mass. We achieve this by using cubic splines between specified knot positions. 1 We do not expect to need many knots because the quantities we study here should be smooth monotonic  Table II presented as the ratio, R 0 , of the value in QCD þ QED to that in pure QCD at fixed valence quark mass in lattice units.  1 We use splines that are monotonic between knots (Steffen splines [23]). functions of M ϕ h in the continuum limit at physical sea quark masses. The fit function for our lattice results also needs to include dependence on the lattice spacing and the mistuning of sea quark masses. Both of these effects can also depend on the heavy quark mass (M ϕ h ) through smooth monotonic functions and so we also include cubic splines in their description.
We use fits of the following form for the pure QCD part of our fit: Fða; M ϕ h Þ are the lattice QCD results in physical units of GeV for the hyperfine splitting and decay constants and GeV −1 for the time moments. A is then a dimensionful number of a size commensurate with the size of the quantity being fitted (this will be given in each section) so that the rest of the fit in square brackets is dimensionless. F 0 is a dimensionless function of M ϕ h that differs between the different quantities we examine (and will also be given in each section). It is a very simple function of M ϕ h that captures the general trend in the mass dependence, on top of which the corrections modeled by splines are relatively minor. G n denotes a cubic spline. All of the splines have different parameters but keep the same positions for the knots. The splines are functions of M ϕ h , where M ϕ h is in GeV, except for the first spline, G 0 , which provides the physical corrections to F 0 . We found that a spline in 1=M ϕ h rather than M ϕ h gave a better χ 2 for quantities like the hyperfine splitting which fall as M ϕ h grows, approximately as the inverse. We allow for two kinds of discretization effects, those that are set by the heavy quark mass am h and those that are independent of the heavy quark mass and instead are set by a fixed scale that we call Λ. We take Λ to be 0.5 GeV. The G ðiÞ 1 splines allow for heavy mass dependence in the am h discretization effects and the G ðjÞ 2 splines allow for heavy mass dependence in the aΛ discretization effects. We also include a mixed term in ðam h Þ 2 ðaΛÞ 2 with spline G 3 , although this has little impact on the fits.
The last line of the fit form in Eq. (4) allow for sea quark mass mistunings. The light quark mass mistuning parameter δ l is defined as m phys s and m phys l are taken to be the same as those used in [9]. The charm quark mass mistuning parameter is defined similarly: We allow for discretization effects and heavy mass dependence set by splines G ðkÞ 4 within the light sea quark mass mistuning term. We allow for heavy mass dependence in the sea charm mass mistuning term through the spline G 5 .
For data that includes quenched QED effects we add extra terms to the fit function in Eq. (4). The full fit takes the form For all the fits considered here we use knots placed at f2.5; 4.9; 10g GeV, taking knots at the beginning and end of the fit range and one in the middle. This is the optimal number of knots according to the Bayes factor. We have checked that varying the sums over discretization effects so that the number of terms included changes by AE1 has no significant effect. Because statistical errors are so small here (see Table II) we employ a singular value decomposition (SVD) cut [24] in the fit to account for tiny effects that are too small to be modeled by our fit form of Eq. (4). The SVD cut value times the maximum eigenvalue of the covariance matrix sets a minimum value for the eigenvalues. Any eigenvalues of the covariance matrix below this minimum are then replaced with the minimum value. This is a conservative move which increases our uncertainties, as can be seen from the impact of the SVD cut in our error budgets.
The prior information given to the fit is central values and widths for the values of the coefficients, c F , in F 0 and for the values of the spline functions at each knot. We use priors of 0 AE 1 for all of these. We use the lsqfit PYTHON module [25] to do the fits, implementing the splines with the GVAR module [26].
To obtain our final results, the fit function is evaluated at lattice spacing equal to zero, sea quark masses tuned to their physical values and M ϕ h equal to M ϒ . This is taken from experiment as 9.4603 GeV [1] with negligible uncertainty.

III. HYPERFINE SPLITTING
The hyperfine splitting, ΔM hyp , is the difference in mass between the ground-state ϕ h and η h mesons. The values for the hyperfine splitting on each ensemble for a variety of heavy quark masses are given in lattice units in column 5 of Table II. The separate η h and ϕ h masses are given in lattice units in columns 3 and 4. The impact of quenched QED at fixed valence quark mass is given in Table III. The effect of QED is similar to that for charmonium [9] but reduced because of the smaller electric charge of the b quark. The direct effect of QED on the hyperfine splitting is to increase it, through a QED hyperfine effect which has the same sign as the QCD hyperfine effect. QED also increases the meson masses, however, and this requires a retuning of the bare quark masses downward to match the same meson mass. This then has an indirect effect, increasing the hyperfine splitting by a very small amount.
Our lattice results are plotted in Fig. 1 as a function of the vector heavyonium mass, M ϕ h . The points include both pure QCD and QCD þ QED values but the QCD þ QED values are indistinguishable from pure QCD on this scale.
To fit our results for the hyperfine splitting using Eq. (4) we take A ¼ 0.1 GeV and the simple form for F 0 , F have prior values 0(1). Note that we multiply the QED correction term in the fit [Eq. (7)] by a factor of 2 because of the size of the QED corrections that we see in the results (Table III and [19]). That these prior widths are very conservative can be judged from the values and variation across Fig. 1. Evaluating the fit result at zero lattice spacing, tuned quark masses and with M ϕ h equal to the ϒ mass, we obtain the physical result for the bottomonium hyperfine splitting using connected correlation functions of This is the QCD þ QED value. For the ratio of the QCD þ QED value to the pure QCD result, we obtain Note that this is the "renormalized" ratio with the b quark mass tuned from the ϒ in both QCD þ QED and QCD. We see no significant impact of quenched QED at the 0.2% level. The fit has a χ 2 =dof of 0.73 using a SVD cut of 5 × 10 −3 . Figure 1 shows our fit curve as a function of M ϕ h in the continuum limit for tuned sea quark masses. This gives useful physical insight into how the hyperfine splitting falls as the quark mass increases. At the high mass end of the plot we mark with a black cross the experimental average value [1] for the bottomonium system. We will discuss the comparison to experiment further in Sec. III A. We note that at the lower mass end of the curve we have results for charmonium. Our fit here does not include all of the charmonium results that went into [9] but gives a value for the charmonium hyperfine splitting that is consistent (within 1σ) with [9] for the pure QCD case. The QED þ QCD result here is too small at the charmonium end of the fit curve because the QED is being included with quark charge 1=3e rather than the correct charm quark charge of 2=3e.
We will discuss in Sec. III A what the impact of quarkline disconnected (but gluon-connected) correlation functions could be on the bottomonium hyperfine splitting. For our charmonium calculation of [9] we included an estimate of the QED quark-line disconnected contribution to the hyperfine splitting coming from cc annihilation to a single photon, which then converts back to cc. The contribution of this to the charmonium hyperfine splitting is 0.7 MeV, which was a little more than half the uncertainty in our result. The equivalent contribution for the ϒ here is much smaller, at 0.17 MeV, because of the smaller electric charge of the b quark. At a size of one tenth of the uncertainty in our result in Eq. (8), this would then have negligible impact and we do not include it.
A complete error budget for the bottomonium hyperfine splitting is given in Table IV. Statistical uncertainties are divided between those arising from our two-point fits and those coming from the lattice spacing determination, both correlated between ensembles (w 0 ) and uncorrelated (w 0 =a). The uncertainty from the two-point fits is further divided in two. As already mentioned, the use of a SVD cut is conservative and increases the uncertainty in the fit output. We can calculate the contribution to an error budget  Tables II and III, with different symbols denoting different ensembles as in the legend. The errors are dominated by uncertainties from a that are correlated between the points. QCD þ QED points are shown in cyan. They are not distinguishable from their pure QCD counterparts, but are visible by being shown on top of these values. Cubic splines are used to fit the heavy mass dependence, as described in the text. The fit evaluated at the physical point and zero lattice spacing is given by the purple dashed line with error band. The experimental average value for the hyperfine splitting [1] is plotted as the black cross at the physical ϒ mass.
of the data both with and without the SVD cut applied to its correlation matrix. In the error budgets of Table IV we give the contribution from the data with the original correlation matrix under the heading "statistics." The additional contribution from the SVD cut is then defined as the square root of the difference of the squared contributions from the data with and without a SVD cut applied. The contributions from various parts of the heavy mass dependence in Eqs. (4) and (7) are shown individually, labeled by the set of spline functions for that contribution.
The fit parameters required to reproduce the physical curve of the hyperfine splitting as a function of M ϕ h plotted in Fig. 1 are given in Table IX of the Appendix. A. Discussion: Hyperfine splitting Our bottomonium hyperfine splitting result of Eq. (8) is compared to earlier lattice QCD results in Fig. 2, going back to the first lattice QCD calculation to include sea quarks [27]. Clearly the use of the heavy-HISQ approach has allowed us to reduce the uncertainty significantly (by a factor of 3) relative to these earlier results. The earlier results all use nonrelativistic actions, or actions with nonrelativistic input such as the Fermilab formalism [28], for the b quarks. This leads to uncertainties from the normalization of relativistic corrections to the Hamiltonian, such as the σ · B term that is responsible for the hyperfine splitting. We avoid this uncertainty with the HISQ action at the cost of having to calculate at multiple heavy quark masses rather than directly at the b quark mass.
As discussed in Sec. II, we have only computed connected correlators. This is also true for the earlier results except for that in [32]. This means that we are neglecting the contribution to the η b mass from its annihilation to gluons. This contribution can be related to the η b hadronic width using nonrelativistic QCD (NRQCD) perturbation theory [4]: Using the total width of the η b of 10(5) MeV [1] gives a shift to the η b mass from the leading-order term of −1.0ð5Þ MeV. This would result in an upward shift in the hyperfine splitting of approximately 1 MeV, which amounts to 0.5σ for our result [Eq. (8)]. We recently showed, for the first time, that this leadingorder analysis fails in the case of the charmonium hyperfine splitting [9] where, with the improved accuracy we were able to achieve, it becomes clear that the lattice QCD þ QED result is significantly higher than the experimental average. Assuming that this difference is the result of the   (8) is given by the top purple hexagon. Previous results (green squares) come from HPQCD/UKQCD using Oðv 4 Þ NRQCD b quarks and 2 þ 1 flavors of asqtad sea quarks [27]; the Fermilab Lattice/ MILC collaborations using the Fermilab formalism for the b quark and 2 þ 1 flavors of asqtad sea quarks [29]; S.Meinel using NRQCD b quarks with Oðv 6 Þ spin-dependent terms and 2 þ 1 flavors of domain-wall sea quarks [30]; the RBC/UKQCD collaborations using the relativistic heavy quark formalism for the b quark and 2 þ 1 flavors of domain-wall sea quarks [31] and HPQCD using radiatively improved NRQCD b quarks with Oðv 6 Þ spin-dependent terms and 2 þ 1 þ 1 flavors of HISQ sea quarks [32]. All of these results come from calculation of connected correlation functions and do not include an uncertainty from missing quark-line disconnected diagrams, except for [32]. Reference [32] includes the effect of these disconnected diagrams through the inclusion of 4-quark operators with coefficients, calculated in perturbation theory through Oðα s Þ. See the text for discussion of the impact on the hyperfine splitting through η b annihilation to gluons. The red band is the PDG experimental average [1]. The result for the hyperfine splitting calculated here shows a clear improvement on previous lattice QCD results, as well as being the first to include QED effects. This improvement is in large part due to the elimination of systematic uncertainties from the use of nonrelativistic actions which arise in previous calculations.
effect of η c annihilation missing from the lattice calculation, it seems that the leading-order perturbative analysis is misleading in this case. Presumably missing higher-order terms in the perturbative analysis or nonperturbative effects from mixing between the η c and other flavor-singlet pseudoscalar mesons [33], or both combined, have a larger effect than the leading-order term and opposite sign. In the bottomonium case the η b is considerably farther from these lighter states and so we may expect a much smaller effect from this. We also expect perturbation theory to be more reliable at the higher energy associated with bottomonium states. We therefore allow an additional 1 MeV uncertainty for the impact of η b annihilation on the hyperfine splitting and give a final result of The first uncertainty is from the lattice calculation and the second from missing quark-line disconnected contributions.
The experimental average value for the bottomonium hyperfine splitting (62.3 AE 3.2 MeV) [1] is shown by a red band on Fig. 2. A more detailed comparison with experimental results is given in Fig. 3. This makes clear the spread in the experimental results, handled in [1] by increasing the uncertainty in the average by a factor of 1.8. In particular it shows that the most recent and most precise result from BELLE [34] is noticeably lower than the others. This BELLE result is in agreement with our determination to within 1σ.
Our result is also in agreement with the PDG average to within 1.5σ. We see no disagreement with the experimental result that would signal a larger contribution from η b annihilation than the 1 MeV that we have allowed above. Indeed a shift upward of our hyperfine splitting result by 1 MeV, as suggested by leading-order perturbation theory, would improve the agreement between lattice QCD and experiment, although the shift would not be significant. In contrast, a shift downward of the bottomonium hyperfine splitting by several MeV, as we found for the charmonium hyperfine splitting, would cause tension with the experimental results.
Finally we note that the high precision we are able to achieve for the bottomonium hyperfine splitting is the result of concentrating on the ground-state mesons with a highly improved relativistic action. For a more complete picture of the bottomonium spectrum, obtained on an anisotropic lattice with the Fermilab heavy quark action and focusing on highly excited states see [38].

IV. ϒ AND η b DECAY CONSTANTS
We define the vector heavyonium meson (ϕ h ) decay constant from the annihilation matrix element as This means that we can determine the decay constant from our fits to the vector meson correlation functions using where A V 0 is the ground-state amplitude from a correlator fit of the form given in Eq. (2). Z V is the renormalization constant required to match the local vector current in lattice QCD to that of continuum QCD at each value of the lattice spacing. We use Z V values calculated in a nonperturbative implementation of the regularization-invariant symmetric momentum subtraction (RI-SMOM) scheme [14,39,40]. The pure QCD results for Z V for the HISQ action are given in [9,14]; we use values at scale μ ¼ 2 GeV. Note that no additional matching factor is required to reach MS from the RI-SMOM scheme and, because Z V has no anomalous dimensions, any μ dependence is purely a discretization effect [14].
The vector meson decay constant is the amplitude for annihilation of the valence quark/antiquark pair, into a photon, for example. It is related to the experimentally measurable leptonic width by where e h is the quark electric charge (1=3 for b). The α QED here is evaluated at the mass of the heavy quark and is equal to 1=132.15 [41] at the b. Note that our result here includes an uncertainty from the effect of η b annihilation missing from our lattice calculation. There is some tension between the different experimental results with our value favoring the most recent result from BELLE [34]. The result labeled CLEO is from [35], BABAR01 from [36] and BABAR02 from [37].
We also compute the decay constant of the pseudoscalar heavyonium meson, f η h . In terms of the parameters of our correlator fit, Eq. (1), this is defined as Because the partially conserved axial current relation holds for HISQ quarks the pseudoscalar decay constant is absolutely normalized and no Z factor is required to match to the continuum regularization of QCD. Since the pseudoscalar meson does not annihilate to a single particle, there is no experimental decay process that gives direct access to the decay constant. Its value is nevertheless of interest for comparison to that of the corresponding vector meson and other pseudoscalar mesons. The values of the decay constants, in lattice units, on each ensemble and for each heavy mass are given in the sixth and seventh columns of Table II. The decay constants converted to GeV units, and renormalized in the case of the vector decay constant, are plotted as a function of the ϕ h mass in Fig. 4. The decay constants increase with increasing ϕ h mass. Discretization effects are clearly visible that cause the lattice results to peel away from the physical curve upward. The same effect was seen previously for both heavy-light and heavyonium mesons [10,11].
We also show results in Fig. 4 that include the effect of quenched QED. Those results are given in Table III as the ratio of values in QCD þ QED to those in pure QCD. For the decay constant of the ϕ h these ratios do not include the impact of QED on the vector current renormalization factor, Z V . This was calculated in [14] for the case of a quark with electric charge 2e=3, again as a ratio of results in QCD þ QED to those in pure QCD. These results are given in Table IV of [14], with further results in Table X of [9]. The ratio is within 0.05% of 1, as expected for an Oðα QED Þ correction to a Z factor that is already very close to 1 for the HISQ action in pure QCD. Here we need results for an electric charge of e=3 so we determine the ratios in QCD þ QED to pure QCD in that case by taking the values from [9,14] (for μ ¼ 2 GeV) and dividing the difference from 1 by a factor of 4.
To fit our decay constant results as a function of lattice spacing and heavy quark mass we again use the fit form of Eq. (4) but we use A ¼ 0.7 GeV, as appropriate for the decay constant values, and a different form for F 0 from that used for the hyperfine splitting case. The dependence of the decay constants on the heavy mass is approximately linear and so we choose F are fit parameters with prior values of 0 AE 1. We fit the ϕ h and η h decay constants simultaneously, including the correlations between them, and take the same F 0 for both since they are so close in value. The spline functions that map out the differences from F 0 in physical heavy quark mass dependence and the dependence on the lattice spacing and sea quark masses take independent values in the two cases. The fit has a χ 2 =dof value of 0.44 using an SVD cut of 1 × 10 −4 . We again evaluate our fits at zero lattice spacing, physical sea quark masses and with M ϕ h ¼ M ϒ to obtain the physical bottomonium results.
We obtain, for the ϒ, For the η b , with Again, QED effects are not discernible within our 0.1% uncertainties. At the charmonium end of our range our results agree within uncertainties with the values we  Table I  obtained in [9], remembering that the calculation done here is for an electric charge that does not match that of the c quark. The error budget for both decay constants is given in Table IV. The fit curves evaluated at zero lattice spacing and physical sea quark masses are plotted as a function of heavy quark mass (given by M ϕ h ) in Fig. 4. The fit parameters required to reproduce these physical curves of the decay constants as a function of M ϕ h are given in Table X of the Appendix. Given that the heavy mass dependence and discretization effects in the vector and pseudoscalar decay constants are similar we can study the ratio of the two as a function of the heavy mass to high precision. Our results for the ratio are shown as a function of M ϕ h in Fig. 5. A slow downward drift of the ratio is seen with increasing M ϕ h from a value slightly above 1 for c quarks to a value slightly below 1 for b quarks.
To obtain a physical result for the ratio we again use the fit form of Eq. (4), now taking F 0 to be a constant, c F , since the ratio is relatively flat, so that the spline functions handle all of the mass dependence. We take the prior value of c F to be 1(1), i.e., with a very conservative width. Since we expect a lot of systematic effects to cancel in this ratio (and Fig. 5 shows that they do) we halve the prior widths on all of the correction terms in Eq. (4), i.e., we take prior values on the function values at the knots of 0.0(5). The fit has a χ 2 =dof of 0.22 and no SVD cut is required. Evaluating the fit function at the physical point gives and The total uncertainty in the ratio for the b is 1%, with a value clearly below 1. The fit curve evaluated at zero lattice spacing and physical sea quark masses is plotted as a function of M ϕ h in Fig. 5. The fit parameters required to reproduce this physical curve are given in Table XI of the Appendix.
A. Discussion: Decay constants Figure 6 compares our result for the ϒ decay constant, f ϒ , to that of an earlier lattice QCD calculation on a subset of the same gluon field configurations used here but using an improved NRQCD action for the b quarks [13]. Clearly, we achieve a considerably improved uncertainty over that of [13]. A large amount of the NRQCD uncertainty arises from the normalization and Oðv 2 =c 2 Þ improvement of the NRQCD vector current, where v is the nonrelativistic quark velocity. Here, since we use the HISQ action which is relativistic and we have performed the vector current renormalization to very high precision previously [14], these sources of uncertainty are effectively eliminated. Figure 6 also compares our result for f ϒ to that obtained from the experimental average for the ϒ leptonic width using Eq. (14). Using Γðϒ → e þ e − Þ ¼ 1.340ð18Þ keV [1] gives where the first uncertainty comes from the experimental uncertainty in Γ and the second allows for an Oðα QED =πÞ uncertainty from higher order in QED terms in Eq. (14) coming, for example, from final-state radiation. Note that using α QED of 1=137 here instead of 1=132.15 would increase the experimental value of f ϒ by 3.7% or 25 MeV. This is several times larger than either the experimental uncertainty or our lattice QCD uncertainty. Figure 6 shows good agreement, within 1σ, between our lattice QCD result and that from experiment [Eq. (22)]. The experimental uncertainty is about half that from our lattice QCD result.
Our result for f ϒ can be converted into a determination of the width for ϒ decay to light leptons in the Standard Model using Eq. (14). This gives where the first uncertainty comes from the lattice QCD result and the second allows for a relative Oðα QED =πÞ correction to Eq. (14) from higher-order QED effects.
Our result for f η b can be compared to an earlier HPQCD lattice QCD result using HISQ quarks and the heavy-HISQ approach on gluon field configurations including the effect of 2 þ 1 flavors of asqtad sea quarks [11]. That work obtained a value f η b ¼ 667ð6Þ MeV, which is significantly lower (by 4σ) than our result here. The discrepancy is most likely to result from a bias in the earlier results from not having values on lattices with spacings as fine as we do here. Another possible source of the discrepancy is the fact that the earlier calculation did not include c quarks in the sea. Having more flavors of quarks in the sea results in a slower running of the strong coupling constant. Hence, using the language of potential models, we expect the Coulomb-like term in the heavy quark potential [of the form −4α s ðrÞ=ð3rÞ] to have a larger value for α s at the short-distance scales to which the η b meson decay constant is sensitive. This corresponds to a deeper potential at short distances and a correspondingly larger "wave function at the origin," which is the quantity in a potential model that translates approximately into the decay constant. This effect could explain some of the discrepancy but is unlikely to be large enough to explain it all. The calculations in [11] also used a different form to fit the lattice results as a function of heavy quark mass (in that case using as proxy M η h ). This consisted of multiple powers of the inverse heavy quark mass multiplied by a leading function of the form ðM=M 0 Þ b where b was allowed to float. We have checked that using that fit form here gives us results for f η b very consistent with our spline fits, so the discrepancy with [11] is not related to the form of the fit used.
The ratio of vector to pseudoscalar decay constants as a function of heavyonium mass provides a test of our understanding of these mesons. In the language of potential models the heavyonium vector and pseudoscalar mesons differ only through spin-dependent relativistic corrections to the central potential [42]. The size of relativistic corrections falls as the heavy quark mass increases and the mean squared velocity of the heavy quarks falls. In the infinite quark mass limit pseudoscalar and vector heavyonium mesons have the same mass and the same wave function at the origin. The decay constants differ, however, by the matching factors that are needed to renormalize temporal axial and spatial vector currents from this nonrelativistic framework to full continuum QCD. The ratio of the vector to pseudoscalar heavyonium decay constants would then be expected to become the ratio of the vector to temporal axial vector matching factors in the heavy quark limit. The matching factors come from high-momentum regions of phase space and so can be calculated in QCD perturbation theory. An Oðα s Þ matching calculation was done in [43] for spin-independent nonrelativistic QCD and gave the result From this we conclude that we would expect the ratio of vector to pseudoscalar decay constants to be below 1 for large heavy quark mass. Equation (24) expects the difference from 1 to be Oð5%Þ, taking α s ≈ 0.25, but this formula will have corrections from higher orders in α s . A value for the ratio of 5% below 1 is very consistent with our results in Fig. 5, however. Very similar behavior is seen for the ratio of vector to pseudoscalar decay constants for heavy-light mesons from lattice QCD calculations. The decay constant of the D Ã s meson is found to be several percent larger than that of the FIG. 7. Upper panel: the physical ratio (evaluated at zero lattice spacing and with sea quark mass mistunings set to zero) of the mass to decay constant for the pseudoscalar heavyonium meson, η h , as a function of the ratio of pseudoscalar to vector heavyonium masses. Lower panel: the physical ratio of decay constant to mass for the vector heavyonium meson, ϕ h (the quantity denoted f V in [52]), plotted against the same ratio of masses. D s [44][45][46] whereas that of the B Ã s is a few percent below that of the B s [45,47]. This behavior can be understood on the same basis as the arguments for heavyonium above. In the heavy-light case an α 3 s calculation of the matching factors is available in the infinite heavy quark mass limit [48]. The corrections to the Oðα s Þ formula for the ratio [which is the same as for heavyonium in Eq. (24)] are sizable but have the same (negative) sign and so do not change the qualitative behavior of the difference of the ratio from 1.
Having performed the fits of the previous subsections we now have physical values for the decay constants not only at the b quark mass but also at the full range of masses between the c and b quark masses. The physical curves as a function of meson mass in Figs. 4 and 5 could be used to tune phenomenological QCD potential models, which often differ markedly on features of heavyonium physics such as details of the wave function even when reproducing the spectrum (see, for example, [49][50][51]).
They may also be useful beyond QCD. In [52] lattice QCD results across a range of masses were collected with the intention of providing useful information for phenomenologists studying strongly coupled beyond the Standard Model (BSM) theories. These theories are often QCD-like but typically with heavier (relative to the confinement scale) fundamental fermions than the light quarks of QCD. Reference [52] makes the point that information from lattice QCD calculations about how (for example) meson masses and decay constants depend on quark masses can be useful to constrain such BSM theories. This then requires lattice QCD results for quark masses not at their physical values, as we have here. The lattice QCD results need to be presented in an appropriate way with dimensionless combinations of decay constants and masses on both axes. A convenient x axis is the ratio of pseudoscalar to vector meson mass. In [52] the square of this quantity was used since the lattice QCD results were concentrated at light quark masses. Here, since we have heavy quarks and the ratio of pseudoscalar to vector meson masses is close to 1, we simply use the ratio.
Dimensionless ratios are readily obtained for our raw lattice results using the values in Table II. Correlations can be ignored because statistical uncertainties are so small. In the following we construct appropriate ratios from our fit functions in the limit of zero lattice spacing and physical sea quark masses and do not include the raw lattice results in the figures, for clarity. One useful quantity [53] is the ratio of the pseudoscalar meson mass and decay constant for a meson made of quarks of degenerate mass (i.e., the "pions" of the BSM model). Using the physical heavy mass dependence of f η h extracted from our fit we display the ratio of M η h and f η h as a function of the ratio of pseudoscalar to vector meson masses in Fig. 7. Our results show values of M η h =f η h around 10, and continuing to rise, as the ratio of pseudoscalar to vector meson masses heads toward 1. Note that our definition of the pseudoscalar decay constant in Eq. (15) corresponds to the normalization f π ≈ 130 MeV.
As discussed in [52] composite models of a dark sector in which a "dark ρ" meson couples to ordinary matter through a dark photon (see, e.g., [54]) need information on the vector meson decay constant for an appropriate range of fermion masses. The ratio of vector meson decay constant to vector meson mass is denoted f V in [52]. In our convention for the vector meson decay constant [Eq. (12)] it is f ϕ h =M ϕ h . We plot f ϕ h =M ϕ h against the pseudoscalar to vector meson mass ratio in Fig. 7. We see that this ratio becomes small as the ratio of pseudoscalar to vector meson masses heads toward 1. It also has relatively strong dependence on the mass ratio, so using an approximately constant value (based, for example, on naive dimensional analysis) would not agree well with our results.
We also note that accurate lattice QCD results are available at the ratio of pseudoscalar to vector meson masses of 0.673, which corresponds to ss mesons when only connected correlation functions are calculated. This means that the pseudoscalar meson is not allowed to annihilate and mix with flavor-singlet mesons made from lighter quarks, and likewise the vector decay to two pseudoscalar mesons incorporating lighter quarks is not included. This is then the scenario that would match that required in a composite BSM scenario. For this case HPQCD calculates M η s =f η s ¼3.801ð16Þ [18] and f ϕ s =M ϕ s ¼ 0.233ð3Þ [55]. These results must connect smoothly to the ones shown in Fig. 7 as the quark mass is increased.
The ground-state vector heavyonium decay constant is determined by the amplitude of the state that dominates the correlator at large times and this can be connected to experiment via the leptonic width, as we have seen. We can also calculate the time moments of the correlator. These depend on the behavior of the correlator at shorter time distances and can also be connected to experimental results [8,56]. The moments of the vector heavyonium correlator are defined by wheret is lattice time symmetrized around the center of the lattice, C ϕ h is the vector two-point correlation function and  Table I  Z V is the renormalization factor for the heavyonium vector current operator used. Results for ðG V n =Z 2 V Þ 1=ðn−2Þ in lattice units on each of our ensembles are given in Table V for n ¼ 4 to 10. The power 1=ðn − 2Þ is taken to reduce all the moments to the same dimension. We take the Z V factor for the vector current to be the same one used for the leptonic width above [14]. Figure 8 then shows the physical results for these moments as a function of M ϕ h .
Results that include quenched QED corrections for a subset of ensembles are given in Table VI. These are given as the ratio of the result in QCD þ QED to that in QCD at fixed valence quark mass in lattice units. The values of R 0 are very slightly below 1, as for charmonium [9]. The difference from 1 is even smaller here because of the smaller quark electric charge. Note that the vector current is not renormalized in these raw results and QED effects in Z V must also be taken into account, as for the decay constant [14]. These results are also plotted in Fig. 8 as the cyan points. The impact of QED is not visible.
To fit the time-moment results as a function of lattice spacing and heavy quark mass we again use the fit of Eq. (4), supplemented with QED effects in Eq. (7). For the time moments we use F 0 ¼ ðc ð0Þ F þ c ð1Þ F ð3 GeVÞ=M ϕ h Þ, and the dimensionful parameter A is taken as 0.5 GeV −1 for every moment. The prior values on c ð0Þ F and c ð1Þ F are taken to be 0 AE 1 for each moment. We fit all moments separately using a SVD cut of 5 × 10 −4 in all cases. The χ 2 =dof of the fits, in order of increasing n, are 0.9, 0.19, 0.26 and 0.4. The curves in Fig. 8 show the fit results evaluated at zero lattice spacing and with tuned sea quark masses. Table VII gives our results for the time moments evaluated at the b quark mass in the continuum limit, with their total uncertainties. The corresponding error budget is given in Table VIII. In the next section we compare these results to earlier lattice analyses and values determined from experimental data for Rðe þ e − → hadronsÞ. We will also use the results to improve the determination of the b quark contribution to the hadronic vacuum polarization term in the Standard Model determination of the anomalous magnetic moment of the muon. Column 3 of Table VII gives the ratio of the QCD þ QED result to that in pure QCD for each moment. Again we are not able to distinguish any impact of QED on the results at the level of our uncertainties (which range from 0.4% down to 0.1%).
A. Discussion: Time moments and a b μ In Fig. 9 we compare our results for the time moments to those of an earlier HPQCD calculation that used NRQCD b quarks [13]. For the NRQCD results the key sources of TABLE VI. Quenched QED corrections, for quark electric charge e=3, to the time moments given for a subset of the results in Table V for the (n − 2)th root of the unrenormalized G V n =Z 2 V . The results are given as the ratio, R 0 , of the value in QCD þ QED to that in pure QCD at fixed valence quark mass in lattice units.  (11) TABLE VII. Results for the time moments of the bottomonium vector current-current correlator obtained from evaluating our fit functions in the continuum limit at the b quark mass. These are given in the second column for moment numbers listed in the first column. The results extracted from experimental data in [57] are given in the third column for comparison. The fourth column gives the quenched QED correction to these moments, as a ratio of the value in QCD plus QED to that in pure QCD with a tuned b quark mass (to reproduce the ϒ mass from experiment) in both cases. All of the ratios are consistent with 1.0. error were from the vector current normalization (using a method based on matching the time moments to continuum perturbation theory) and from the lattice spacing dependence effects in the NRQCD action. Our uncertainties here are a considerable improvement (by over a factor of 2) on the NRQCD results, because we have a very accurate vector current normalization and have results over a large range of lattice spacing values to control the lattice spacing dependence. Figure 9 shows that our results agree within 2σ with the values extracted for the q 2 -derivative moments, M k (n ¼ 2k þ 2), of the b quark vacuum polarization using experimental values for R e þ e − ¼ σðe þ e − → hadronsÞ=σ pt [57]. The appropriate normalization of these results for the comparison to ours is Our results from lattice QCD have considerably larger uncertainties than those of the experimental values but together these results provide a further test of QCD at the level of 2%. We may also use these time moments to extract the b quark connected contribution to the leading-order hadronic vacuum polarization contribution to the anomalous magnetic moment of the muon. This was done in [13] and, given that we have improved on the time moments of that work, we provide an update here. We obtain This agrees with the value in [13] with an improvement in uncertainty of a factor of 2.5. Since the b quark is so heavy, this is not a significant contribution to the anomalous magnetic moment of the muon [58].

VI. CONCLUSIONS
We have used the fully relativistic HISQ action to calculate the masses and decay constants of ground-state bottomonium mesons in lattice QCD including the effects of u, d, s and c quarks in the sea. We have used very fine lattices and a range of heavy quark masses at each lattice spacing to control the discretization effects as a function of heavy quark mass along with the physical dependence on the heavy quark mass of the quantities being studied. We have used a fit function with completely generic dependence on the heavy quark mass in each of its component pieces, capturing this dependence through cubic spline functions. Values for bottomonium are obtained by evaluating the fit function at zero lattice spacing with tuned sea quark masses and a valence quark mass tuned to that of the b, defined to be the point at which the ϒ mass agrees with experiment. We have also included an analysis of the impact of the electric charge of the valence b quarks on the quantities being studied. The results given are from the QCD þ QED fit but, in all cases, we find the impact of QED to be negligible at the level of our uncertainties.
Our results yield the most precise, to date, lattice calculation of the bottomonium hyperfine splitting. We obtain the value [repeating Eq. (11)] The first uncertainty is from our fit results (see error budget in Table IV) and the second uncertainty is from an estimate of missing quark-line disconnected contributions that would affect the mass of the η b meson. Our result is in agreement with, but on the low side of, the experimental average value [1]. It tends to favor the most recent experimental result obtained by the BELLE Collaboration [34], although uncertainties (both ours and from the experiment) are still too large to draw strong conclusions from this. We also provide the most precise lattice QCD determination of the ϒ decay constant, which can be used to determine the ϒ leptonic width. Our uncertainty of 1.5% is 3 times better than the previous lattice QCD calculation of [13]. The big advantage of using a relativistic formalism, as we do here, is that the vector current can be normalized very accurately and nonperturbatively [14]. Our result [repeating Eq. (16)] is with the error budget in Γðϒ → e þ e − Þ ¼ 1.292ð37Þð3Þ keV: The first uncertainty is from our result for f ϒ and the second from possible Oðα QED =πÞ corrections to the FIG. 9. Comparison of different determinations of the four lowest time moments of the bottomonium vector current-current correlator. The three determinations are, from the top, this work, the previous calculation by the HPQCD collaboration using NRQCD b quarks [13] and the values obtained from experimental data on Rðe þ e − → hadronsÞ in [57].
formula connecting decay constant and leptonic width [Eq. (14)]. This is to be compared with the current experimental average of 1.340(18) keV [1]. We see that our result is in good agreement with experiment and our uncertainty is just twice as large.
The decay constant of the η b can also be accurately calculated with our approach. There is no experimental decay rate that can be directly compared to this determination, but the value of f η b is important for our phenomenological understanding of the relationships between decay constants for different mesons. We obtain [repeating Eq. (18)] In particular, repeating Eq. (20), we find that i.e., less than 1. This is in contrast to the charmonium case where f J=ψ =f η c is larger than 1 [9]. Figure 5 shows how the ratio of the decay constants for vector and pseudoscalar heavyonium mesons varies with heavy quark mass. This is qualitatively similar to the behavior seen for the decay constants of heavy-light mesons [47]. Finally, in Fig. (7) we plot the ratios of mass to decay constant for pseudoscalar and vector mesons as a function of the ratio of pseudoscalar to vector meson masses. These may provide useful information to constrain these ratios in QCD-like beyond the Standard Model scenarios. The low time moments of the bottomonium vector current-current correlator provide a further opportunity to compare lattice QCD results to experiment, where the matching inverse-s moments of the b quark contribution to Rðe þ e − → hadronsÞ can be determined. Our results for the fourth, sixth, eighth and tenth time moments are given in Table VII, where they can be compared to the results obtained from experiment. Our uncertainties are 2% so they provide the most accurate test to date for these quantities. The time moments can be used to determine the b quark contribution to the anomalous magnetic moment of the muon. We find [repeating Eq. (27)] Together these results demonstrate how the properties of low-lying bottomonium states can be determined in a fully relativistic calculation in lattice QCD and the gains in precision that such an approach makes possible. The results given here also allow us to improve the fully nonperturbative determination of the ratio of quark masses, m b to m c . We will present this analysis separately.

ACKNOWLEDGMENTS
We are grateful to the MILC Collaboration for the use of their gluon field configurations and for the use of MILC's QCD code. We have modified the code to generate quenched U(1) gauge fields and incorporate those into the quark propagator calculation as described here. We are grateful to B. Galloway for contributions to this project at a very early stage, and to R. Horgan, C. McNeile and J. Rosner for useful discussions. Computing was done on the Darwin supercomputer at the University of Cambridge High Performance Computing Service as part of the DiRAC facility, jointly funded by the Science and Technology Facilities Council, the Large Facilities Capital Fund of BIS and the Universities of Cambridge and Glasgow. We are grateful to the Darwin support staff for assistance. Funding for this work came from the Science and Technology Facilities Council and the National Science Foundation.

APPENDIX: RECONSTRUCTING THE HEAVY QUARK MASS DEPENDENCE
We give here the fit parameters that enable our fit results for the dependence on heavy quark mass of the hyperfine splitting, decay constants and ratio of decay constants to be reconstructed. The pieces of Eq. (4) that give the physical curves in the continuum limit are F 0 ðM ϕ h Þ and G 0 ð1=M ϕ h Þ, multiplied by dimensionful constant, A (absent for the case of the ratio of decay constants). We ignore here the QED pieces of the fit; these have a negligible effect in all cases.
F 0 is a simple function with at most two parameters, c ð0Þ F and c ð1Þ F . G 0 is a Steffen spline function [23] with three knots at 2.5, 4.9 and 10.0 GeV in M ϕ h , so that in 1=M ϕ h they are at 1=2.5, 1=4.9 and 1=10.0. Tables IX, X and XI give the mean and standard deviation of c  values at the three knots of G 0 : G k1 0 , G k2 0 and G k3 0 . This is followed underneath by the correlation matrix between these parameters. The parameters are strongly correlated and this is why we give the values to four significant figures. The splines can easily be implemented using the GVAR PYTHON module [26]. TABLE X. Fit parameters for F 0 and G 0 for the fit of Eq. (4) to the vector decay constant of the vector heavyonium ϕ h meson (upper set) and the decay constant of the pseudoscalar heavyonium η h meson (lower set), both as a function of the mass M ϕ h . The dimensionful constant, A, is 0.7 GeV in these cases and F 0 ¼ c ð0Þ F þ c ð1Þ F × M ϕ h =3 GeV. The top row of each set gives the mean and standard deviations for c ð0Þ F and c ð1Þ F and the values at the three knot positions for G 0 . The correlation matrix for these five parameters is given underneath. These results enable the red fit curves of both plots of Fig. 4 to be reconstructed within the errors given.