How much do charm sea quarks affect the charmonium spectrum?

TL;DR

This study assesses the influence of charm sea quarks on the charmonium spectrum and related properties using lattice QCD with two degenerate charm quarks, finding effects below 5% and emphasizing the need for fine lattice spacings.

Contribution

It provides the first detailed lattice QCD analysis of charm sea quark effects on charmonium, including mass shifts and content, in a model with two degenerate charm quarks.

Findings

Charm sea quarks affect hyperfine splitting below 2%.

Charm quark mass is influenced by about 5%.

Charm content of eta_c changes by approximately 3%.

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-c-tau 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 0.023 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 the quark mass and 3% in the value of the charm quark content of the eta_c 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 <0.06 fm are needed for safe continuum extrapolations of the charmonium spectrum with O(a) 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.

Tables

Table 1. Table 1 : Simulation parameters of our ensembles. The columns show the lattice sizes, the gauge coupling β = 6 / g 0 2 𝛽 6 superscript subscript 𝑔 0 2 \beta=6/{g_{0}^{2}} , the critical hopping parameter, the twisted mass parameter μ 𝜇 \mu , the pseudoscalar mass in t 0 subscript 𝑡 0 t_{0} units, the hadronic scale t 0 / a 2 subscript 𝑡 0 superscript 𝑎 2 t_{0}/a^{2} defined in [ 39 ] and the total statistics in molecular dynamics units. Note that even though the number of sites in the temporal direction is even, the temporal extent T 𝑇 T is an odd multiple of a 𝑎 a due to the open boundaries. The links pointing out of the lattice volume are absent.

$N_f$	ID	$\frac{T}{a}\times {\left(\frac{L}{a}\right)}^3$	$\beta$	$\kappa$	$a\mu$	$\sqrt{t_0}{m}_P$	$t_0/a^2$	MDUs
2	E	$95\times {24}^3$	5.300	0.135943	0.36151	1.79303(55)	1.23907(82)	8000
	N	$119\times {32}^3$	5.500	0.136638	0.165997	1.8048(15)	4.4730(93)	8000
	O	$191\times {48}^3$	5.600	0.136710	0.130949	1.7656(14)	6.561(12)	8000
	P	$119\times {32}^3$	5.700	0.136698	0.113200	1.7931(28)	9.105(35)	17184
	S	$191\times {48}^3$	5.880	0.136509	0.087626	1.8130(29)	15.621(60)	23088
	W	$191\times {48}^3$	6.000	0.136335	0.072557	1.8075(43)	22.39(12)	22400
0	qN	$119\times {32}^3$	6.100	–	–	–	4.4329(38)	64000
	qP	$119\times {32}^3$	6.340	–	–	–	9.037(30)	20080
	qW	$191\times {48}^3$	6.672	–	–	–	21.925(83)	73920
	qX	$191\times {64}^3$	6.900	–	–	–	39.41(14)	160200

Table 2. Table 2 : Typical interpolators for meson states and relations between physical and twisted basis. The particle name is the closest relative in nature.

State	$J^{PC}$	Particle	Physical basis	Twisted basis
Scalar	$0^{++}$	${\chi }_{c0}$	$S^{1,2}=\overline{\psi }\frac{{\tau }^{1,2}}{2}\psi$	$\overline{\chi }\frac{{\tau }^{1,2}}{2}\chi$
Pseudoscalar	$0^{-+}$	${\eta }_c$	$P^{1,2}=\overline{\psi }{\gamma }_5\frac{{\tau }^{1,2}}{2}\psi$	$\overline{\chi }{\gamma }_5\frac{{\tau }^{1,2}}{2}\chi$
Vector	$1^{--}$	$J/\psi$	$V_i^{1,2}=\overline{\psi }{\gamma }_i\frac{{\tau }^{1,2}}{2}\psi$	$\pm \overline{\chi }{\gamma }_i{\gamma }_5\frac{{\tau }^{2,1}}{2}\chi$
Axial vector	$1^{++}$	${\chi }_{c1}$	$A_i^{1,2}=\overline{\psi }{\gamma }_i{\gamma }_5\frac{{\tau }^{1,2}}{2}\psi$	$\pm \overline{\chi }{\gamma }_i\frac{{\tau }^{2,1}}{2}\chi$
Tensor333The notation refers to the $\gamma$-structure of the operator.	$1^{+-}$	$h_c$	$T_{ij}^{1,2}=\overline{\psi }{\gamma }_i{\gamma }_j\frac{{\tau }^{1,2}}{2}\psi$	$\overline{\chi }{\gamma }_i{\gamma }_j\frac{{\tau }^{1,2}}{2}\chi$

