How much do charm sea quarks affect the charmonium spectrum?

# How much do charm sea quarks affect the charmonium spectrum?

[    [    [
###### Abstract

The properties of charmonium states are or will be intensively studied by the B-factories Belle II and BESIII, the LHCb and PANDA experiments and at a future Super-- Factory. Precise lattice calculations provide valuable input and several results have been obtained by simulating up, down and strange quarks in the sea. We investigate the impact of a charm quark in the sea on the charmonium spectrum, the renormalization group invariant charm-quark mass and the scalar charm-quark content of charmonium. The latter is obtained by the direct computation of the mass-derivatives of the charmonium masses. We do this investigation in a model, QCD with two degenerate charm quarks. The absence of light quarks allows us to reach very small lattice spacings down to fm. By comparing to pure gauge theory we find that charm quarks in the sea affect the hyperfine splitting at a level below 2%. The most significant effects are 5% in and 3% in the value of the charm quark content of the meson. Given that we simulate two charm quarks these estimates are upper bounds for the contribution of a single charm quark. We show that lattice spacings fm are needed for safe continuum extrapolations of the charmonium spectrum with O() improved Wilson quarks. A useful relation for the projection to the desired parity of operators in two-point functions computed with twisted mass fermions is proven.

ucy,wup] Salvatore Calì wup]Francesco Knechtli wup]Tomasz Korzec

WUB/19-03

How much do charm sea quarks affect the charmonium spectrum?

Department of Physics, University of Cyprus, P.O. Box 20537,1678 Nicosia, Cyprus

Department of Physics, Bergische Universität Wuppertal, Gaussstr. 20, 42119 Wuppertal, Germany

## 1 Introduction

The charmonium system is frequently characterized as the “hydrogen atom” of meson spectroscopy owing to the fact that it is non-relativistic enough to be reasonably well described by certain potential models [1]. It is a perfect testing ground for a comparison of theory with experiment. Over the last years, there has been a renewed interest in spectral calculations with charmonia because of the experimental discovery of many states which are not predicted by potential models [2], e.g. the so-called X,Y,Z states, like the state [3] or the pentaquark candidates [4]. More exciting experimental data are expected from the B-factories Belle II [5] and BESIII [6], the LHCb experiment [7], the PANDA experiment at FAIR [8] and at a future Super-- Factory [9].

Simulations of QCD on the lattice are a first-principle tool for precision computations of charmonium states below the open charm thresholds ( etc.) [10, 11, 12, 13, 14, 15], see also [16]. States above the open charm thresholds decay strongly and multi-hadron channels need to be included for a full treatment. The masses of these resonances can be computed in the approximation that they are treated as stable and are accurate up to the hadronic width [10, 11].

For the computation of the charmonium spectrum the relevant quarks to include in the lattice simulations are , , , and . The question which we address in this work is the necessity to include the charm quark in the sea, i.e. as a dynamical quark which contributes through loops and not only as a valence quark. QCD with dynamical quarks is cheaper to simulate than QCD with dynamical quarks. Adding a dynamical charm quark requires finer lattices than they are needed for the lighter quarks and complicates the tuning of the parameters.

For processes at energies which are much smaller than the charm-quark mass the charm quark decouples [17]. It can be integrated out and its effects are absorbed in the renormalization of the gauge coupling and light quark masses, and in small corrections proportional to inverse powers of . In [18, 19, 20] the decoupling of the charm quark at low energies was studied in a model, namely QCD with degenerate heavy quarks of mass . Simulations at very small lattice spacings down to in physical volumes comparable to those used in the Yang–Mills theory allow to control the continuum limit. Concerning the renormalization the model study confirmed that treating decoupling in perturbation theory only introduces small non-perturbative corrections which can be estimated through the model calculation. Denoting a low energy scale of mass dimension one by , the mass-scaling function defined by

 ηM = MS∂S∂M, (1.1)

where is the renormalization group invariant mass of the heavy quark, is universal (i.e. it does not depend on the specific scale chosen) up to non-perturbative corrections . In [20] the conclusion was that for the charm quark in QCD. We emphasize that eq. (1.1) corresponds to the charm quark content of the nucleon and is needed to compute the cross-section of the scalar interaction of dark matter with nucleons.

In this work we extend the model study of charm loop effects to observables which explicitly depend on a valence charm quark. The paper is organized as follows. In Section 2 we introduce the model. Section 3 explains our lattice setup based on twisted mass fermions at maximal twist, the observables and the computation of their derivative with respect to the quark mass. In Section 4 we present our results for the charm loop effects, specifically in the charmonium spectrum and the renormalized charm-quark mass. We also compute the generalization of the mass-scaling function in eq. (1.1) to describe the charm-quark mass-dependence of the charmonium states. All our results are evaluated after continuum extrapolation and we discuss the size of lattice artifacts. Section 5 contains the summary of this work. Appendix A shows how to construct two-point functions which project to definite parity states with twisted mass fermions. In Appendix B the charmonium masses obtained on our ensembles are listed.

