Splittings of low-lying charmonium masses at the physical point

We present high-precision results from lattice QCD for the mass splittings of the low-lying charmonium states. For the valence charm quark, the calculation uses Wilson-clover quarks in the Fermilab interpretation. The gauge-field ensembles are generated in the presence of up, down, and strange sea quarks, based on the improved staggered (asqtad) action, and gluon fields, based on the one-loop, tadpole-improved gauge action. We use five lattice spacings and two values of the light sea quark mass to extrapolate the results to the physical point. An enlarged set of interpolating operators is used for a variational analysis to improve the determination of the energies of the ground states in each channel. We present and implement a continuum extrapolation within the Fermilab interpretation, based on power-counting arguments, and thoroughly discuss all sources of systematic uncertainty. We compare our results for various mass splittings with their experimental values, namely, the 1S hyperfine splitting, the 1P-1S splitting and the P-wave spin-orbit and tensor splittings. Given the uncertainty related to the width of the resonances, we find excellent agreement.


I. INTRODUCTION
Over the past decade, the experimental study of the products of B-meson decays has led to the discovery of a wealth of excited charmonium states. Many of them present interesting challenges for theoretical interpretation. Because lattice quantum chromodynamics (QCD) is an ab initio method for studying hadron spectroscopy, in principle, it should provide a guide to the interpretation of these states [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. To address these questions with confidence, it is important that lattice discretization (cutoff) effects be under control. The more limited objective of the present work is to carry out a high-precision study of the splittings of the lowlying charmonium states-particularly the 1S and 1P states-and, thus, lay the foundation for further calculations of excited states. Spin-dependent mass splittings are expected to be extremely sensitive to the charm-quark mass and to heavy-quark discretization effects. Reproducing these delicate splittings can therefore serve as another demonstration that systematic uncertainties are under excellent control.
Our effort follows a previous analysis campaign on the same gauge configurations [2]. Preliminary results have been reported [17] and some additional details about our quark sources can be found in Ref. [18]. Our new results supersede the results in those publications. Although other groups have reported partial results for the low-lying charmonium spectrum [1,[19][20][21][22], there are no systematic, high-precision studies for any action. Thus, to our knowledge, our campaign is the first that includes precise tuning of the charm-quark mass, precise determination of the lattice scale, and a controlled extrapolation to physical light sea-quark masses and zero lattice spacing. Our paper is organized as follows. In Sec. II, we further describe the objectives of our current work, while we describe our lattice setup in detail in Sec. III. Section IV shows our results for the splittings among low-lying charmonium states including a chiral and continuum extrapolation of results and a full error budget. We summarize our findings in Sec. V where we also provide a brief outlook.

II. THEORETICAL BACKGROUND
In this paper, our main objective is the QCD determination of the masses of the 1S and 1P states in the charmonium spectrum. From the ground state masses in the quantum number channels corresponding to the 1S and 1P states, we calculate the hyperfine splitting between the 1S triplet and singlet states the spin-average 1P-1S splitting

3)
M 1S = 1 4 (M ηc + 3M J/ψ ), (2.4) and the spin-orbit, tensor, and 1P hyperfine splittings among the P-wave states ∆M spin-orbit = 1 9 (5M χ c2 − 3M χ c1 − 2M χ c0 ), (2.5)  It is these splittings, extrapolated to zero lattice spacing and physical sea-quark masses, that we compare with their experimental values. These combinations are of phenomenological interest in constructing the heavy quarkonium potential, since they correspond to separate terms in the potential derived from the heavy-quark limit [23,24]. (2.9) Thus, their dependence on the lattice spacing provides useful information about discretization effects in each of the relevant terms, as discussed in detail in Ref. [2]. Table I lists the 1S and 1P states along with their masses and widths, as determined from experiment [25]. While some of these states are extremely narrow, both the η c and χ c0 have a non-negligible hadronic decay width, resulting from charm-anticharm annihilation. In lattice QCD, this effect comes from disconnected diagrams, which our current simulation omits. That is, we treat all low-lying charmonium states as stable. It is therefore not a priori clear whether we will obtain good agreement with the η c and χ c0 masses. This shortcoming complicates the comparison with experiment, in particular for the 1S hyperfine splitting. 1 We further comment on this issue when comparing our results with previous results in Sec. IV F.

III. METHODOLOGY
This section presents the methodology for the lattice determination of the charmonium masses described in Sec. IV. In addition to our lattice setup, our procedures to deal with uncertainties from the mistuning of the charm-quark mass, our strategy for the chiralcontinuum fits, and the systematic uncertainty arising from the determination of the lattice spacing are discussed in this section. TABLE II. Parameters of the MILC ensembles used in this study. Listed are the lattice spacing a, the ratio of the sea-quark masses m l /m s used in the simulation, and the lattice size L 3 × T , Also included are the number of source time slices used in the calculation N src , the tuned charm-quark hopping parameter κ c , the charm-quark hopping parameter of our simulation, κ c , and a citation for the ensemble. The first uncertainty in κ c is statistical, and the second is from the uncertainty in the lattice scale. A. Gauge configurations We use the (2+1)-flavor gauge configurations generated by the MILC collaboration [42] with the asqtad fermion action for sea quarks. The ensembles used in this work are listed in Table II. The use of five different lattice spacings a and two different light sea-quark masses (given in the table as a fraction of the strange quark mass in the simulation) enables us to perform a controlled chiral-continuum extrapolation. Four source time slices per gauge configuration are used, for a total of approximately 2000 to 4000 sources per ensemble. We use the Fermilab prescription [43] for the charm quarks, which suppresses heavy-quark discretization effects in mass splittings [44]. The charm-quark hopping parameter κ c has been tuned by demanding that the D s kinetic mass be equal to the physical D s meson mass in the way described in Ref. [45]. The resulting κ c and the (sometimes slightly different) simulation value κ c are also given in Table II. Note that we refer to the quark masses used in the simulation as m l and m s while we denote the physical light-and strange-quark masses by m l and m s . When calculating observables, we need to take into account this difference in our chiral-continuum extrapolations.

B. Calculation of observables
We calculate a matrix of correlators C(t) using quark-antiquark interpolators with the J P C quantum numbers of the states in question, where J is the total spin and P and C are parity and charge conjugation quantum numbers. We opt for a basis built from interpolators with derivatives and use interpolating operators similar to those suggested by Liao and Manke [46], which have also been used by Dudek et al. [1]. A subset of similar interpolators has also been used in Ref. [47] and similar interpolators using displacements or full pla-quettes rather than derivatives have previously been considered in Ref. [48]. Disconnected contributions, where a valence charm-anticharm-quark pair annihilates, are omitted when calculating the correlators.
Our operators are constructed from stochastic wall sources, including covariant Gaussian smearing. Stochastic sources consist of a four-component-spinor field on a single time slice with random color orientation, but definite spin: where r labels the stochastic source, β its spin, and a and α are the 12 Dirac color and spin components. Averaged over a sufficiently large number of stochastic sources N r , we have With both charm and anticharm quarks originating from the same source, or with one source modified by Gaussian smearing, the stochastic average gives the effect of charmonium sources composed of local or smeared bilinears of the form where the smearing operators are included in the definition of O i . All links appearing in the Gaussian smearing operators and in the covariant derivatives below are smeared with a fixed number of APE-smearing [49] steps. Gaussian smearing is implemented by acting with a smearing operator M on the stochastic sources S to obtain Gaussian sources: where ∆ is a covariant 3D Laplacian, N is just a normalization factor, and σ/a and N are chosen such that M approximates a Gaussian with (physical) standard deviation σ in coordinate space. Thus, lim Table III lists the smearing parameters for both the gauge link smearing and for the Gaussian quark sources.
In the constructions discussed in Appendix A, we use the following operators: Here M is the Gaussian smearing operator defined in Eq. (3.4a), and P i is a derivative-type operator on a given time slice t,  where r is kept of roughly the same length in physical units and W i (x, t; x + rî, t) denotes the shortest Wilson line connecting (x, t) and (x + rî, t). In Eqs. (3.6), the P i act to the right. The continuum version of operator B i has a relation to the chromomagnetic parts of the field strength tensor To avoid an (anti)symmetrization of the derivatives, which would require more sources, we first apply derivatives and then the Gaussian smearing. A detailed discussion of this approach can be found in Ref. [50]. We use the variational method [51][52][53][54], solving the generalized eigenvalue problem with reference time slice t 0 . The ground state mass can be extracted from the large time behavior of the largest eigenvalue. For this we use (multi)exponential fits to the eigenvalues in the interval [t min , t max ], taking into account correlations in time separation. At fixed t 0 , ∆E k is formally given by while for the special case of t ≤ 2t 0 and a basis of N correlators [55] ∆E k is given by We investigate the dependence of our results on t 0 and find that in practice a rather small value of t 0 provides the best compromise between excited-state contaminations and statistical uncertainty. Here and elsewhere, the statistical uncertainties are computed from a singleelimination jackknife. In our analysis, the reference time t 0 and the lower boundary of the fit window t min are kept roughly constant in fm for the 1S and 1P states respectively. 2 The upper boundary of the fit-window t max is chosen such that the eigenvectors ψ (k) remain stable within statistics in the whole fit range, which in some cases results in a somewhat shorter fit windows than just considering plateaus in the effective masses. For the P-wave states on the coarsest lattice spacing, where the tuning of the quark-smearing was performed, remaining excited-state contaminations are extremely small, and we need to use loose priors on the mass splittings between the ground state and the lowest excitations in order to avoid clearly unphysical fit results with two almost mass-degenerate ground states. In some cases increasing the size of the basis used in the variational method leads to no improvement in the ground state but adds statistical noise. For our final results we therefore opted to suitably prune the interpolator basis, and we list our choices of basis in Appendix A.

C. Charm-quark-mass corrections
For some of the ensembles listed in Table II the charm-quark hopping parameter of the simulation κ c differ slightly from the physical charm-quark hopping parameter κ c . The raw splittings on these ensembles have to be corrected for this mistuning. To determine the needed correction, we compute the derivative of each mass splitting with respect to κ c on one ensemble, namely the one with a = 0.114 fm and m l /m h = 0.1 [32]. We assume that once the slope is expressed in terms of physical quantities, it remains the same for that mass splitting for all ensembles. Since the adjustments are small, any residual lattice spacing dependence in the slopes should be negligible.  To be explicit, for mass splitting ∆M i , we assume that the following derivative is the same for all ensembles: where m 2 (κ c ) is the kinetic mass of the charm quark. For a given κ c we estimate that mass from the ensemble's critical hopping parameter κ crit and tadpole factor u 0 using the tree-level expressions [Eq. (4.9) of [43]]: The correction to the mass splitting, resulting from a shift dκ c is then given in r 1 units [56] by Values of κ crit , u 0 , and A for each ensemble are listed in Table IV. A quantitative estimate for the uncertainty from this procedure is provided in Sec. IV E.

D. Chiral and continuum fits
We perform a combined extrapolation to the continuum values and to physical light-and strange-quark masses. Our data indicate a clear sea-quark mass dependence for some of the observables, 3 which means that we also need to take into account the effect of mistuned [58], and the notation for the terms follows that reference. Values of α s consistent with those in Table V have been used. The terms arising from mass mismatches are denoted in the plot by the masses in the short-distance coefficients. In addition a rotational symmetry breaking term (with coefficient w 4 ) is important for the 1P-1S splitting. Expressions for the short-distance coefficients can be found in Ref. [58].
strange sea-quark masses. Our model for the lattice spacing dependence is based on the Oktay-Kronfeld [58] analysis of the Fermilab prescription, which provides NRQCD powercounting [57] estimates of various heavy-quark discretization effects in quarkonium. They are parameterized as mass mismatches, leading to functions f i (a) of the lattice spacing that are determined separately for each observable. In addition to the terms for the heavy-quark discretization effects, we also add a term linear in α s a 2 as appropriate for the asqtad sea quarks. For our combined sea-quark mass and continuum fit we use the Ansatz 18) as our fit model. The values for m l , m l , m s and m s are given in Table V along with the values of the renormalized coupling in the V scheme [59] α s at scale 2/a used in the analysis of discretization effects. For each observable we determine the most important mass mismatches arising at O(v 4 ) and/or O(v 6 ) in NRQCD power-counting. Figure 1 shows the expected discretization uncertainties from power counting estimates for the splitting indicated in the respective figure. The plotted curves correspond to c i = 1. For the 1P-1S splitting, this includes a term from rotational symmetry breaking (w 4 term). In our default fits we use Bayesian priors centered around 0 with a prior uncertainty of 1 as a constraint for all terms originating from heavy-quark discretization effects. As part of our systematic variations described in Sec. IV E, this prior uncertainty is varied. In addition to these terms we also allow for a generic α s a 2 term (without prior) characteristic of light-quark discretization effects. We discuss the relevant mass mismatches for a given splitting when we present our results in Sec. IV. For each observable we compare continuum extrapolations with just two terms (α s a 2 and the leading heavy-quark discretization term) and with three terms (the α s a 2 term and the leading and subleading heavy-quark discretization terms). We further check the variation from replacing the α s a 2 term by an a 2 term. While a single leading shape is usually enough to get a good fit of the data, including further possible shapes leads to a larger and more realistic uncertainty estimate. The fit variations described above are among the fit variations shown in Sec. IV E, where our error budget is also discussed.

