{'fraction_non_alphanumeric': 0.02840459435758133, 'fraction_numerical': 0.01499369024447629, 'mean_word_length': 4.245596319982049, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scaling of the Upsilon Spectrum in Lattice NRQCD

We present results for the spectrum of b-bbar bound states in the quenched approximation for three different values of the lattice spacing, in the range 0.05fm to 0.15fm. We find our results for spin-independent splittings in physical units to be independent of the lattice spacing, indicating the absence of systematic errors from discretisation effects. Spin-dependent splittings are more sensitive to the lattice spacing and higher order corrections to the action; we discuss the size of these effects and what can be done to arrive at a physical result.


I. INTRODUCTION
Accurate calculations of the hadron spectrum in lattice QCD require control of systematic errors. This has become a very important issue now that statistical errors have been reduced in recent years to the point where systematic errors can dominate the results reported.
A major source of systematic error is that arising from the use of a space-time lattice with finite lattice spacing. All operators in the continuum Lagrangian must be replaced with discrete versions and discretisation errors consequently appear. This means that physical results (for example a mass in GeV) depend on the value of the lattice spacing. This is obviously wrong. Since the lattice is simply a regulator for the theory, physical results must not depend upon its value. One approach has been to extrapolate to zero lattice spacing for a 'continuum' result. This is difficult numerically, especially if the variation with lattice spacing is severe. However, recent progress in understanding discretisation errors and how to formulate an improved action [1,2], has meant that we can obtain essentially continuum results at finite values of the lattice spacing. The lattice spacing dependence is reduced to such a low level that extrapolation is unnecessary.
The spectrum of bottomonium bound states is one of the most accurate calculations that can be done on the lattice [3]. Since the b quarks are non-relativistic in these systems (v 2 /c 2 ≈ 0.1) a non-relativistic action can be used [4,5]. This allows a b quark propagator to be calculated on one sweep through the gluon field configuration with low computational cost. Multiple sources can be used on a single configuration because the bound states are much smaller than the volume of a typical lattice. Also sources for both ground and excited states can be used, allowing multi-exponential fits to hadron correlators and improving the confidence in the fitted masses. These techniques mean that very small statistical errors can be obtained and the improvement of systematic errors becomes a priority. In this paper we discuss the issue of discretisation errors for the bottomonium spectrum [6].
The approach that we use, Non-relativistic QCD (NRQCD) [4,5], is an effective field theory. Its Lagrangian is suited to a description of non-relativistic quarks since operators are classified according to the powers of v 2 /c 2 they contain, where v is the velocity of the heavy quark. The number of operators to be included can then be truncated at a fixed order in v 2 /c 2 and this is clearly a sensible thing to do if v 2 /c 2 ≪ 1. The renormalisability of QCD is lost in this process but physical results are still obtained by putting an explicit momentum cut-off into NRQCD. This cut-off should exclude relativistic momenta and thus be of the same order, or smaller, than the heavy quark mass. On the lattice this cut-off is provided by the lattice spacing, with Ma > ∼ 1. The excluded momenta cause renormalisation of the coefficients of the NRQCD operators when, say, lattice NRQCD is matched to full continuum QCD. The coefficients will be well-behaved and essentially cut-off independent provided that the cut-off is not too large. Any attempt to take the cut-off to infinity (lattice spacing to zero) will cause them to diverge as the non-renormalisability of the theory becomes apparent. Thus, no continuum extrapolation can be done for lattice NRQCD. However, as discussed above, a continuum extrapolation is not necessary for a suitably improved action. All that is necessary is to demonstrate lattice spacing independence of physical results. For the bottomonium spectrum from NRQCD this should be possible in a region of lattice spacing M b a > ∼ 1, and this is what we show in this paper. The size of discretisation errors will vary from one quantity to another. In general it is to be expected that the coefficient of the dependence on a should represent some typical momentum scale appropriate to that quantity [7]. For the light hadron spectrum this would then be a few hundred MeV. For heavy hadrons the scale of discretisation errors is likely to be larger. The scale is not set by the heavy quark mass since this is an irrelevant scale to the dynamics of the bound states. It is set rather by typical momenta exchanged inside the hadrons. For bottomonium these momenta are of order 1 GeV and so discretisation errors might be expected to present a problem on coarse lattices if the action is not improved. Here we report results with leading order (a 2 ) discretisation errors removed from spin-independent terms, but not from the spin-dependent terms (which are of lower order in the non-relativistic expansion).
Section 2 describes the lattice calculations and results at three different values of the lattice spacing. Section 3 discusses the scaling behaviour of spin-independent and spindependent splittings. Section 4 contains our conclusions.

II. NRQCD CALCULATIONS AND RESULTS
Quark propagators in lattice NRQCD are determined, in a single pass through the gaugefield configuration, from evolution equations that specify the propagator for t > 0 in terms of its value at t = 0. We use here: [3,8] H 0 is the kinetic energy operator, the lowest order (in v 2 /c 2 ) term in the Hamiltonian: The correction terms to the Hamiltonian that we include in δH are O(v 4 /c 4 ). They comprise relativistic corrections to the spin-independent H 0 as well as the first spin-dependent terms that give rise to spin-splittings in the spectrum.
The last two terms in δH come from finite lattice spacing corrections to the lattice laplacian and the lattice time derivative respectively [5]. ∆ is the symmetric lattice derivative and ∆ (4) is a lattice version of the continuum operator D 4 i . We used the standard traceless cloverleaf operators for the chromo-electric and magnetic fields, gE and gB. The parameter n is introduced to remove instabilities in the heavy quark propagator caused by the highest momentum modes of the theory.
We tadpole-improve [9] our lattice action by dividing all the gauge fields, U, that appear in E, B, and the covariant lattice derivatives fields by u 0P , the fourth root of the plaquette. This is most easily done as the U µ 's are read by the code that evolves propagators. Tadpoleimprovement of the action allows us to work with tree-level values for the c i 's in δH (i.e. 1) without, we believe, having to worry about large renormalizations [10]. Hence our lattice action depends only on two parameters, the bare mass M 0 b and the QCD coupling constant, g. Table I shows the parameters used in the calculations at 3 different values of β. The configurations were all generated using the standard unimproved Wilson plaquette action and generously made available to us by the UKQCD collaboration [11] and by Kogut et al [12]. The results described here at β = 6.0 agree with our previous results [3] but generally have higher precision, because of an increased number of sources on different time slices and the increased length of the lattice in the time direction.
Once the quark propagators have been calculated it is straightforward to obtain antiquark propagators and meson correlation functions. We used the standard interpolating operators described in [3] with source and sink 'smearing functions'. We worked in Coulomb gauge and took wavefunctions for smearing functions, either from a Richardson potential (β = 6.0 and 6.2) or from a Coulomb potential (with modifications, β = 5.7). We took a ground state wavefunction and 2(1) radial excitations for S states at β = 6.0 and 6.2 (β = 5.7). For P states we used a ground state wavefunction and 1(0) radial excitations. We also used local sources and sinks which were delta functions for S states and combinations of delta functions for higher orbital excitations. In addition we looked at S-state mesons with small non-zero momenta. In the following discussion correlation functions at zero momentum will be denoted ( 2S+1 L J ) ab where a is the source smearing function and b the sink smearing function with l for a delta function (local operator), 1 for the ground state, 2 for the first excited state and so on.
At β = 5.7 we summed over both initial quark spins. At β = 6.0 and 6.2 we saved CPU time by fixing the initial quark spin to +1/2, since the spin-flip operators in the Hamiltonian are suppressed with respect to H 0 . We then used the strong correlations between different polarisations to obtain reduced errors on the spin splittings for P states [3].
As described in [3] we used multi-exponential fits to the multiple correlation functions obtained by different combinations of source and sink. This allowed us to obtain ground state energies and one or two excited state energies. Two different types of fit were employed; the 'matrix' fit and the 'row' fit. The matrix fit used the matrix of correlators obtained with ground and excited state sources and sinks. The row fit used the row of correlators with ground and excited state sources and local operator sinks. We found the correlators with both local source and sink to be of very little use in fitting. We had a large number of measurements in every case and so did not run into problems with our covariance matrix, even for multi-exponential fits with several parameters. The two different fits gave consistent results within the errors that we quote.
Final fitted values were chosen by monitoring the quality of the fits (Q) for given ranges of fitting time, t min to t max , as well as the stability of the fitted parameters. For a given type of fit, Q generally increases sharply with t min until it reaches a plateau. The first fit for which this happens is taken as the preferred value. In general the nth excited energy is taken as reliable from a fit to n + 1 exponentials. Table II shows the quality of our fitted results for a 2 exponential fit to 2 correlators for the 3 S 1 ; ( 3 S 1 ) 1l and ( 3 S 1 ) 2l , at each β value. Notice how Q increases from small values of t min as contamination from a third state dies away. Notice also how stable the fitted ground state energies are for a very large range of t min values. It is interesting to study how the noise in the meson correlators changes with β. We expect the ground state meson correlators ( 1 S 0 ) to have noise governed by the same mass as the signal since it is the lightest mass available [14]. This means that the errors in an effective mass plot will not grow with lattice time. Figure 1 shows this clearly for the effective masses from the ( 1 S 0 ) 1l correlator. The size of the errors at different values of β reflects partly the different statistics available for the different sets of configurations (see Table I). If we multiply the errors at β of 5.7 by √ 2 and at 6.0 by 2, for the different number of time origins (assuming these are independent), then the errors are in the ordering 5.7 > 6.0 > 6.2. The error depends on the overlap between the squared correlation function and two 1 S 0 particles [14]. On coarser lattices, the local sink will provide a better overlap with two 1 S 0 particles than on finer lattices and so we would expect the error to be larger. If instead we compare correlation functions in which the ground state smearing is applied at both source and sink, ( 1 S 0 ) 11 , then the errors at all three values of β are very similar when adjustments for statistics are made as above. This reflects the fact that the noise should not change if the physical overlap with two 1 S 0 states doesn't change.
For higher states than the ground state the noise grows exponentially with time according to the splitting between that state and the 1 S 0 . Figure 2 shows this effect for the 1 P 1 correlator with ground state smearing at the source and a local sink. We expect the doubling time for the error to be ln(2)/(1P − 1S) = 1.6 GeV −1 , and this is roughly true at all three values of β. Again the absolute size of the error at fixed physical time is very similar between all 3 β values when adjustments for the different statistics are made as above. However, on the coarser lattice many fewer lattice time points occur before the noise grows overwhelmingly large. Table III shows energies in dimensionless units obtained from our fits at each value of β. 3 exponential fits were used in general so the values for the 3S states should be used with some caution. 3 P is the spin average of the 3 P 0,1,2 states defined by This is obtained by measuring spin splittings (see below) between the 3 P states and the 1 P 1 , and the energy of the 1 P 1 . Because the quark mass term is missing from our Hamiltonian the zero of energy becomes shifted so that the energies measured in the simulation and given in Table III cannot be directly converted to hadron masses. Differences in energy can be converted directly to physical units using a value for the lattice spacing, but to obtain absolute masses we need to know the energy shift. It is sufficient to calculate an absolute mass for one meson only, and the one for which the most accurate calculation can be done and compared to experiment is the 1 3 S 1 , the Υ.
To calculate the absolute mass of the Υ we measure the dispersion relation from the energy of meson correlation functions at small, non-zero momenta, and fit to a non-relativistic energy-momentum form: E N R is the energy at zero momentum normally measured (and given in Table III). C 1 is a constant to be obtained from the fit. aE(p) − aE N R is obtained accurately by a single exponential fit to the bootstrapped ratio of correlators at finite and zero momentum. We use the lowest 1 or 2 non-zero momenta in the fit. Table III shows the kinetic masses in lattice units, aM Υ , obtained at the 3 different values for β for the bare quark masses given in Table I. The value at β = 6.0 is taken from Ref. [3] and has not been recalculated on the Kogut et al configurations. P-wave spin splittings can also be obtained most accurately from ratio fits. Single exponential fits are performed to the ratio of appropriate polarisation components to maximise the correlations, as discussed above. Table III shows the results for different splittings in lattice units at the 3 different values of β. Table IV gives a more detailed breakdown for different polarisation components of the 1P fine structure at β = 6.0 and 6.2. For a given splitting there is no significant difference between different polarisation components so we average to get a final value and allow for variations in the error. We see no significant difference between 3 P 2 T states and 3 P 2 E states. This has been checked explicitly by taking ratios of those correlators. For the hyperfine splitting (Υ − η b ) we are able to extract both ground and excited splittings (for the first time in a lattice calculation) from simultaneous fits to the 3 × 3 matrix of correlators for the 1 S 0 and 3 S 1 . These results are also given in Table III.
The wavefunction at the origin is calculated from the ratio of amplitudes of row and matrix fits as described in ref. [3].

A. Setting the scale
One of the useful features of the spectrum of heavy quark bound states is that the splittings between radial and orbital excitations, spin-averaged, are to a good approximation independent of quark mass in the region between bottom and charm. Since not all the bottomonium fine structure has been seen experimentally, this statement relies to some extent on estimates of the spin splittings that have not been measured. However, since spin splittings are very small for bottomonium systems (≈ 10% of spin-averaged splittings), we still expect little quark mass dependence for radial and orbital splittings in the region of M b when non-spin-averaged splittings are used. This allows us to set the scale from lattice calculations, independently of the requirement to tune the bare lattice quark mass to get the right kinetic mass for the Υ.
In Table V we show values for the lattice spacing, obtained by fixing various radial and orbital splittings to experiment (spin-averaging where possible), at the 3 different values of β. These lattice spacing determinations are very accurate ones for setting the scale in the determination of α s [15]. Notice that different splittings at a given value of β give slightly different values for a −1 . This is a feature of the quenched approximation which we return to below. First we describe how the determination of a −1 is done and how errors are assigned.
Since the NRQCD action of eq. 3 is corrected through O(a 2 ) for discretisation errors we would also like to remove other O(a 2 ) errors that come from using gluon fields generated with the simple plaquette action. Fortunately, these errors can be corrected for after the calculation. Perturbatively the correction appears as a shift to energies and is related to the wavefunction at the origin. It can then be written in terms of the hyperfine splitting for s states (for p states the shift is zero). We use [15] with M b set to 5 GeV. The resulting shift to the splittings of Table V is given in column 4. This shift is added to the splittings in lattice units before they are divided into the physical splitting to obtain a −1 . For the shift for 2S states we use the ratio of 2S to 1S hyperfines determined at β = 6.0 and given in Table III. The statistical error in the splitting is then inflated by a∆M g /2 before calculating the statistical error in a −1 given in column 5. This determination of a −1 also has systematic errors, as in all lattice determinations. We attempt here to quantify the errors relevant to our calculation. There are two sources, physical and unphysical. The errors from higher order relativistic corrections which have been ignored are physical and will give the same percentage error at all values of the lattice spacing. The unphysical systematic errors come from higher order discretisation corrections that have not been included. These are all much larger at the coarsest lattice spacing than elsewhere.
As a guide to estimating the size of these errors we have estimated, using a potential model [16], the size of the shifts in energy caused by the relativistic and discretisation corrections of δH that we have included and compared those to the results of NRQCD calculations [3]. The potential model estimates for the sum of relativistic and discretisation corrections in δH at β = 6.0 give a resulting shift to the bb 1P − 1S splitting of −10 MeV, made up of −10 MeV from the relativistic corrections (10 MeV from the p 4 term and −20 MeV from the D · E term) and cancelling contributions each of around 10 MeV in magnitude from the two discretisation corrections [16]. The 2S − 1S splitting has smaller shifts because the expectation values on which the shifts depend are more similar for the 2S and 1S than the 1P and 1S. The expected result from adding all the terms is still −10 MeV. The NRQCD results show a 1P − 1S splitting that is 20(30) MeV larger without δH than with, and a 2S − 1S splitting that is 5(25) MeV larger. This is in good agreement with expectations, albeit with statistical errors that are too large to show a clear effect. However, if the shifts had been much larger than the estimates, they would have been visible above the noise.
From this result we can extrapolate to the size of relativistic and discretisation corrections that we have not included. The terms that we can estimate most readily are those involving powers of quark momenta since these are easy to relate to lower order terms. The terms involving chromo-electric and magnetic fields and those terms with a structure that appears for the first time at higher order are much harder and we have not estimated these. There is no reason, however, to suppose from our study above (which compares the p 4 and Darwin terms), that these terms should be any larger than the ones we can easily estimate.
Higher order relativistic corrections would appear as v 6 terms in δH, i.e. O(v 4 ) relative to the leading terms. The percentage error we expect is then naïvely 10% 2 = 1%. On the other hand, the estimates above using potential models [16], of the v 4 spin-independent relativistic corrections that we have included show them each to be less than half of the 10% (= 50 MeV) naïvely expected. In addition we actually need the difference between the corrections for, say 1P and 1S, to get the shift in the splitting. This indicates that higher order corrections could be smaller than 1% too. Another type of similar higher order correction is that from radiative corrections to the c i coefficients beyond tadpole-improvement. These should appear at the level of α s at an ultra-violet scale times the v 4 relativistic corrections, giving 0.5% at β = 5.7 and less at higher β values. To encompass both these higher order physical corrections a 1% error is given as the second error in column 5 of Table V.
The NRQCD action, eq. 3, includes the leading a and a 2 corrections which appear in power counting form as p 2 a 2 v 2 and Kav 2 relative to the leading order v 2 term (H 0 ). p and K are a typical momentum and kinetic energy associated with the bound state. Estimates of these terms using potential models [16] yield shifts of 40 MeV at β = 5.7, 10 MeV at β = 6.0 and 5 MeV at β = 6.2, in the 1P − 1S splitting. Higher order discretisation errors not included could be radiative corrections to those included (beyond the tadpole-improvement of these terms which has been done) and we can estimate these as α s (π/a) [15] times the leading errors. This gives 7 MeV at β = 5.7, 1.5 MeV at β = 6.0 and 0.7 MeV at β = 6.2. Higher order terms in a (such as p 4 a 4 v 2 terms) would give a percentage effect roughly the square of the leading order terms, i.e 4 MeV at β = 5.7, 0.2 MeV at β = 6.0 and essentially 0 at β = 6.2. We should also consider the first discretisation corrections to the first relativistic corrections i.e. terms of order p 2 a 2 v 4 . These are very similar looking terms to the p 4 a 4 v 2 terms and so we can use this to estimate their size. Including powers of Ma and numerical factors we get a similar size correction at β = 5.7. At higher values of β these corrections are actually more important than those which are higher order in a but not suppressed by powers of v 2 . However, all discretisation corrections become smaller at higher β so they are still negligible; 0.5 MeV at β = 6.0 and 0 at β = 6.2. The third error given in column 5 of Table V is then conservatively estimated by the sum of the three errors given above from discretisation errors in the NRQCD action. It is interesting to note that the discretisation corrections are smaller than the anticipated higher order relativistic corrections except at β = 5.7. Note also that the statistical errors are generally larger than the systematic -to see any effect from including higher order terms we would have to improve our statistical error significantly. Higher order discretisation corrections from the gluon action are anticipated to be negligible given the size of the O(a 2 ) correction in Table V.

B. Determining the quark mass
As discussed above, experimental spin-averaged radial and orbital splittings are very insensitive to the value of the quark mass. However, in the quenched approximation there is some mass dependence for these splittings on the lattice [18], and this is increased if spinaveraging is not done. It is therefore true that a tuned bare quark mass is necessary to get the right radial and orbital splittings. The spin splittings are much more sensitive to the quark mass (roughly as its inverse), and it is essential to tune the quark mass to get these correct. The way in which we tune the quark mass is to adjust it until the kinetic mass of the Υ agrees with experiment. Table VI shows the kinetic mass of the Υ (given in lattice units in Table III) in GeV, for each of the 3 values of β, using a −1 from Table V. There is reasonable agreement with the experimental result 9.46 GeV in each case, provided that the value for a −1 is taken from the Υ ′ − Υ splitting. The systematic errors in the determination of aM kin are at the 1% level from the same sources as systematic errors in the determination of a −1 . This is smaller in every case than the statistical error, dominated by the statistical uncertainty in a −1 .
The difference between the kinetic mass of a meson and its energy at zero momentum is calculable in perturbation theory [19]. The formula which relates E N R and the kinetic mass M is : where Z m is the mass renormalisation and E 0 the energy shift. Table VII gives values for Z m and E 0 appropriate to the different bare masses used at each value of β. Tadpole-improved lattice perturbation theory has been used for these parameters and the scale of α s set using the BLM scheme [19]. The values obtained for M Υ from the perturbative expression, (7), are given in lattice units in the sixth column and should be compared with the results for aM kin in Table VI. There should be agreement at all values of β independent of whether the quark mass is well tuned to that appropriate to the b or not, and we see that there is. The perturbative error is at O(α 2 s ) and is taken here as the square of the O(α s ) term. Note that the relationship in eq. 7 is well defined perturbatively. Non-perturbative values for the shift between E N R and M can also be measured on the lattice for, say, Υ and used, divided by 2, for B physics [20].
We can use the information above to determine the b quark mass [21] in two independent ways and ask whether the renormalised b quark mass we obtain scales from one value of β to the next. These results will be presented elsewhere [22].

C. Radial and Orbital splittings
From the results of Table III we can investigate the scaling of spin-averaged splittings across the range of β values we have used. Our analysis of systematic errors in section 3.1 already implied that we do not expect to see violations of scaling in these quantities. Figure 3 shows with circles the ratio of our χ b −Υ splitting to the ρ mass, obtained by the UKQCD collaboration at the same three values of β [11,23]. The ρ mass was calculated using an action with O(a) errors removed using a clover term with tadpole-improved coefficient; i.e. a similar philosophy to that used here to remove discretisation errors from the Υ splittings. The remaining errors in the ρ determination may then be of O(a 2 ). For the 1P −1S splitting we have removed the O(a 2 ) gluonic errors so, as discussed earlier, the remaining errors are higher order in a and/or suppressed by powers of v 2 . The plot is very flat showing that good scaling is obtained. We can rule out scaling violations in this ratio with a scale larger than a few hundred MeV, and it would be very unlikely for the ρ mass and the Υ system to have identical and cancelling large scaling violations. The plot also shows in the comparison to the GF11 results, that if an unimproved calculation is done of the ρ mass [24], an absence of scaling is quite evident (in this case linear in a).
The ratio ∆(χ b − Υ)/m ρ shows a big discrepancy with experiment. This we believe is an error from the quenched approximation. The scales intrinsic to a light hadron system and the Υ are quite different and the coupling constant does not run correctly between these scales in the quenched approximation. Again, these errors are masked by discretisation errors for an unimproved light hadron spectrum. Figure 3 also shows the ratio of the χ b − Υ splitting to the Λ parameter of QCD in the V scheme [9]. Again good scaling is seen over this limited range of lattice spacings. Figure 4 shows the dimensionless ratios of various splittings within the Υ system (spinaveraged where possible) as a function of lattice spacing. The ratios shown are for the Υ ′′ − Υ, Υ ′ − Υ and the h ′ b − Υ (experiment for the latter case using χ . The ratios are constant as a function of lattice spacing, although the statistical errors are rather large for the 3S and 2P cases. This means that we can interpret our results as continuum results. Note that our current value for the energy of the 2 3 S 1 at β = 6.0 has dropped by 1σ from our previous calculation [3] with slightly smaller statistical errors. Our 1 1 P 1 energy has not changed and the error has fallen by a factor of 2. Our best value for the ratio (2S − 1S)/(1P − 1S) is now 1.35 (6) and the discrepancy with experiment remains (1 − 2)σ, controlled by the uncertainty in the 2 3 S 1 energy. The discrepancy with experiment for the ratios involving the 3S and 2P could be affected by finite volume errors as well as quenching errors.

D. Spin splittings
Spin splittings in lattice units are given in Table III. This includes, for the first time, the radially excited hyperfine splittings, Υ ′ − η ′ b and Υ ′′ − η ′′ b at β = 6.0 where our results are most precise. No hyperfine splitting has yet been seen experimentally for the Υ system, and a radially excited hyperfine splitting is still awaited for charmonioum. We find the 2S hyperfine splitting to be about half that of the ground state here. A proper analysis however, has to study the lattice spacing dependence and the effect of unquenching on this ratio.
In this section we investigate the dependence in physical units of the ground state spin splittings on the lattice spacing. We also compare to recent UKQCD results [25] in which a systematic study at next order in relativistic and discretisation corrections has been done for spin-dependent terms, i.e. those which affect the spin splittings, only. By comparing the scaling of our results with their results we can untangle to some extent the difference between the discretisation errors at fixed a, which are unphysical, and the relativistic corrections, which are physical. Whether or not the spin-dependent relativistic corrections (i.e. terms of O(Mv 6 )) are sizeable or not is important since these are terms outside the scope of potential models.
To investigate scaling we show in Figures 5 and 6 various spin splittings in MeV setting the scale from the Υ ′ − Υ splitting as a function of a 2 . The Υ ′ − Υ splitting is used since this is the one for which our quark masses are closest to being tuned and this is important for spin splittings [26]. The results from Table III are given as squares. We include in the Figures the direct uncertainty in the spin splittings from statistical errors in the lattice spacing. There is an additional error, which we have not included, from the uncertainty in the kinetic mass from the error in the lattice spacing. This uncertainty leads to an uncertainty in the tuned quark mass which affects the spin splittings. If the spin splittings vary as the inverse quark mass (roughly true at least for the hyperfine splitting), then this second uncertainty is actually equal in size, correlated with and in the same direction as the error that we have included. This would lead to an approximate doubling of the error bars in the Figures.
Clear dependence on the lattice spacing is visible for the hyperfine splitting, Υ − η b . For the 1P fine structure it is less clear because statistical errors, particularly at β = 6.2, are so much higher. Lattice spacing dependence in our results is not surprising since the Hamiltonian we have used has only leading order spin-dependent terms with no corrections for discretisation errors in these terms. These errors are expected to be O(a 2 ) relative to the leading term and therefore in relative terms significantly larger than for our spinindependent splittings. In addition, in a potential model picture the fine structure is provided by potentials which are generally of much shorter range than the central potential so we would expect the fine structure to have a harder scale for discretisation errors than the radial and orbital splittings. If we parameterize the scaling violations in the hyperfine splitting by [7]: we find the parameter µ to take a value around 0.8 GeV for the squares. For the 1P fine structure we have the advantage that we can compare to experimental results as well as comparing results at different values of the lattice spacing. In Figure 6 we notice a clear disagreement with experiment for the squares in that both the overall scale of splittings is too small and the ratio is too large. Similar results were noticed in our charmonium spectrum calculation [27]. The disagreement with experiment over ρ was suggested there as a discretisation effect and results here tend to confirm that discretisation errors do have an effect in this ratio. Our value for ρ at beta=5.7 is 1.4(3) and at 6.0, 1.1 (3). The experimental value is 0.66. The comparison between the squares at β = 5.7 and 6.0 shows that scaling violations are just about visible for the 3 P 0 and 3 P 2 above the statistical errors. The slope of the scaling violations is similar to that for the hyperfine splitting (but with large errors). Also in agreement with the hyperfine splitting we can see that adding discretisation corrections would increase the splittings and this would make the 1P fine structure closer to experiment.
On Figures 5 and 6 we also show results from Manke et al [25] at β = 6.0 and 5.7 with diamonds. They have used the same bare quark masses as us at both values of β but included higher order spin-dependent relativistic (O(v 6 )) and discretisation (O(a 2 v 4 )) corrections to the action [5], all tadpole-improved with u 0P . We take their spin splittings in lattice units and convert to physical units using our Υ ′ − Υ splitting since, as explained earlier, this is the a −1 for which the quark masses are best tuned. Our result for this splitting should be the same as theirs since they have added only spin-dependent terms [28] and ours is more precise. The lattice spacing dependence of their results should be much reduced over ours, and we see that it is. By comparing their results and ours, we can see that the relativistic corrections act in the opposite direction to the discretisation corrections, in agreement with the findings of [29].
By comparing the results of [25] with ours we can attempt to untangle the size of v 6 relativistic corrections to the fine structure. First, however, we should also include radiative corrections to the leading order terms, beyond tadpole improvement. These are of O(α s v 4 ) and therefore of a similar size to the relativistic v 6 corrections. The c i in equation 3 take the form 1 + c (1). Preliminary perturbative results are available [30] for c (1) 4 , for the σ · B term in the action, and confirm this expectation. In perturbation theory the hyperfine splitting then picks up a factor of the square of this correction, so it is quite sensitive to c (1) 4 . Figure 5 shows with fancy squares our results for the hyperfine splitting rescaled by ( Values for α P were taken from [15]. c (1) 4 increases as Ma increases, which gives a larger correction at low β. Since this is also the direction in which α P (3.4/a) increases the correction at β = 5.7 is large, 50%. At β = 6.2 the correction from c For the hyperfine splitting, we will assume that the discretisation corrections at β = 6.0 in the (v 6 , a 2 v 4 ) action [25] are of O(a 4 ) and therefore negligible already at β = 6.0. For the (v 4 ) action (this work) we will take the result at β = 6.2 as being the least affected by discretisation errors (we don't expect more than a 5% effect here). The difference then between the result in physical units between β = 6.2 (v 4 ) and β = 6.0 (v 6 , a 2 v 4 ) represents the effect of relativistic corrections. From Figure 5 with fancy squares and diamonds this gives a 15% decrease in the hyperfine splitting and from the plain squares and diamonds a 20% decrease, consistent with the expectation of around 10% from power counting in v 2 /c 2 .
For the p fine structure, the radiative corrections are not available and they would in any case be much harder to apply without doing the full lattice calculation including them. We therefore just compare the results with (diamonds) and without (squares) the v 6 and a 2 v 4 corrections in Figure 6. The lattice spacing dependence is apparently reduced by the a 2 v 4 corrections, although the statistical errors make it hard to quantify this result. The relativistic corrections again act to reduce the splittings, now to something like 50% of experiment. The size of the relativistic corrections is significant, possibly several times the naïve estimate of 10%. They could be offset by radiative corrections, however, in a fully consistent calculation that includes the next order in v 2 , a 2 and α over our calculation.
Radiative corrections which increase the spin splittings, as for the hyperfine case, would improve the agreement with experiment, but unless these are surprisingly large, the splittings will remain too small in the quenched approximation. Initial unquenched results for the v 4 action [15] and the v 6 , a 2 v 4 action [32] have indicated an increase in the p spin splittings, but currently statistical errors are too large to clearly determine the size of unquenching effects.
Ref. [25] points out that the (v 6 , a 2 v 4 ) action has a much improved value of ρ (eq. 9) over the (v 4 ) action. They quote 0.56 (19) at β = 6.0 in the quenched approximation, now in agreement with experiment, confirming the discretisation effects for this ratio discussed above. In fact, we do not expect this ratio to agree exactly with experiment in the quenched approximation. An analysis of ρ for the simple Cornell potential [33] shows that for a quenched α s which is too small, the value of ρ is reduced. With better statistics it may become clear that the quenched value for ρ is actually smaller than the experimental result.
Another way to study the sensitivity of the fine structure to higher order terms in the action is to change the value of the tadpole coefficient, u 0 . Several calculations [1] have indicated that the average link in Landau gauge, u 0L is to be preferred over the fourth root of the plaquette, u 0P , that we have used here. Changing u 0 is equivalent to resumming various terms to all orders in the c i . It is again easy to rescale the hyperfine splitting after the simulation to account for a different value of u 0 , by counting the powers of u 0 that appear in front of the σ · B term. The hyperfine splitting is very sensitive to u 0 , scaling as 1/u 6 0 , when allowing for u 0 factors in H 0 which change the renormalised quark mass. The expected scaling is confirmed by calculations at different values of u 0 [3] (although here a fixed bare quark mass was used). The hyperfine splitting increases most on using u 0L at β = 5.7. There the ratio of u 0P /u 0L is 1.045; at β = 6.0 it is 1.020 and at β = 6.2, 1.013 [34]. The hyperfine splitting with u 0L then has somewhat reduced lattice spacing dependence. Note that c 4 will be different for the action with u 0L , since a different amount of the radiative correction has been absorbed by the tadpole factors. c (1) 4 in fact is larger, by 0.24, being the cube of the ratio of the u 0 factors to O(α). The size of relativistic and discretisation corrections will change by differing amounts on changing u 0 because they have different effective powers of u 0 in their coefficients.
From the above we can conclude that systematic errors in the fine structure still exceed 10%, even in the quenched approximation. To improve the accuracy the v 6 and a 2 v 4 action must be used, with radiative corrections included for the leading spin-dependent terms. For the hyperfine splitting the sensitivity to the c i is such that the O(α 2 ) term also needs to be known (or at least bounded). The determination of the c i is in progress [30]. More accurate meson kinetic masses and lattice spacing determinations are also required so that the quark mass can be tuned more accurately, since the spin splittings are sensitive to this.
The wavefunction at the origin is another short distance quantity but one which does not vanish as the spin-dependent terms are switched off. In Figure 7 we show results for ψ(0) for the Υ and Υ ′ as a function of a 2 , again setting the scale from the Υ − Υ ′ splitting since the wave function is sensitive to quark mass. The lattice spacing dependence does not seem severe but the statistical errors are large. They are dominated by those from a −1 which here is raised to the 3/2 power. The bursts show experimental values based on the leptonic widths and the naïve Van Royen-Weisskopf formula, where M v is the vector meson mass. There will be radiative corrections to this formula of O(α s ). For potential model wavefunctions these radiative corrections have coefficient 16/3π but this does not have to be the same for NRQCD. Calculations of the renormalisation factors have only been done for lower order actions than the one used here [35]. There will also be terms of O(1/M b , α s /M b ) from mixing of the vector current with higher dimension operators (in a similar way to the mixing for f B which has been worked out in some detail [20]). A lot of the radiative corrections cancel in the ratio of leptonic widths for the 1S and 2S states. Experimentally this ratio is 2.5 (2). Our result at β = 6.0 simply from the ratio of |Ψ(0)| 2 is 1.7(4). We expect errors in the quenched approximation which would suppress the 1S wavefunction at the origin relative to the 2S, so that it is not surprising that our result is somewhat smaller than experiment. We cannot yet say whether the reduction observed is reasonable.
To obtain results for the real world from simulations in the quenched approximation and at 2 flavours of dynamical quarks, extrapolations in N f (= the number of flavours of dynamical quarks) will be required. For this it is important to have scaling quantities, as we have discussed here. It is also necessary to use quantities which are sensitive to the same number of dynamical fermions. This means that for a given splitting we should set the scale using a quantity which has similar internal momentum scales. For example, for the fine structure it may not be optimal to use the softer radial and orbital splittings to set the scale. Typical momentum scales for radial and orbital splittings are around 1 GeV so they can be expected to 'see' 3 flavours of dynamical fermions in the real world [15]. The fine structure may instead 'see' 4 flavors and this would imply that an extrapolation in N f can only be done in terms of a quantity which also sees 4 flavors.
In Figure 8 we show a plot of the hyperfine splitting in MeV against the square of the lattice spacing using the 1 3 P 2 − 1 3 P 0 splitting to set the scale. Statistical errors are now very large, despite using the most accurate of the P fine structure splittings. Much higher accuracy will be needed to use this ratio for N f extrapolations.

IV. CONCLUSIONS
We have calculated the spectrum of bottomonium bound states at 3 different values of the lattice spacing using NRQCD in the quenched approximation. Ratios of spin-independent splittings are constant in this region, independent of the lattice spacing. This is a necessary requirement for the results to make physical sense. There is no necessity to extrapolate to vanishing lattice spacing (impossible anyway for NRQCD). The constancy of the results implies that they can now be used in extrapolations to real-world N f values when combined with results on dynamical configurations [15,32].
For spin-dependent splittings there is some visible (unphysical) lattice spacing dependence. This is not surprising since our spin splittings are only accurate to leading order, given the terms in the Hamiltonian that we have used. A comparison to a calculation [25] which does include higher order relativistic and discretisation corrections to spin-dependent terms shows that the relativistic spin-dependent corrections could be as large as 20% for the hyperfine splitting where they can be quantified. They act to reduce spin-dependent splittings in all cases, and this increases the disagreement with experiment for the P fine structure in the quenched approximation.
For extrapolations to real world values of N f we may need to employ different techniques for fine structure to that for radial and orbital splittings since different momentum scales are involved and this may lead to a different number of dynamical flavors being 'seen'. Statistical errors are currently too high to use a P spin splitting to set the scale in such a programme. We thank Sara Collins, Joachim Hein, Thomas Manke and Howard Trottier for useful discussions. CTHD is grateful to the Institute for Theoretical Physics, UC Santa Barbara for hospitality and to the Leverhulme Trust and the Fulbright Commission for funding while this work was written up.  (Table V).
FIG. 2. Effective mass plots for the ( 1 P 1 ) 1l correlator at all three values of β, in order with β = 5.7 at the top. The time axis has been converted to physical units of GeV −1 using the χ b − Υ splitting to set the scale (Table V). FIG. 3. Dimensionless ratios of the χ b − Υ splitting to the parameter Λ V (diamond) and to the UKQCD ρ mass [11] (circle) against the lattice spacing in fm (set by the χ b − Υ splitting). Results using the GF11 ρ mass [24] are given by squares. The burst represents the experimental value for ∆(χ b − Υ)/m ρ . Plain squares indicate our results from Table III. The diamonds indicate the results from ref. [25] using a higher order action. They are at matching values of β but offset slightly for clarity. Fancy squares indicate our results, rescaled by the square of c 4 calculated to O(α) in [30]. The fancy diamonds indicate the result from ref. [25] shifted by the same amount as our results to account for radiative corrections to the σ · B term. The error bars shown are statistical only and include some of the error from the uncertainty in the scale (see text). The x axis errors from uncertainty in the scale are shown only for the squares for clarity.