## 2 Model

Consider QCD with quarks , . We denote their Dirac operators by . Our goal is to estimate the contribution of charm-quark loops in physical observables , where represents the gauge field. The expectation value of the observable is

 ⟨A[qi,U]⟩ = 1Z∫D[U]⎛⎝∏j=u,d,sdetDj⎞⎠detDc~A[D−1i,U]e−S[U] (2.2)

The charm-quark loop effects1 1 1 Notice that here we mean non-perturbative effects due to quark loops on arbitrary gauge backgrounds. stem from the determinant . Quenching the charm, i.e. setting means neglecting the charm loops. This approximation is made in the computations of the charmonium spectrum of Refs.  [10, 11, 12, 13, 14]. In order to assess how good this approximation is, one would need a comparison in the continuum limit with simulations where a dynamical charm quark is added. Assuming this was possible, the comparison would be superfluous since one would stick with the more complete theory anyhow. But adding a dynamical charm quark means a significant increase in the complexity and costs of the simulations. This is so because of the additional tuning of the charm quark mass and the combination of small lattice spacings, which are required by the large charm-quark mass and the large physical volumes, which are needed to accommodate the light mesons. So the really interesting question is if it is possible to decide whether a dynamical charm quark is necessary before doing the simulations.

This is why the study of a model, QCD with just degenerate charm quarks, is appealing. Observables in this model are defined in terms of a doublet of charm quarks and their expectation value is

 ⟨A[qc,U]⟩ = 1Z∫D[U](detDc)2~A[D−1c,U]e−S[U] (2.3) ≡ ⟨~A[D−1c,U]⟩gauge. (2.4)

After matching this theory with a Yang–Mills (or pure gauge) theory, the difference in physical observable will be a direct measure of the effects of charm-quark loops. There are two differences with respect to a comparison between QCD with four (, , , and ) and three (, , ) quarks: in the model we miss the effects of the light quarks and we double the number of sea charm quarks. Since what we are interested in is a comparison of a theory with and without charm quarks in the sea we do not expect the light quarks to affect the difference of the same quantity computed in the two theories much. The extra charm quark in the sea will make the effects larger. For a low energy quantity, where the theory of decoupling applies, the effects scale proportionally to the number of quarks [20], so they are overestimated by a factor of two in the model. For a quantity with a valence charm quark decoupling does not apply in the obvious way2 2 2 Decoupling might apply for differences of masses or for binding energies. and we consider the effects computed with two charm quarks in the sea as an upper bound for those with only one charm quark.

## 3 Simulations

In this section we introduce the lattice setup used for this work and all the observables under investigation. We mainly focus on quantities with an explicit charm-quark dependence, like charmonium masses, the hyperfine splitting and the renormalization group invariant quark mass.

### 3.1 Actions and algorithms

We use relatively simple and theoretically well understood lattice actions for our simulations. For the ensembles the standard Wilson plaquette action [21] is employed. In the case, a doublet of twisted mass Wilson fermions is added [22, 23]. In massless schemes, theories with standard and with twisted-mass fermions share the same renormalization factors, as long as other details of the action are the same. We therefore also include a clover term [24] in our action. Although not necessary for improvement of physical quantities [23] (at maximal twist), it has been shown to reduce artifacts in some cases [25], and more importantly gives us access to the wide range of renormalization factors that have been determined non-perturbatively in the past. In particular we benefit from the knowledge of the critical mass  [26, 27] and the axial current and pseudoscalar density renormalization factors  [28, 29, 30] and  [31, 26].

Since one of our goals is a detailed understanding of charm related lattice artifacts, we simulate also at very fine lattice spacings, much finer than what is currently feasible in simulations that include light quarks. Problems related to deficient sampling of topological sectors are avoided by the implementation of open boundary conditions in the time directions [32]. The spatial dimensions are kept periodic.

To summarize, our action is , with gauge action

 Sg=1g20∑pw(p)tr[1−U(p)], (3.5)

where the summation is over all oriented plaquettes on the lattice, weighted by which is one everywhere except for spatial plaquettes on the temporal boundary time-slices, where it is . is the product of four SU(3) gauge fields around the elementary plaquette . Gauge fields are periodic in spatial directions and absent on temporal links sticking out of the lattice (i.e. open boundaries). The free parameter of the gauge action is the bare coupling . In case of the simulations, a fermionic action is added

 Sf=∑xa4¯χ(x)[Dχ](x), (3.6)

where is a flavor doublet of quarks and the Dirac operator is

 D=Dw+Dsw+m0+iμγ5τ3, (3.7)

with bare mass and twisted bare mass . The third Pauli matrix in the twisted mass term acts in flavor-space, all other terms of the operator are flavor diagonal.

 Dw=3∑μ=012(γμ[∇∗μ+∇μ]−∇∗μ∇μ) (3.8)