E. Scale-setting uncertainty
For the figures presented in Sec. IV, we use MILC's version of the Sommer scale, r 1 [56]. The values of r 1 /a for the asqtad ensembles and an explanation for our value r 1 = 0.31174(216) fm can be found in Ref. [60]. This value was determined in the "massindependent" scale-setting scheme, the one adopted here. To estimate the scale-setting error, for each observable we first determine the result using the central value for both r 1 and κ c and then repeat the procedure, shifting r 1 by one standard deviation while simultaneously shifting the tuned κ c by an amount that results from the same shift in r 1 . The scale-setting uncertainty for each observable is discussed in Sec. IV E and tabulated in Table VII.

IV. RESULTS
In this section, results for the mass splittings from Sec. II are presented. After discussing each splitting in turn, the systematic uncertainties associated with the determination are quantified and the resulting values are compared with the results from previous determinations.

A. 1S hyperfine splitting
Like all spin-dependent splittings, the 1S-hyperfine splitting is highly sensitive to heavyquark discretization and charm-quark tuning effects. As such, it is an important benchmark quantity for lattice-QCD calculations of the charmonium spectrum.
For the 1S hyperfine splitting, autocorrelations in the Markov chain of gauge configurations are significant and need to be taken into account. To do so, we estimate the integrated  Due to the mistuning of the strange quark in the sea, which differs from ensemble to ensemble, the data points appear away from the curves. To illustrate that the data are well described by the fit, the black crosses show the fit results evaluated at the lattice parameters of the gauge ensemble. The magenta symbol indicates the result in the combined chiral and continuum limit.
autocorrelation time using two methods. Method one is to determine the autocorrelation time from the jackknife sample using the method and software of Wolff [61]. An alternative consists of constructing binned data from the jackknife estimates of the unbinned set and extrapolating the results for bins of sizes 1 to 5 to infinite binsize using the expected scaling. We determined the integrated autocorrelation time using both methods and check the results for consistency. The two methods agree excellently and, as the uncertainty estimate on different ensembles is independent, we use the second method to inflate the statistical uncertainties on a single ensemble appropriately. 4 Figure 2 shows the results for the 1S hyperfine splitting, along with a chiral-continuum extrapolation of the results. Where needed, the data have already been shifted for mistuning of the charm-quark hopping parameter, as outlined in Sec. III C. Note that significant contributions from charm-annihilation diagrams to this observable are expected [62]. When comparing our final results with the experimental value in Table IX we use the determination of −1.5 to −4 MeV from Ref. [62] as an estimate for the uncertainty from neglecting disconnected contributions.
The leading heavy-quark discretization effects contributing to the hyperfine splittings come from mismatches of m B and m 2 . Following Ref. [58], we use NRQCD power counting with v 2 = 0.3 and m c = 1400 MeV along with the tree level formulas from Ref. [43] to estimate the expected size of all heavy-quark discretization effects. The relevant formula for m B is Eq. (4.22) of Ref. [58], and the shape of the resulting mismatch is plotted in the first pane of Fig. 1. Note that our fermion action includes a clover term [63] with the tadpoleimproved tree-level value c B = c E = u −3 0 , where u 0 is the average link from the plaquette. A same as "default" but using sea-quark discretization effects of order a 2 rather than α s a 2 B results when omitting the lattice data at the coarsest lattice spacing C results when omitting the lattice data at the finest lattice spacing D result using just terms of order α s a 2 and a single shape for the heavy-quark discretization effects E heavy-quark discretization effects with priors for c i half of the default width (0 ± 0.5) F heavy-quark discretization effects with priors for c i double the default width (0 ± 2) G 1σ variation of the κ c slope used to shift data to physical κ c This contribution is therefore suppressed relative to m c v 2 (the kinetic energy of the meson) by a factor 1 2 α s v 2 c . The sign of the contribution is, however, not known. The next largest heavy-quark discretization effects come from mismatches of m B and m 2 , where the relevant formula for m B is given by Eq. (4.23) of Ref. [58]. Again, the resulting estimate of discretization effects from the mismatch is plotted in the first pane of Fig. 1. Note that at tree-level and with only terms up to dimension 5 in the action, this mismatch is the same as the one from the difference between m 4 and m 2 (see below) but it is of a higher order in the NRQCD power counting and therefore suppressed by 1 8 v 4 with respect to the kinetic energy. For our final fits we use both of these mass mismatches with priors for the coefficients c i from Eq. (3.18) given by 0 ± 1 as well as an unconstrained α s a 2 term. The expected shapes for the mismatches are plotted in Fig. 1.
Finally, the stability of results with regard to systematic variations of the chiralcontinuum fit needs to be assessed. Table VI describes a number of important fit variations (A-G), and their effect on the 1S hyperfine splitting can be seen in the first pane of Fig. 7 in Sec. IV E, below. One of these variations (D) consists of limiting the continuum extrapolations to just two shapes (leading heavy-quark mismatch and sea-quark term). From the difference between the default value and D it can be seen that the central value is largely unaffected while the uncertainty estimate from the fit with leading and subleading shapes is more conservative. Note also that the fit results are stable when the default prior widths are doubled (variation F in Table VI and Fig. 7). While the results are stable when omitting the finest lattice spacing (variation C) there is a somewhat significant shift when excluding the coarsest lattice spacing (variation B). Therefore we take the difference between B and the default fit model as an additional systematic uncertainty. For our final uncertainty estimate provided in Table VII the uncertainties from the scale determination (direct and through the uncertainty in the charm-quark mass) and from the correction of the data for simulation at unphysical charm-quark mass are non-negligible.
B. 1P-1S splitting Figure 3 shows our result for the splitting between the spin-averaged P-and S-wave states ∆M 1P-1S . As in the 1S hyperfine splitting, significant effects from mistuned strangequark masses are visible in our data. Unlike the hyperfine splitting, there is no statistically significant autocorrelation in the Monte Carlo chain, and we therefore treat the data as uncorrelated. In this case, we find large discretization effects, emphasizing the need for several lattice spacings. Having normalized the kinetic energy correctly, we expect leading heavy-quark discretization effects of order v 4 in NRQCD power counting and we plot the expected shapes of the discretization effects in Fig. 1. The terms from the mismatch of m 4 and m 2 and the rotational symmetry breaking term arising at order p 4 are of about equal size. The relevant formulae for w 4 and m 4 are Eqs. (4.4) and (4.5) of Ref. [58]. We also consider the discretization effects from the mismatch of m E and m 2 , where m E is given by Eq. (4.17) of Ref. [58], and we evaluate m E for c E = c B = 1. Again we use Bayesian priors with default value 0 ± 1 for the coefficients c i in Eq. (3.18) associated with heavy-quark discretization effects.
The chiral-continuum fits are stable under all variations shown in Table VI. The effect of these variations is illustrated in the second pane of Fig. 7. In particular the central values do not change when the prior width is increased. As for all other splittings, our final result based on the default fit model also takes into account possible discretization effects of order α s a 2 . The largest variation with respect to this fit model occurs when replacing this term by an a 2 term, which is not motivated by the sea quark action used.