Table 3. Table 3 : The meson masses are determined from effective masses in the region t low < x 0 + a / 2 < t high subscript 𝑡 low subscript 𝑥 0 𝑎 2 subscript 𝑡 high t_{\rm low}<x_{0}+a/2<t_{\rm high} . The table shows t low / a − t high / a subscript 𝑡 low 𝑎 subscript 𝑡 high 𝑎 t_{\rm low}/a-t_{\rm high}/a for the different ensembles and channels.

ID	$m_P$	$m_V$	$m_S$	$m_T$
E	21-35	21-35	-	-
N	30-58	30-58	26-46	26-46
O	34-71	34-71	-	-
P	37-71	37-71	25-51	30-51
S	47-101	47-101	-	-
W	55-101	55-101	31-41	31-41
qN	32-58	32-58	26-42	26-42
qP	39-71	39-71	28-51	28-41
qW	60-101	60-101	35-44	30-43
qX	104-173	104-173	-	-

Table 4. Table 4 : Results for various quantities in the continuum limit for both the N f = 0 subscript 𝑁 f 0 {N_{\rm f}}=0 and the N f = 2 subscript 𝑁 f 2 {N_{\rm f}}=2 theory.

Quantity	$N_f=2$	$N_f=0$	sea effects [%]
$m_V/m_P$	1.05405(60)	1.05274(46)	0.124(71)
$m_S/m_P$	1.258(14)	1.224(20)	2.7(1.9)
$m_T/m_P$	1.271(38)	1.321(33)	3.9(4.1)
${\eta }_P$	0.6996(81)	0.67553(42)	3.4(1.1)
${\eta }_V$	0.666(31)	0.6060(13)	9.0(4.2)
$M_c/m_P$	0.4764(74)	0.4528(51)	5.0(1.8)

Table 5. Table 5 : Lattice spacings in physical units on our quenched and dynamical ensembles, determined in two different ways.

Ensemble	$a^{L_1}$ [fm]	$a^{t_0}$ [fm]
E	0.066	0.104
N	0.049	0.054
O	0.042	0.045
P	0.036	0.038
S	0.028	0.029
W	0.023	0.024
qN	-	0.054
qP	-	0.038
qW	-	0.024
qX	-	0.018

Table 6. Table 6 : Meson masses and t 0 subscript 𝑡 0 t_{0} in lattice units. For N f = 2 subscript 𝑁 f 2 {N_{\rm f}}=2 simulations the first line contains the values extrapolated to the tuning point μ ⋆ superscript 𝜇 ⋆ \mu^{\star} and the second line the values at the simulated parameters. For N f = 0 subscript 𝑁 f 0 {N_{\rm f}}=0 ensembles the first line contains the values interpolated to μ ⋆ superscript 𝜇 ⋆ \mu^{\star} , and the following three lines contain the values measured at different valence quark masses.

Ensemble	$a\mu$	$a{m}_P$	$a{m}_V$	$a{m}_T$	$a{m}_S$	$t_0/a^2$
E	0.3809(54)	1.667(16)	1.757(18)	-	-	1.172(25)
	0.36151	1.61079(13)	1.69494(36)	-	-	1.23907(82)
N	0.16647(28)	0.8551(12)	0.8985(12)	1.005(51)	0.977(15)	4.468(12)
	0.166	0.85337(17)	0.89687(34)	1.000(51)	0.975(15)	4.4730(93)
O	0.13714(31)	0.7117(13)	0.7488(15)	-	-	6.445(25)
	0.13095	0.68929(11)	0.72742(26)	-	-	6.561(12)
P	0.11482(32)	0.6001(14)	0.6317(15)	0.707(19)	0.7337(53)	9.070(42)
	0.1132	0.59421(24)	0.62579(40)	0.703(20)	0.7269(53)	9.105(35)
S	0.08717(25)	0.4570(10)	0.4814(12)	-	-	15.641(71)
	0.087626	0.45870(18)	0.48309(34)	-	-	15.621(60)
W	0.072557	0.38200(16)	0.40219(31)	0.481(10)	0.4680(37)	22.39(11)
qN	0.17632(11)	0.85846(37)	0.90106(52)	0.994(33)	1.018(18)	4.4329(38)
	0.16	0.806508(100)	0.85216(48)	0.933(38)	0.972(22)	-
	0.17	0.838522(96)	0.88224(41)	0.970(35)	1.000(20)	-
	0.18	0.870096(93)	0.91204(36)	1.008(32)	1.028(18)	-
qP	0.12235(26)	0.60125(100)	0.6328(11)	0.762(43)	0.710(14)	9.037(30)
	0.11	0.56074(32)	0.59454(81)	0.728(51)	0.669(15)	-
	0.12	0.59373(31)	0.62562(69)	0.756(44)	0.703(14)	-
	0.13	0.62619(30)	0.65634(61)	0.783(38)	0.736(12)	-