is the massless Wilson operator, containing the usual covariant forward and backward finite difference operators and . Finally, the operator in the Sheikholeslami-Wohlert term acts as

 Dswχ(x)=csw3∑μ,ν=0i4σμν^Fμν(x)χ(x). (3.9)

A symmetric discretization of the field strength tensor , as e.g. in [33], is used. The fermionic fields are periodic in spatial directions and satisfy on the first and last time-slice of the lattice. The fermionic part of the action has dimensionless simulation parameters and . The above choice for the actions corresponds to setting the gluonic and fermionic boundary improvement terms to their tree level values.

Both and theories are simulated with a Hybrid Monte Carlo (HMC) [34] algorithm. The molecular dynamics equations are integrated using a fourth order Omelyan-Mryglod-Folk integrator. In the case with fermions a multi-level variant is employed, with fermionic forces being integrated with a coarser step size than the forces deriving from the gauge action. In addition the quark determinant is factorized into two factors which are then represented by two separate path integrals over pseudo fermion fields [35].

The costs are dominated by solutions of the Dirac equation. The relatively high quark masses in our simulations, mean that a standard conjugate gradient algorithm is often more efficient than more complicated preconditioned variants. On the finer lattices however, SAP preconditioning [36] of the equations involving the light Hasenbusch mass is beneficial. Our simulations are carried out using a variant of openQCD [37]. A minor change allows us to choose a different twisted mass parameter in the SAP preconditioner than in the simulation [20, 38]. Table 1 summarizes our simulation parameters.

The simulation algorithm performs very well. In particular no increased critical slowing down due to deficient sampling of topological sectors can be observed. The scaling of the exponential auto-correlation time with the lattice spacing is compatible with the expected behavior [32]. The expected scaling of autocorrelation times with open boundary conditions has been shown in Figure 8 of Ref. [20].

### 3.2 Observables

#### 3.2.1 t0

The Wilson-flow equation [39, 40]

 ∂Vμ(x,t)∂t=−g20(∂x,μSg[V])Vμ(x,t),Vμ(x,0)=Uμ(x), (3.10)

relates a “smeared” gauge field at flow time to the original gauge field , that is integrated over in the path integral. is a gauge action of the smeared fields, in our case the Wilson plaquette action, and the link differential operator is defined in the usual way [41, 39]. It has been shown that correlators constructed from gauge fields at are automatically renormalized [42]. Among other things, this allows to define the low-energy length scale  [39] as the flow time at which

 t2⟨E(t)⟩=0.3. (3.11)

In this equation denotes the Yang-Mills action density at flow-time , away from the temporal boundaries. A different discretization than the one used in the simulations may be used. We follow [33] and use a symmetrized clover definition

 E(x,t) = 14GaμνGaμν, (3.12)

where are the Lie algebra components of the lattice field strength tensor.

#### 3.2.2 Isovector meson masses

We study mesons that are ground states in the channels that are excited by operators . Twisted mass fermions at maximal twist, and , are related to the fields in the physical basis by

 ψ = 1+iγ5τ3√2χ, (3.13) ¯ψ = ¯χ1+iγ5τ3√2. (3.14)

This means that some operators take an unusual form. For flavor components and , the relations are summarized in Table 2.

Meson masses can be extracted from zero momentum correlation functions of the form

 fO1O2(x0,y0)=a6∑x,y⟨O1(x)O†2(y)⟩, (3.15)

with various choices of the operators . We work with definite flavor assignments, e.g. . Then, integrating over the fermions leaves us with a single connected diagram of the form , where are matrices related to the operators in the correlation function, and () is the inverse Dirac operator with positive (negative) twisted mass term. Spatial translation invariance could be exploited to eliminate one of the sums, which would allow to compute a correlator at the cost of 12 solutions of the Dirac equation per choice of . The signal however is highly improved, by keeping the two sums. The trace can then be efficiently estimated stochastically. We use time-dilution with 16 noise sources per time-slice, which amounts to 16 inversions per value and Dirac structure.

An improved signal and exact symmetries are achieved by defining the averages

 fP(x0−a) ≡ 12(fPP(x0,a)+fPP(T−x0,T−a)), (3.16) fA(x0−a) ≡ 12(fPA(x0,a)−fPA(T−x0,T−a)), (3.17) fV(x0−a) ≡ 163∑k=1(fVkVk(x0,a)+fVkVk(T−x0,T−a)), (3.18) fS(x0−a) ≡ 12(fSS(x0,a)+fSS(T−x0,T−a)), (3.19) fT(x0−a) ≡ 16∑j>i(fTijTij(x0,a)+fTijTij(T−x0,T−a)). (3.20)

Enforcing the continuum time reflection symmetries prevents opposite parity operators from mixing, as explained in Appendix A. From the exponential decay of these correlators at , meson masses are extracted. First, effective masses are computed

 ameff(x0+a/2)≡ln(f(x0)f(x0+a)), (3.21)