C. Spin-dependent P-wave splittings
The P-wave spin-orbit splitting-shown in Fig. 4-exhibits only small discretization uncertainties. The leading heavy-quark discretization effects come from the mismatch of m E and m 2 and are of order v 4 in NRQCD power counting. Subleading effects of oder v 6 come from the mismatch of m EE and m 2 and from terms of mass-dimension eight not considered in Ref. [58]. Our default fit employs Bayesian priors given by 0 ± 1 for all relevant shapes from Fig. 1.
The results for the spin-orbit splitting are very stable with respect to the variations of the chiral-continuum fit in Table VI pane of Fig. 7. The P-wave tensor splitting (Fig. 5) receives heavy-quark discretization effects from the same mass mismatches as the 1S hyperfine splitting, and the observed total discretization effects in the 1P tensor splitting are of the same absolute size as those in the 1S hyperfine splitting. Variations of our fit model are displayed in the fourth pane of Fig. 7.
As in the case of the hyperfine splitting, we take the difference between variation B (omitting the coarsest ensembles) and the default fit model as an additional systematic uncertainty. The 1P hyperfine splitting defined in Eq. (2.7) is expected to be very small and, indeed, experiments measure a value compatible with zero. Our results are shown in Fig. 6. Our data for this quantity are rather noisy. We find a central value slightly more than 1σ away from zero, but we do not believe this extrapolation to be fully under control, and the strong cancellation may make this combination sensitive to charm-anticharm annihilation.