qW	0.07798(19)	0.38602(73)	0.40607(100)	0.4960(98)	0.4713(63)	21.925(83)
	0.07	0.35943(14)	0.38103(77)	0.472(11)	0.4465(72)	-
	0.08	0.39294(14)	0.41247(57)	0.5022(94)	0.4775(60)	-
	0.09	0.42573(13)	0.44361(44)	0.5322(79)	0.5086(52)	-
qX	0.05771(13)	0.28792(53)	0.30296(56)	-	-	39.41(14)
	0.056	0.28218(10)	0.29755(20)	-	-	-
	0.058	0.288920(99)	0.30389(19)	-	-	-
	0.06	0.295621(98)	0.31022(19)	-	-	-

Table 7. Table 7 : Ratios of masses and the mass scaling functions. All values are at the tuning point μ ⋆ superscript 𝜇 ⋆ \mu^{\star} .

Ensemble	$m_V/m_P$	$m_S/m_P$	$m_T/m_P$	$\overline{m}/m_P$	${\eta }_P$	${\eta }_V$
E	1.0542(17)	-	-	0.4407(46)	0.664(14)	0.698(33)
N	1.05084(36)	1.143(18)	1.175(60)	0.3756(38)	0.6977(75)	0.652(26)
O	1.0520(15)	-	-	0.3717(38)	0.6984(74)	0.632(28)
P	1.05258(55)	1.2226(86)	1.177(31)	0.3691(38)	0.6974(47)	0.660(14)
S	1.05342(64)	-	-	0.3679(38)	0.7002(65)	0.657(32)
W	1.05285(82)	1.225(12)	1.259(29)	0.3664(37)	-	-
qN	1.04968(42)	1.186(21)	1.158(39)	0.3945(20)	0.65290(20)	0.5856(16)
qP	1.0526(10)	1.182(22)	1.269(71)	0.3939(20)	0.66593(59)	0.5972(29)
qW	1.0523(15)	1.222(17)	1.286(26)	0.3933(20)	0.66936(59)	0.6003(42)
qX	1.05226(42)	-	-	0.3903(20)	0.67342(50)	0.6037(12)

Equations

$${\eta }^{\mathrm{M}}$$

$$⟨A[q^i,U]⟩$$

$$⟨A[q^c,U]⟩$$

$$\equiv$$

$$S_{\mathrm{g}}=\frac{1}{g_0^2}\sum _{p}w(p)\mathrm{tr}\left[1-U(p)\right],$$

$$S_{\mathrm{f}}=\sum _x{a}^4\bar{\chi }(x)[D\chi ](x),$$

$$D=D_{\mathrm{w}}+D_{\mathrm{sw}}+m_0+i\mu {\gamma }_5{\tau }^3,$$

$$D_{\mathrm{w}}=\sum _{\mu =0}^3\frac{1}{2}\left({\gamma }_{\mu }\left[{\nabla }_{\mu }^{\ast }+{\nabla }_{\mu }\right]-{\nabla }_{\mu }^{\ast }{\nabla }_{\mu }\right)$$

$$D_{\mathrm{sw}}\chi (x)=c_{\mathrm{sw}}\sum _{\mu ,\nu =0}^3\frac{i}{4}{\sigma }_{\mu \nu }{\hat{F}}_{\mu \nu }(x)\chi (x).$$

$$\frac{\partial V_{\mu }(x,t)}{\partial t}=-g_0^2\left({\partial }_{x,\mu }{S}_{\mathrm{g}}[V]\right)V_{\mu }(x,t),V_{\mu }(x,0)=U_{\mu }(x),$$

$$t^2⟨E(t)⟩=0.3.$$

$$E(x,t)$$

$$\psi$$

$$\bar{\psi }$$

$$f_{{\mathcal{O}}_1{\mathcal{O}}_2}(x_0,y_0)=a^6\sum _{\mathbf{x},\mathbf{y}}⟨{\mathcal{O}}_1(x){\mathcal{O}}_2^{†}(y)⟩,$$

$$f_P(x_0-a)$$

$$f_A(x_0-a)$$

$$f_V(x_0-a)$$

$$f_S(x_0-a)$$

$$f_T(x_0-a)$$

$$a{m}^{\mathrm{eff}}(x_0+a/2)\equiv \ln \left(\frac{f(x_0)}{f(x_0+a)}\right),$$

$$m=\frac{\sum _{x_0=t_{\mathrm{low}}}^{t_{\mathrm{high}}}w(x_0+a/2)m^{\mathrm{eff}}(x_0+a/2)}{\sum _{x_0=t_{\mathrm{low}}}^{t_{\mathrm{high}}}w(x_0+a/2)}.$$

$${\partial }_{\mu }{\hat{A}}_{\mu }=2m_{\mathrm{PCAC}}\hat{P}.$$