and the meson mass is then given as a weighted plateau average

 m=thigh∑x0=tloww(x0+a/2)meff(x0+a/2)thigh∑x0=tloww(x0+a/2). (3.22)

The start of the plateau, , is chosen such that excited state contributions are completely negligible, and the weights are given by the inverse squared errors of the corresponding effective masses. All masses that we extract are those of iso-vector mesons. In the light sector these would be called pions or kaons (), - or -mesons (), , , (), , (). However, since both our quarks have the mass of a charm-quark, the meson masses that we obtain are more comparable to the charmonia masses , , respectively. The difference being, that these are iso-scalars and the determination of their masses would require the computation of disconnected (charm annihilation) diagrams.

#### 3.2.3 PCAC mass

Partial conservation of the axial current is an operator relation

 ∂μ^Aμ=2mPCAC^P. (3.23)

On the lattice it holds up to lattice artifacts, when inserted into any correlation function, as long as and are at a different positions than all other operators in the correlator. These lattice artifacts depend on the exact choice of correlation function, and can be quite large. We extract the bare PCAC quark mass from

 mPCAC=~∂0fA2fP, (3.24)

where denotes the symmetric finite difference operator. The lattice artifacts in this quantity increase, when the correlators are evaluated close to the boundary (small ). We form an average value from the time-slices in the plateau region away from boundaries.

For us the main use of the PCAC mass is to find the critical value of the bare mass , i.e. the maximal twist condition. It is given by the value at which . Instead of determining it ourselves, we use very precise critical masses obtained in [26, 27]. These were computed from slightly different correlation functions in a finite volume, and differ from ours by an lattice artifact. We thus do not expect the PCAC masses that we determine to be zero, but to be small and to vanish when the continuum limit is approached. By computing them, we put this expectation to a test. Figure 1 demonstrates that indeed, up to lattice artifacts, we are at maximal twist. If we consider the usual definition of the twist angle

 ω=arctan(μZAmPCAC), (3.25)

the largest deviation from maximal twist () that we encounter in our simulations is around in the ensemble E, whilst the smallest deviation is around in the ensemble W.

#### 3.2.4 RGI quark mass

At maximal twist a renormalized quark mass is given by , and depends on the scale and scheme in which was computed. Away from maximal the more general relation

 ¯¯¯¯¯m=Z−1P√μ2+Z2Am2PCAC (3.26)

holds. We neglect the (very small) contribution due to non-vanishing in our determination, after veryfying that it is compatible with being of . The axial current renormalization factor is scale independent. It has been determined non-perturbatively in the theory with our action, by exploiting current algebra relations in a massless Schrödinger functional [28]. The same technique has also been applied to the theory [29]. In this case also a more precise determination based on universal relations between correlators in a chirally rotated Schrödinger functional exists [30], and these are the values that we use here. The pseudoscalar renormalization factor depends on the renormalization scheme and scale. It is known non-perturbatively in the SF scheme in both  [31] and  [26] theories for a wide range of bare couplings, albeit at slightly different scales. The renormalized charm quark mass can thus be computed in the continuum limit, in this particular scheme. To be able to compare the two theories, also the scales should match. We go one step further and compute directly the RGI masses, which are scale and scheme independent. The necessary relations between renormalized and RGI masses are well known for the scales and schemes used above, namely in the theory [31], and in the theory with two dynamical quarks [26].

#### 3.2.5 Twisted mass derivatives

We also computed the derivatives of all the observables above, with respect to the twisted mass parameter . The twisted mass derivative of a primary observable is given by

 (3.27)

Most quantities we are interested in, are non-linear functions of various primary observables (e.g. , which depends on the correlator at various distances in the plateau region). For these the chain rule dictates

 df(⟨A1⟩,…,⟨AN⟩,μ)dμ=∂f∂μ+N∑i=1∂f∂⟨Ai⟩d⟨Ai⟩dμ. (3.28)

None of the observables that we consider have an explicit dependence, so the last term in eq. (3.27) is absent. The derivative of the action is , and this is all that is needed to compute the twisted-mass derivatives of purely gluonic observables. More precisely,

 ⟨dSdμA[U]⟩ = ia4∑x⟨(¯c1(x)γ5c1(x)−¯c2(x)γ5c2(x))A[U]⟩ (3.29) = ia4∑x⟨tr[γ5(S1(x,x)−S2(x,x))]A[U]⟩gauge (3.30) = −2μa8∑x,y⟨tr[S†1(x,y)S2(x,y)]A[U]⟩gauge. (3.31)

The last line is a consequence of the twisted-mass relation and allows for a more precise stochastic determination of the trace [43]. We found that 64 noise vectors are enough for the errors in the determination of the derivative to be dominated by gauge-noise, rather than the noise from the stochastic trace evaluation.