D. 2S-1S splitting
Beyond the mesons listed in Table I, the only known charmonia below theDD threshold are the ψ(2S) and η c (2S). With our interpolator basis these states are not well determined. Furthermore, we do not include theDD scattering states in our basis and the threshold states therefore cannot be cleanly separated from the close-to-threshold 2S bound states. As a result, the energy values we obtain depend strongly on the lower boundary t min of the fit range, as demonstrated previously in [17]. This issue is not seen in a recent simulation of the ψ(3770) resonance using a more sophisticated basis of both quark-antiquark andDD interpolators [13], where the QCD bound state corresponding to the ψ(2S) can be obtained to a good statistical precision. Note, however, that Ref. [13] was limited to just two sets of gauge configurations, so that the chiral and continuum limits could not be taken.

E. Uncertainty estimates
To obtain final best estimates for the calculated mass splittings, we need to assess the relevant systematic uncertainties associated with our procedures. In total, we consider uncertainties arising from correlator fits, from the charm-quark mass tuning procedure, from the correction to physical charm-quark mass described in Sec. III C, from the chiralcontinuum fit, and from our limited knowledge of the lattice scale. For the charm-quark tuning procedure, the main uncertainty in the continuum limit arises from the effect of the lattice-scale determination on the charm-quark tuning. We account for this uncertainty as part of our scale-setting uncertainty below. All relevant uncertainties are tabulated in Table VII. We now discuss them in turn.