$$m_{\mathrm{PCAC}}=\frac{{\tilde{\partial }}_0{f}_A}{2f_P},$$

$$\omega =\arctan \left(\frac{\mu }{Z_A{m}_{\mathrm{PCAC}}}\right),$$

$$\overline{m}=Z_P^{-1}\sqrt{{\mu }^2+Z_A^2{m}_{\mathrm{PCAC}}^2}$$

$$\frac{\mathrm{d}⟨A⟩}{\mathrm{d}\mu }=-⟨\frac{\mathrm{d}S}{\mathrm{d}\mu }A⟩+⟨\frac{\mathrm{d}S}{\mathrm{d}\mu }⟩⟨A⟩+⟨\frac{\mathrm{d}A}{\mathrm{d}\mu }⟩.$$

$$\frac{\mathrm{d}f(⟨A_1⟩,\dots ,⟨A_N⟩,\mu )}{\mathrm{d}\mu }=\frac{\partial f}{\partial \mu }+\sum _{i=1}^N\frac{\partial f}{\partial ⟨A_i⟩}\frac{\mathrm{d}⟨A_i⟩}{\mathrm{d}\mu }.$$

$$⟨\frac{\mathrm{d}S}{\mathrm{d}\mu }A[U]⟩$$

$$=$$

$$i{a}^{10}\sum _{\mathbf{x},\mathbf{y},z}{⟨\mathrm{tr}\left[{\gamma }_5{\mathcal{S}}_2(z,y){\Gamma }_2{\mathcal{S}}_1(y,x){\Gamma }_1{\mathcal{S}}_2(x,z)\right]-\mathrm{tr}\left[{\gamma }_5{\mathcal{S}}_1(z,x){\Gamma }_1{\mathcal{S}}_2(x,y){\Gamma }_2{\mathcal{S}}_1(y,z)\right]⟩}^{\mathrm{gauge}}$$

$${\eta }_x\equiv \frac{\mu }{m_x}\frac{\mathrm{d}{m}_x}{\mathrm{d}\mu }=\frac{M}{m_x}\frac{\mathrm{d}{m}_x}{\mathrm{d}M},$$

$${\eta }_x=\frac{1}{m_x}⟨x|M_c(\bar{c}c{)}_{\mathrm{RGI}}|x⟩.$$

Full text

1. 1 Introduction

2. 2 Model

3. 3 Simulations

4. 3.1 Actions and algorithms

5. 3.2 Observables

6. 3.2.1 $t_{0}$

7. 3.2.2 Isovector meson masses

8. 3.2.3 PCAC mass

9. 3.2.4 RGI quark mass

10. 3.2.5 Twisted mass derivatives

11. 3.2.6 Mass scaling functions

12. 3.3 Parameters, tuning and mis-tuning corrections

13. 3.4 Data Analysis

14. 4 Results

15. 4.1 Raw results

16. 4.2 Continuum extrapolations

17. 4.3 Dynamical charm effects

18. 4.4 Lattice Artifacts

19. 5 Conclusions

20. A Parity and time-reflection symmetries

21. B Tables

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 $X(3872)$ state [3] or the $P_{c}$ 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-$c$-$\tau$ Factory [9].