If the observables depend on fermionic fields too, the first term of eq. (3.27) gives rise to new contractions that have to be computed. These are different for every fermionic observable. In the case of our two-point functions eq. (3.15) we find contractions of the form , that can be immediately computed because both traces have already been estimated for the evaluation of the correlator and of respectively, and new terms

 ia10∑x,y,z⟨tr[γ5S2(z,y)Γ2S1(y,x)Γ1S2(x,z)]−tr[γ5S1(z,x)Γ1S2(x,y)Γ2S1(y,z)]⟩gauge (3.32)

that require some attention. When evaluated stochastically together with the correlator itself, the number of necessary inversions is increased by a factor of 3. While the terms quantify the dependence on the sea-quarks, this last term gives the valance quark mass dependence of the correlator, which is generally much stronger - especially with heavy quarks.

#### 3.2.6 Mass scaling functions

At last, we also investigate the mass scaling functions

 ηx≡μmxdmxdμ=MmxdmxdM, (3.33)

where denotes the mass of a meson in a generic channel (scalar, pseudoscalar, vector, axial vector, tensor) and is the renormalization group invariant quark mass. Note that is a renormalized quantity and its continuum limit can be easily extracted from the measurements performed at different lattice spacings, without the need of any renormalization factor. Notice that by the Hellman-Feynman-Theorem [44], is proportional to the scalar charm quark density between meson states , i.e. the -term

 ηx=1mx⟨x|Mc(¯cc)RGI|x⟩. (3.34)

Once the twisted mass derivatives of the meson correlators are known, the determination of amounts to the evaluation of eq. (3.28) with a particular function . Since the action of QCD does not depend on , the calculation is greatly simplified in this case. Eq. (3.27) receives a single contribution of the form . In the theory on the other hand, also the -derivative of the action must be taken into account.

### 3.3 Parameters, tuning and mis-tuning corrections

Apart from the lattice size, the bare parameters of the simulations are the inverse bare coupling , the bare mass and the bare twisted mass . The choice of corresponds to a choice of the lattice spacing. We choose to simulate at which spans a wide range of lattice spacings, see Table 5 and allows for very controlled continuum extrapolations.

The bare mass is set to its critical value . To achieve this, the values in [26] are fitted to a Padé function, as described in [20]. This puts us to maximal twist, up to . In this situation the physical quark mass is given by the twisted mass parameter . On our finest lattice at we choose

 aμ=McΛ¯¯¯¯¯¯¯MS×ZP(L−11)×¯¯¯¯¯mc(L1)Mc×Λ¯¯¯¯¯¯¯MSL1×aL1, (3.35)

where the ratio of the RGI charm quark mass and the two flavor parameter is set to 4.87, the pseudoscalar renormalization factor at scale and in the SF scheme is  [26], the relation between a renormalized quark mass in the SF-scheme at scale and the RGI quark mass is known in the continuum  [26], and  [45]. Finally in lattice units is obtained by an interpolation of vs data from [26] to . A quadratic fit of as a function of , describes the data very well and yields . The quite substantial errors mean, that our simulated mass corresponds to the charm quark mass only up to about . This is however fully sufficient for us, as long as the relative mass differences between the different ensembles are under better control. To achieve this, we do not use eq. (3.35) at the other lattice spacings. Instead we proceed as follows: the dimensionless, renormalized quantity

 √t0mP=1.807463 (3.36)

is determined on the ensemble with the finest lattice spacing. On the coarser lattices is tuned such that the same value of is obtained. This condition determines the bare twisted mass parameter very precisely and ensures that all ensembles have the same renormalized quark mass up to . Finally, the clover coefficient is set to its non-perturbatively determined value [46].

The tuning of the twisted mass parameter can be only carried out to a limited precision - at most to within the statistical errors. To account for the mis-tuning, a correction is applied to all observables, based on the computed twisted mass derivatives. First a target tuning point is determined

 μ⋆=μ+(√t0mP−1.807463)(d√t0mPdμ)−1, (3.37)

and afterwards all quantities, denoted by below, are corrected

 Φ(μ⋆)=Φ(μ)+(μ⋆−μ)dΦdμ. (3.38)

The error of the tuning point is propagated to the value of taking all correlations into account. It is assumed that the initial tuning was precise enough for the omitted quadratic terms to be negligible, compared to the statistical precision. Figure 2 demonstrates the procedure.

A comparison with direct simulations indicates that even for large shifts of in the linear approximation works well. The true shifts, that are needed are all much smaller, at most . Note that the -shifts could also be computed using the mass reweighting, as explained in [47, 48].

The simulations are carried out at . The valence quarks have  [49], non-perturbative from [49] and three values of the twisted mass parameter, chosen such that a short interpolation to the value of given in eq. (3.36) can be performed. An example of this procedure is shown in Figure 3. Since decoupling applies to , the condition eq. (3.36) means that the quark mass in the theory is the same as in the theory with two flavors, up to and tiny power corrections [19, 20].