Variations of the correlator fits
We have investigated many variations of the correlator fits, including correlator basis variations, variations of the fit interval, fit shape, etc., and found that our results are stable under sensible variations of the fitting procedure. Our final choices, which are displayed in full in Tables XI and XII (in Appendix B), are reasonably conservative and encompass almost all stable fit choices with a reasonable goodness of fit. For our final uncertainty estimate, we therefore do not include an additional uncertainty for these variations in correlator fits.

Uncertainty from the determination of the slope in κ c
Variation "G" from Table VI illustrates the results when the slope in κ c used for the charm-quark mass corrections is varied by one standard deviation. This variation is small and straightforward to quantify.

Variations of the chiral-continuum fits
Beyond our default fit results, Table VI lists several variations of the chiral-continuum fit we performed. The results associated with these variations are shown together with our default results in Fig. 7. Among these, variations A-D vary the fit forms used, while variations E and F test whether the results are sensitive to the prior widths selected for the coefficients of the heavy-quark-discretization shapes from Fig. 1. In general the variations among the different fits are rather mild. Significant variations have been discussed for each observables in the previous subsections.
For the S-wave hyperfine splitting and the P-wave tensor splitting, we assess the systematic uncertainty of the chiral-continuum fit by taking the difference between the default fit and variation "B", which results from omitting our data at the coarsest lattice spacing. For all other splittings, the variations are insignificant compared with the statistical uncertainty of the fit. While wider priors leads to a slightly increased uncertainty estimate, there is no significant variation in our best estimates for the splittings.