Simulations of QCD on the lattice are a first-principle tool for precision computations of charmonium states below the open charm thresholds ($D\bar{D}$ 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 $u$, $d$, $s$, and $c$. The question which we address in this work is the necessity to include the charm quark $c$ in the sea, i.e. as a dynamical quark which contributes through loops and not only as a valence quark. QCD with ${N_{\rm f}}=2+1\;(u,d,s)$ dynamical quarks is cheaper to simulate than QCD with ${N_{\rm f}}=2+1+1\;(u,d,s,c)$ 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 $E$ which are much smaller than the charm-quark mass $M_{\rm c}$ 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 $M_{\rm c}$. In [18, 19, 20] the decoupling of the charm quark at low energies was studied in a model, namely QCD with ${N_{\rm f}}=2$ degenerate heavy quarks of mass $1.2\,M_{\rm c}\gtrsim M\gtrsim M_{\rm c}/8$. Simulations at very small lattice spacings down to $a=0.023\,{\rm fm}$ 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 ${\cal S}$, the mass-scaling function defined by

[TABLE]

where $M$ 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 $1/M^{2}$ corrections $\Delta\eta^{\mathrm{M}}_{\mathrm{NP}}$. In [20] the conclusion was that $\Delta\eta^{\mathrm{M}}_{\mathrm{NP}}<0.014$ 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 $q^{i}$, $i=\{u,d,s,c\}$. We denote their Dirac operators by $D_{i}$. Our goal is to estimate the contribution of charm-quark loops in physical observables $A[q^{i},U]$, where $U$ represents the gauge field. The expectation value of the observable is

[TABLE]

The charm-quark loop effects111 Notice that here we mean non-perturbative effects due to quark loops on arbitrary gauge backgrounds. stem from the determinant $\det D_{c}$. Quenching the charm, i.e. setting $\det D_{c}=1$ 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 ${N_{\rm f}}=2$ degenerate charm quarks, is appealing. Observables in this model are defined in terms of a doublet of charm quarks $q^{c}=(c_{1}~{},c_{2})$ and their expectation value is

[TABLE]

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 ($u$, $d$, $s$, and $c$) and three ($u$, $d$, $s$) 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 way222 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 ${N_{\rm f}}=0$ ensembles the standard Wilson plaquette action [21] is employed. In the ${N_{\rm f}}=2$ 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 $O(a)$ improvement of physical quantities [23] (at maximal twist), it has been shown to reduce $O(a^{2})$ 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 $m_{\rm cr}$ [26, 27] and the axial current and pseudoscalar density renormalization factors $Z_{A}$ [28, 29, 30] and $Z_{P}$ [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 $S=S_{\rm g}+S_{\rm f}$, with gauge action

[TABLE]

where the summation is over all oriented plaquettes $p$ on the lattice, weighted by $w(p)$ which is one everywhere except for spatial plaquettes on the temporal boundary time-slices, where it is $1/2$. $U(p)$ is the product of four SU(3) gauge fields $U_{\mu}(x)$ around the elementary plaquette $p$. 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 $g_{0}^{2}\equiv 6/\beta$. In case of the ${N_{\rm f}}=2$ simulations, a fermionic action is added

[TABLE]

where $\chi=(c_{1},c_{2})^{\top}$ is a flavor doublet of quarks and the Dirac operator is

[TABLE]

with bare mass $m_{0}$ and twisted bare mass $\mu$. The third Pauli matrix $\tau^{3}$ in the twisted mass term acts in flavor-space, all other terms of the operator are flavor diagonal.

[TABLE]

is the massless Wilson operator, containing the usual covariant forward and backward finite difference operators $\nabla_{\mu}\chi(x)=U_{\mu}(x)\chi(x+\hat{\mu})-\chi(x)$ and $\nabla^{*}_{\mu}\chi(x)=\chi(x)-U^{\dagger}_{\mu}(x-\hat{\mu})\chi(x-\hat{\mu})$. Finally, the operator in the Sheikholeslami-Wohlert term acts as

[TABLE]

A symmetric discretization of the field strength tensor $\hat{F}_{\mu\nu}$, as e.g. in [33], is used. The fermionic fields are periodic in spatial directions and satisfy $D\chi=0$ on the first and last time-slice of the lattice. The fermionic part of the action has dimensionless simulation parameters $\kappa\equiv\frac{1}{2am_{0}+8},a\mu$ and $c_{\rm sw}$. The above choice for the actions corresponds to setting the gluonic and fermionic boundary improvement terms to their tree level values.

Both ${N_{\rm f}}=0$ and ${N_{\rm f}}=2$ 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 $\tau_{\rm exp}\propto a^{-2}$ 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 $t_{0}$

The Wilson-flow equation [39, 40]

[TABLE]

relates a “smeared” gauge field $V_{\mu}(x,t)$ at flow time $t$ to the original gauge field $U_{\mu}(x)$, that is integrated over in the path integral. $S_{\rm g}[V]$ is a gauge action of the smeared fields, in our case the Wilson plaquette action, and the link differential operator $\partial_{x,\mu}$ is defined in the usual way [41, 39]. It has been shown that correlators constructed from gauge fields at $t>0$ are automatically renormalized [42]. Among other things, this allows to define the low-energy length scale $t_{0}$ [39] as the flow time $t$ at which

[TABLE]

In this equation $E(t)$ denotes the Yang-Mills action density at flow-time $t$, 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

[TABLE]

where $G_{\mu\nu}^{a}(x,t)$ 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 $\bar{\psi}\Gamma\tau^{a}\psi$. Twisted mass fermions at maximal twist, $\bar{\chi}$ and $\chi$, are related to the fields in the physical basis by

[TABLE]

This means that some operators take an unusual form. For flavor components $\tau^{1}$ and $\tau^{2}$, the relations are summarized in Table 2.

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

[TABLE]

with various choices of the operators $\mathcal{O}_{i}$. We work with definite flavor assignments, e.g. $P^{+}\equiv P^{1}+iP^{2}=\bar{c}_{1}\gamma_{5}c_{2}$. Then, integrating over the fermions leaves us with a single connected diagram of the form $-\sum_{\mathbf{x},\mathbf{y}}\langle{\rm tr}[\Gamma_{1}{\mathcal{S}}_{1}(x,y)\Gamma_{2}{\mathcal{S}}_{2}(y,x)]\rangle^{\rm gauge}$, where $\Gamma_{i}$ are $4\times 4$ matrices related to the operators in the correlation function, and ${\mathcal{S}}_{1}$ (${\mathcal{S}}_{2}$) 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 $y_{0}$. 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 $U(1)$ noise sources per time-slice, which amounts to 16 inversions per $y_{0}$ value and Dirac structure.

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

[TABLE]

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 $x_{0}\gg a$, meson masses are extracted. First, effective masses are computed

[TABLE]

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

[TABLE]

The start of the plateau, $t_{\rm low}$, is chosen such that excited state contributions are completely negligible, and the weights $w$ 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 ($f_{P},f_{A}$), $\rho$- or $K^{*}$-mesons ($f_{V}$), $a_{0}$, $f_{0}$, $K_{0}^{\star}$ ($f_{S}$), $h_{1}$, $b_{1}$ ($f_{T}$). 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 $\eta_{c}$, $J/\psi$,$\ldots$ 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

[TABLE]

On the lattice it holds up to lattice artifacts, when inserted into any correlation function, as long as $A$ and $P$ 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

[TABLE]

where $\tilde{\partial}_{\mu}$ denotes the symmetric finite difference operator. The lattice artifacts in this quantity increase, when the correlators are evaluated close to the boundary (small $x_{0}$). 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 $m_{0}$, i.e. the maximal twist condition. It is given by the value at which $m_{\rm PCAC}=0$. 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 $O(a)$ 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

[TABLE]

the largest deviation from maximal twist ($\omega=\pi/2$) that we encounter in our simulations is around $6\%$ in the ensemble E, whilst the smallest deviation is around $2\%$ in the ensemble W.

3.2.4 RGI quark mass

At maximal twist a renormalized quark mass is given by $\overline{m}=Z_{P}^{-1}\mu$, and depends on the scale and scheme in which $Z_{P}$ was computed. Away from maximal the more general relation

[TABLE]

holds. We neglect the (very small) contribution due to non-vanishing $m_{\rm PCAC}$ in our determination, after veryfying that it is compatible with being of $O(a)$. The axial current renormalization factor $Z_{A}$ is scale independent. It has been determined non-perturbatively in the ${N_{\rm f}}=0$ 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 ${N_{\rm f}}=2$ 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 $Z_{P}$ depends on the renormalization scheme and scale. It is known non-perturbatively in the SF scheme in both ${N_{\rm f}}=0$ [31] and ${N_{\rm f}}=2$ [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 $M/\overline{m}=1.157(12)$ in the ${N_{\rm f}}=0$ theory [31], and $M/\overline{m}=1.308(16)$ 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 $\mu$. The twisted mass derivative of a primary observable $A$ is given by

[TABLE]

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

[TABLE]

None of the observables that we consider have an explicit $\mu$ dependence, so the last term in eq. (3.27) is absent. The derivative of the action is ${\rm d}S/{\rm d}\mu=\sum_{x}\bar{\chi}i\gamma_{5}\tau^{3}\chi$, and this is all that is needed to compute the twisted-mass derivatives of purely gluonic observables. More precisely,

[TABLE]

The last line is a consequence of the twisted-mass relation $D_{1}-D_{2}=2i\gamma_{5}\mu$ and allows for a more precise stochastic determination of the trace [43]. We found that 64 $U(1)$ 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 $-a^{10}\sum_{\mathbf{x},\mathbf{y},z}\langle{\rm tr}[\Gamma_{1}{\mathcal{S}}_{2}(x,y)\Gamma_{2}{\mathcal{S}}_{1}(y,x)]{\rm tr}[\gamma_{5}({\mathcal{S}}_{1}(z,z)-{\mathcal{S}}_{2}(z,z))]\rangle^{\rm gauge}$, that can be immediately computed because both traces have already been estimated for the evaluation of the correlator and of ${\rm d}S/{\rm d}\mu$ respectively, and new terms

[TABLE]

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 ${\rm d}S/{\rm d}\mu$ 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

[TABLE]

where $m_{x}$ denotes the mass of a meson in a generic $x$ channel (scalar, pseudoscalar, vector, axial vector, tensor) and $M$ is the renormalization group invariant quark mass. Note that $\eta_{x}$ 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], $\eta_{x}$ is proportional to the scalar charm quark density between meson states $x$, i.e. the $\sigma$-term

[TABLE]

Once the twisted mass derivatives of the meson correlators are known, the determination of $\eta_{x}$ amounts to the evaluation of eq. (3.28) with a particular function $f$. Since the action of ${N_{\rm f}}=0$ QCD does not depend on $\mu$, the calculation is greatly simplified in this case. Eq. (3.27) receives a single contribution of the form $\langle{\rm d}\tilde{A}[D^{-1},U]/{\rm d}\mu\rangle$. In the ${N_{\rm f}}=2$ theory on the other hand, also the $\mu$-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 ${N_{\rm f}}=2$ simulations are the inverse bare coupling $\beta$, the bare mass $am_{0}$ and the bare twisted mass $a\mu$. The choice of $\beta$ corresponds to a choice of the lattice spacing. We choose to simulate at $\beta\in\{5.3,5.5,5.6,5.7,5.88,6.0\}$ 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 $m_{0}=m_{\rm cr}$. 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 $O(a^{2})$. In this situation the physical quark mass is given by the twisted mass parameter $\overline{m}=Z_{P}^{-1}\mu$. On our finest lattice at $\beta=6.0$ we choose

[TABLE]

where the ratio of the RGI charm quark mass and the two flavor $\Lambda$ parameter is set to 4.87, the pseudoscalar renormalization factor at scale $L_{1}^{-1}$ and $\beta=6.0$ in the SF scheme is $Z_{P}=0.5184(33)$ [26], the relation between a renormalized quark mass in the SF-scheme at scale $L_{1}^{-1}$ and the RGI quark mass $M$ is known in the continuum $M/\overline{m}(L_{1}^{-1})=1.308(16)$ [26], and $\Lambda_{\overline{\rm MS}}L_{1}=0.649(45)$ [45]. Finally $L_{1}$ in lattice units is obtained by an interpolation of $L_{1}/a$ vs $\beta$ data from [26] to $\beta=6.0$. A quadratic fit of ${\rm log}(L_{1}/a)$ as a function of $\beta$, describes the data very well and yields $L_{1}/a_{\beta=6.0}=17.27(70)$. The quite substantial errors mean, that our simulated mass corresponds to the charm quark mass only up to about $10\%$. 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

[TABLE]

is determined on the ensemble with the finest lattice spacing. On the coarser lattices $a\mu$ is tuned such that the same value of $\sqrt{t_{0}}m_{P}$ 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 $O(a^{2})$. Finally, the clover coefficient $c_{\rm sw}$ 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 $\mu^{\star}$ is determined

[TABLE]

and afterwards all quantities, denoted by $\Phi$ below, are corrected

[TABLE]

The error of the tuning point $\mu^{*}$ is propagated to the value of $\Phi(\mu^{*})$ 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 $\approx 15\%$ in $a\mu$ the linear approximation works well. The true shifts, that are needed are all much smaller, at most $5.40\%$. Note that the $\mu$-shifts could also be computed using the mass reweighting, as explained in [47, 48].

The ${N_{\rm f}}=0$ simulations are carried out at $\beta\in\{6.1,6.34,6.672,6.9\}$. The valence quarks have $m_{0}=m_{\rm cr}$ [49], non-perturbative $c_{\rm sw}$ from [49] and three values of the twisted mass parameter, chosen such that a short interpolation to the value of $\sqrt{t_{0}}m_{P}$ given in eq. (3.36) can be performed. An example of this procedure is shown in Figure 3. Since decoupling applies to $t_{0}$, the condition eq. (3.36) means that the quark mass in the ${N_{\rm f}}=0$ theory is the same as in the theory with two flavors, up to $O(a^{2})$ and tiny $O(\Lambda^{2}/M_{c}^{2})$ power corrections [19, 20].

3.4 Data Analysis

We use the $\Gamma$-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 $\mu^{\star}$ 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 $\mu$ is required. Furthermore, the vector mass depends on the $\mu$-derivatives of these correlators, on the $\mu$-derivative of the action and on the $\mu$-derivative of the action times the correlators. Combinations like $\sqrt{t_{0}}m_{V}$ 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 $m_{T}$ and $m_{S}$ which were measured only on a subset and the mass derivatives, which were not measured on the $W$ 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

[TABLE]

where $\Delta_{1}$ ($\Delta_{2}$) is the distance between $m$ and the first (second) excited state. In a first preliminary fit we determine $\Delta_{1}$ and $c$. 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 $W$. 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 $m_{V}/m_{P}$, $m_{S}/m_{P}$ and $m_{T}/m_{P}$, or the mass-scaling functions $\eta_{P}$ and $\eta_{V}$. One last quantity is the renormalized quark mass. We take the continuum limit of the dimensionless ratio of $\overline{m}$ and $m_{P}$. All fits are restricted to a region where the data can be well described by the expected leading scaling violations of order $a^{2}$. This means, neglecting data with lattice spacings coarser than $a^{2}/t_{0}>0.25$.

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 ${N_{\rm f}}=2$ theory with those in the ${N_{\rm f}}=0$ 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 $1.6\sigma$ are found in the ratios $m_{V}/m_{P}$ and $m_{S}/m_{P}$. The relative differences between the central values of the first ratio are only $([m_{V}/m_{P}]^{{N_{\rm f}}=2}-[m_{V}/m_{P}]^{{N_{\rm f}}=0})/[m_{V}/m_{P}]^{{N_{\rm f}}=2}=0.12(7)$%. For the hyperfine splitting $(m_{V}-m_{P})/m_{P}$ this means a charm quark effect of around 2%. In the $m_{S}/m_{P}$ ratio the central values deviate by $2.7(1.6)\%$.

A clearer difference between the $N_{f}=0$ and $N_{f}=2$ theories can be observed in the mass-scaling functions and in the RGI quark mass. The values of $\eta_{P}$ and the quark mass differ by almost $3\sigma$. The relative differences are $(\eta_{P}^{{N_{\rm f}}=2}-\eta_{P}^{{N_{\rm f}}=0})/\eta_{P}^{{N_{\rm f}}=2}=3.4(1.1)\%$ and $([M_{c}/m_{P}]^{{N_{\rm f}}=2}-[M_{c}/m_{P}]^{{N_{\rm f}}=0})/[M_{c}/m_{P}]^{{N_{\rm f}}=2}=5.0(1.8)\%$. An even larger (but less significant) difference is found in $\eta_{V}$.

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. $m_{V}/m_{P}$ 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 $a\approx 0.049$ 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 $O((a\mu)^{2})$ not only affects observables like $m_{V}/m_{P}$, but also the value of the lattice spacing $a$ itself. Since it is obtained by determining some hadronic length scale $L^{\rm had}/a$ in lattice units at finite lattice spacing and dividing it by the continuum value in fm, i.e. $a=a/L^{\rm had}\times L^{\rm had,cont}$, its value depends on the lattice artifacts present in $L^{\rm had}$. In our case one possibility to compute the lattice spacings is through the scale $L^{\rm had,cont,1}\equiv L_{1}=0.40(1)$ fm. Its values in lattice units are known for our bare couplings and the resulting lattice spacings are between $a^{L_{1}}=0.023$ fm on ensemble $W$ and $a^{L_{1}}=0.066$ fm on ensemble $E$. Alternatively, one could determine the lattice spacing through $L^{\rm had,cont,2}\equiv\sqrt{t_{0}(M)}=0.1131(38)$ 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 $a^{t_{0}}\approx 0.1$ fm, i.e. we observe a 37% lattice artifact in $a$! Since $t_{0}/a^{2}$ is also determined on the quenched ensembles, we can determine their lattice spacings using the decoupling relation $\sqrt{t_{0}(M)}^{{N_{\rm f}}=2}=\sqrt{t_{0}}^{{N_{\rm f}}=0}+O(M^{-2})$. Note that lattice spacings determined by using the $N_{f}=0$ 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 $M$ 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 ${N_{\rm f}}=2$ charm quarks. By comparing to the ${N_{\rm f}}=0$ 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 ${N_{\rm f}}=2+1$ light quarks in the sea. We also demonstrate in figure 8 that lattice spacings $a<0.06~{}$fm are needed for safe continuum extrapolations of the charmonium spectrum when using O($a$) improved Wilson quarks.

We also computed the effects of sea charm quarks in the mass-scaling function $\eta$ of the charmonium masses eq. (3.33) and in the renormalization group invariant charm-quark mass $M_{\rm c}$. Table 4 lists the comparison in the continuum limit of these quantities in the ${N_{\rm f}}=2$ and ${N_{\rm f}}=0$ theory. The effects of the charm sea quarks are clearly resolved and their size is 3% for $\eta_{P}$ and 5% for $M_{\rm c}$. 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 $B_{c}$ 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

[TABLE]

Time reflection

[TABLE]

They are symmetries of the twisted mass action only if $\mu=0$. 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

[TABLE]

With standard Wilson fermions ($\mu=0$), these equations can be used to show that $\langle\mathcal{O}_{1}\mathcal{O}_{2}\rangle=0$ if the operators $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ have opposite parity, i.e. if $\mathcal{P}[\mathcal{O}_{1}]\mathcal{P}[\mathcal{O}_{2}]=-\mathcal{O}_{1}\mathcal{O}_{2}$. 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 $\mathcal{TP}$ transformation is a symmetry of the twisted mass action

[TABLE]

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

[TABLE]

It is now easy to see that the averaged correlator vanishes, if $\mathcal{P}[\mathcal{O}_{1}]\mathcal{P}[\mathcal{O}_{2}]=-\mathcal{O}_{1}\mathcal{O}_{2}$. So, by enforcing the continuum time reflection symmetry of the correlator, automatically the mixing of opposite parity operators is prohibited.

Appendix B Tables