### 3.4 Data Analysis

We use the -method [50] for the determination of statistical uncertainties. Observables like the effective mass Eq. (3.21) are non-linear functions of “primary observables”, and their errors are determined as described in [51]. When incorporating the mis-tuning corrections of Section 3.3 the necessary nonlinear functions can become quite unwieldy. For instance, the vector meson mass at depends on the vector correlator in the plateau region, but also on the pseudoscalar correlator in its plateau region, to determine how big a shift in is required. Furthermore, the vector mass depends on the -derivatives of these correlators, on the -derivative of the action and on the -derivative of the action times the correlators. Combinations like depend on even more primary data.

## 4 Results

### 4.1 Raw results

We measured all observables described in the previous section on all ensembles, except of and which were measured only on a subset and the mass derivatives, which were not measured on the ensemble. A somewhat delicate issue is the proper choice of the plateau regions over which the effective masses are averaged. The leading correction to a constant effective mass is given by

 ameff(x0+a/2)=am+c e−Δ1x0+O(e−2Δ1x0)+O(e−Δ2x0), (4.39)

where () is the distance between and the first (second) excited state. In a first preliminary fit we determine and . We are then in the position to choose the plateau region such, that the influence of the excited states on the plateau average eq. (3.22) is negligible compared to its statistical uncertainty. The thus determined plateau regions are collected in Table 3. Figure 4 demonstrates the procedure for the case of ensemble . The effective masses in the axial-vector channel become too noisy, before a clean plateau is reached and are hence excluded from the tables.

The results for the plateau averages are summarized in Table 6 in Appendix B, which shows the results at the simulated parameters, as well as the values corrected for small mis-tunings in the twisted mass parameter.

### 4.2 Continuum extrapolations

We perform continuum extrapolations of dimensionless quantities. These are either ratios of meson masses, namely , and , or the mass-scaling functions and . One last quantity is the renormalized quark mass. We take the continuum limit of the dimensionless ratio of and . All fits are restricted to a region where the data can be well described by the expected leading scaling violations of order . This means, neglecting data with lattice spacings coarser than .

Figure 5-Figure 7 and Table 4 summarize our findings. The data entering the fits are collected in Table 7.

The results in the continuum limit are collected in Table 4

### 4.3 Dynamical charm effects

The comparison of continuum results in the theory with those in the theory directly quantifies the typical size of the effects, that the inclusion of dynamical charm quarks have on observables with valence charm quarks.

Although they were determined very precisely, no significant effect can be seen in the meson mass spectrum. The most significant deviations of around are found in the ratios and . The relative differences between the central values of the first ratio are only %. For the hyperfine splitting this means a charm quark effect of around 2%. In the ratio the central values deviate by .

A clearer difference between the and theories can be observed in the mass-scaling functions and in the RGI quark mass. The values of and the quark mass differ by almost . The relative differences are and . An even larger (but less significant) difference is found in .

### 4.4 Lattice Artifacts

Having access to very fine lattice spacings is crucial for reliable continuum extrapolations. Although our fermionic action, i.e. twisted mass fermions with an additional clover term, is known to have relatively mild lattice artifacts, the continuum value of e.g. would be significantly underestimated if we had access only to our two coarsest lattices (E and N). The finer of the two has a lattice spacing of fm, which is comparable to the finest lattice spacings typically achievable in large-volume simulations with light quarks. The situation is depicted in Figure 8.

The presence of large lattice artifacts of not only affects observables like , but also the value of the lattice spacing itself. Since it is obtained by determining some hadronic length scale in lattice units at finite lattice spacing and dividing it by the continuum value in fm, i.e. , its value depends on the lattice artifacts present in . In our case one possibility to compute the lattice spacings is through the scale fm. Its values in lattice units are known for our bare couplings and the resulting lattice spacings are between fm on ensemble and fm on ensemble . Alternatively, one could determine the lattice spacing through fm [20]. While the two lattice spacing determinations agree well on the fine ensembles, the difference is quite substantial on the coarsest one, where we find fm, i.e. we observe a 37% lattice artifact in ! Since is also determined on the quenched ensembles, we can determine their lattice spacings using the decoupling relation . Note that lattice spacings determined by using the theory as an effective theory for our massive two flavor theory differ from those determined by using it as an (uncontrolled) approximation to full QCD. In particular these lattice spacings depend on the value of in the fundamental theory. Table 5 summarizes our scale setting.

## 5 Conclusions

In this work we presented a determination of the effects of charm quarks in the sea based on a simulation of a model, QCD with charm quarks. By comparing to the pure gauge theory at the matching point defined in eq. (3.36) we can compute the size of these effects. We find that they are below 2% for the hyperfine splitting of charmonium. These are good news for lattice QCD computations of charmonium based on simulations of light quarks in the sea. We also demonstrate in figure 8 that lattice spacings fm are needed for safe continuum extrapolations of the charmonium spectrum when using O() improved Wilson quarks.