Uncertainty from the determination of the lattice scale
The procedure for our determination of the scale-setting uncertainty is described above in Sec. III E. For the 1S-hyperfine and 1P-1S splittings this uncertainty is of the same size as the statistical uncertainty from the chiral-continuum fit. In particular, the indirect uncertainty stemming from the uncertainty of the determination of the charm-quark hopping parameter κ c on the scale setting is quite large and this uncertainty has been neglected in some of the  Table VI. literature. For the spin-orbit splitting the uncertainty from scale setting is somewhat smaller than the statistical uncertainty after extrapolation to the physical point. The scale-setting uncertainties for the other splittings are small.

F. Comparison with previous calculations
The Fermilab Lattice and MILC collaborations have previously reported results for the mass splittings in the low-lying charmonium spectrum [2]. Our current results use the same library of gauge configurations. Compared with the previous study we make use of finer lattice spacings, a better determination of the physical quark masses (in particular an improved determination of the charm-quark hopping parameter κ c [45]) and of the lattice spacings used in the simulation [60]. All these ingredients allowed us to perform a more sophisticated chiral-continuum extrapolation. The new results supersede those of Ref. [2].  1P hyperfine −6.1 ± 4.2 ± 0.1 more precise. For the 1P tensor splitting our more elaborate chiral-continuum extrapolation leads to a significant increase in the estimate of the associated uncertainty; the previously quoted uncertainty was probably underestimated.
For the 1S hyperfine splitting there have been several lattice simulations aimed at a full control of systematic uncertainties in the QCD calculation of the connected contribution to the hyperfine splitting [2,20,21,64,65]. Preliminary results with simulations using the HISQ action for charm quarks have also been presented in Ref. [22]. In particular, all these references quote results for physical sea-quark masses in the continuum limit. Figure 8 shows a visual comparison of these calculations. The results from various collaborations are quite consistent. Unfortunately, all these results neglect effects from annihilation of the valence charm-quarks. These have previously been estimated from lattice QCD [62] and from perturbation theory [64]. Note that these results disagree in the sign of the annihilation effects.
More importantly, charm-anticharm annihilation in the physical system results in a substantial total hadronic width of the η c , dominating the total width of 32.0(8) MeV [25]. The PDG lists many decays both into hadronic resonances and into stable final states [25]. While lattice QCD studies of hadronic resonances using Lüscher's finite volume method [66,67] are continuing to make considerable progress (for a review see Ref. [68]) and are now being applied to states close to double open charm thresholds [11][12][13]16], a rigorous study of the For comparison we also show the 1S hyperfine splitting from the PDG [25]. Note that the lattice calculations neglecting disconnected contributions and treating the η c as stable need not result in the same value as the PDG. η c on the lattice is currently out of reach. 5 In our current calculation, the uncertainty from neglecting disconnected contributions is now the largest uncertainty in the error budget for the 1S hyperfine splitting.

V. CONCLUSIONS AND OUTLOOK
In this paper, we have presented results for the splittings of low-lying charmonium states. Table IX shows a comparison of our results with the experimental values [25]. Within our uncertainty estimates, which are described in detail in Sec. IV E, the lattice QCD postdictions are in excellent agreement with experiment, demonstrating that heavy-quark discretization effects for charmonium are well controlled in our setup. Our results improve upon previous results by the Fermilab Lattice and MILC collaborations presented in Ref. [2], which are now superseded.
While our determination of the 1S hyperfine splitting uses the estimate for the charmannihilation contribution from Ref. [62], all current lattice determination including the results presented here neglect effects from charm-anticharm annihilation. For the 1S hyperfine splitting this is now the largest source of uncertainty. A possible direction of further research in this context would be a precision study and prediction of spin-splittings in the B c system, where the contributions from annihilation diagrams are absent. Note that the the hyperfine splitting between the B * c and B c mesons has already been predicted from lattice QCD in [70][71][72], while the B * c has not yet been seen in experiment.

Appendix B: Ground-state energy levels
Tables XI and XII list the determined ground state masses for each ensemble and quantum number combination along with the interpolator basis, reference timeslice t 0 of the variational method, fit range and fit type used to obtain the result.