We also computed the effects of sea charm quarks in the mass-scaling function of the charmonium masses eq. (3.33) and in the renormalization group invariant charm-quark mass . Table 4 lists the comparison in the continuum limit of these quantities in the and theory. The effects of the charm sea quarks are clearly resolved and their size is 3% for and 5% for . We notice that our results are upper bounds for the effects of a charm sea quark in QCD since in our model we have doubled their number.

Further analysis to compute charm loop effects in decay constants and finestructure of mesons is in progress. So far the disconnected contributions due to charm annihilation [52] have been neglected since we computed isovector charmonium masses in our model. Work on these contributions is under way.

## Acknowledgments

We thank our colleagues in the ALPHA collaboration for access to data analysis tools. We gratefully acknowledge the computer resources granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JUROPA at Jülich Supercomputing Centre (JSC) and by the Gauss Centre for Supercomputing (GCS) through the NIC on the GCS share of the supercomputer JUQUEEN at JSC, with funding by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF). S.C. acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 642069.

## Appendix A Parity and time-reflection symmetries

Following transformations can be considered as a change of variables in the lattice path integral:

Parity

 P:⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩U0(x0,→x)→U0(x0,−→x)Uk(x0,→x)→U†k(x0,−→x−a^k),k=1,2,3χ(x0,→x)→γ0χ(x0,−→x)¯χ(x0,→x)→¯χ(x0,−→x)γ0 (A.40)

Time reflection

 T:⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩U0(x0,→x)→U†0(T−x0−a,→x)Uk(x0,→x)→Uk(T−x0,→x),k=1,2,3χ(x0,→x)→γ0γ5χ(T−x0,→x)¯χ(x0,→x)→¯χ(T−x0,→x)γ5γ0 (A.41)

They are symmetries of the twisted mass action only if . In general these transformations lead to relations between expectation values in theories with positive and with negative twisted masses. E.g. for the two-point functions like eq. (3.15) one finds

 ⟨O1O2⟩ = ⟨P[O1]P[O2]⟩−μ (A.42) ⟨O1O2⟩ = ⟨T[O1]T[O2]⟩−μ. (A.43)

With standard Wilson fermions (), these equations can be used to show that if the operators and have opposite parity, i.e. if . As a consequence, an operator with a definite parity will only excite states with the same parity. This property is lost in the twisted mass formulation, and in general operators will excite states with both parities.

The combined transformation is a symmetry of the twisted mass action

 ⟨O1O2⟩=⟨TP[O1]TP[O2]⟩ (A.44)

For the averaged correlators eq. (3.16)-eq. (3.20) this means

 12⟨O1O2+T[O1]T[O2]⟩=12⟨TP[O1]TP[O2]+P[O1]P[O2]⟩. (A.45)

It is now easy to see that the averaged correlator vanishes, if . So, by enforcing the continuum time reflection symmetry of the correlator, automatically the mixing of opposite parity operators is prohibited.

## References

• [1] E. Eichten, K. Gottfried, T. Kinoshita, J. B. Kogut, K. D. Lane and T.-M. Yan Phys. Rev. Lett. 34 (1975) 369–372. [Erratum: Phys. Rev. Lett.36,1276(1976)].
• [2] S. L. Olsen, XYZ Meson Spectroscopy, in Proceedings, 53rd International Winter Meeting on Nuclear Physics (Bormio 2015): Bormio, Italy, January 26-30, 2015, 2015.
• [3] S. K. Choi et. al. Phys. Rev. Lett. 91 (2003) 262001.
• [4] R. Aaij et. al. Phys. Rev. Lett. 115 (2015) 072001.
• [5] W. Altmannshofer et. al., The Belle II Physics Book, (2018)
• [6] M. Ablikim et. al. Phys. Rev. Lett. 118 (2017), no. 9 092002.
• [7] R. Aaij et. al., Physics case for an LHCb Upgrade II - Opportunities in flavour physics, and beyond, in the HL-LHC era, (2018)
• [8] G. Barucca et. al. Eur. Phys. J. A55 (2019), no. 3 42.
• [9] S. Eidelman Nucl. Part. Phys. Proc. 260 (2015) 238–241.
• [10] L. Liu, G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, P. Vilaseca, J. J. Dudek, R. G. Edwards, B. Joo and D. G. Richards JHEP 07 (2012) 126.
• [11] G. K. C. Cheung, C. O’Hara, G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas and D. Tims JHEP 12 (2016) 089.
• [12] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier and K. Wong Phys. Rev. D75 (2007) 054502.
• [13] C. DeTar, A. S. Kronfeld, S.-h. Lee, D. Mohler and J. N. Simone Phys. Rev. D99 (2019), no. 3 034509.
• [14] M. Padmanath, S. Collins, D. Mohler, S. Piemonte, S. Prelovsek, A. Schäfer and S. Weishaeupl Phys. Rev. D99 (2019), no. 1 014513.
• [15] M. Kalinowski and M. Wagner Phys. Rev. D92 (2015), no. 9 094508.
• [16] F. Knechtli EPJ Web Conf. 202 (2019) 01006.
• [17] S. Weinberg Phys. Lett. B91 (1980) 51–55.
• [18] M. Bruno, J. Finkenrath, F. Knechtli, B. Leder and R. Sommer Phys. Rev. Lett. 114 (2015), no. 10 102001.
• [19] F. Knechtli, T. Korzec, B. Leder and G. Moir Phys. Lett. B774 (2017) 649–655.
• [20] A. Athenodorou, J. Finkenrath, F. Knechtli, T. Korzec, B. Leder, M. K. Marinković and R. Sommer Nucl. Phys. B943 (2019) 114612.
• [21] K. G. Wilson Phys. Rev. D10 (1974) 2445–2459. [,319(1974)].
• [22] R. Frezzotti, P. A. Grassi, S. Sint and P. Weisz JHEP 08 (2001) 058.
• [23] R. Frezzotti and G. C. Rossi JHEP 08 (2004) 007.
• [24] B. Sheikholeslami and R. Wohlert Nucl. Phys. B259 (1985) 572.
• [25] P. Dimopoulos, H. Simma and A. Vladikas JHEP 07 (2009) 007.
• [26] P. Fritzsch, F. Knechtli, B. Leder, M. Marinkovic, S. Schaefer, R. Sommer and F. Virotta Nucl. Phys. B865 (2012) 397–429.
• [27] P. Fritzsch, N. Garron and J. Heitger JHEP 01 (2016) 093.
• [28] M. Lüscher, S. Sint, R. Sommer and H. Wittig Nucl. Phys. B491 (1997) 344–364.
• [29] M. Della Morte, R. Hoffmann, F. Knechtli, R. Sommer and U. Wolff JHEP 07 (2005) 007.
• [30] M. Dalla Brida, T. Korzec, S. Sint and P. Vilaseca Eur. Phys. J. C79 (2019), no. 1 23.
• [31] A. Jüttner, Precision lattice computations in the heavy quark sector. PhD thesis, Humboldt U., Berlin, 2004.
• [32] M. Lüscher and S. Schaefer JHEP 07 (2011) 036.
• [33] M. Lüscher, Advanced lattice QCD, in Probing the standard model of particle interactions. Proceedings, Summer School in Theoretical Physics, NATO Advanced Study Institute, 68th session, Les Houches, France, July 28-September 5, 1997. Pt. 1, 2, pp. 229–280, 1998.
• [34] S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth Phys. Lett. B195 (1987) 216–222.
• [35] M. Hasenbusch Phys. Lett. B519 (2001) 177–182.
• [36] M. Lüscher Comput. Phys. Commun. 156 (2004) 209–220.
• [37] M. Lüscher and S. Schaefer Comput. Phys. Commun. 184 (2013) 519–528.
• [38] C. Alexandrou, S. Bacchio, J. Finkenrath, A. Frommer, K. Kahl and M. Rottmann Phys. Rev. D94 (2016), no. 11 114509.
• [39] M. Lüscher JHEP 08 (2010) 071. [Erratum: JHEP03,092(2014)].
• [40] R. Narayanan and H. Neuberger JHEP 03 (2006) 064.
• [41] M. Lüscher Commun. Math. Phys. 293 (2010) 899–919.
• [42] M. Lüscher and P. Weisz JHEP 02 (2011) 051.
• [43] K. Jansen, C. Michael and C. Urbach Eur. Phys. J. C58 (2008) 261–269.
• [44] R. P. Feynman Phys. Rev. 56 (1939) 340–343.
• [45] M. Della Morte, R. Frezzotti, J. Heitger, J. Rolf, R. Sommer and U. Wolff Nucl. Phys. B713 (2005) 378–406.
• [46] K. Jansen and R. Sommer Nucl. Phys. B530 (1998) 185–203. [Erratum: Nucl. Phys.B643,517(2002)].
• [47] J. Finkenrath, F. Knechtli and B. Leder Nucl. Phys. B877 (2013) 441–456. [Erratum: Nucl. Phys.B880,574(2014)].
• [48] B. Leder and J. Finkenrath PoS LATTICE2014 (2015) 040.
• [49] M. Lüscher, S. Sint, R. Sommer, P. Weisz and U. Wolff Nucl. Phys. B491 (1997) 323–343.
• [50] U. Wolff Comput. Phys. Commun. 156 (2004) 143–153. [Erratum: Comput. Phys. Commun.176,383(2007)].
• [51] S. Schaefer, R. Sommer and F. Virotta Nucl. Phys. B845 (2011) 93–119.
• [52] L. Levkova and C. DeTar Phys. Rev. D83 (2011) 074504.