Charm sea effects on charmonium decay constants and heavy meson masses

We estimate the effects on the decay constants of charmonium and on heavy meson masses due to the charm quark in the sea. Our goal is to understand whether for these quantities $N_f=2+1$ lattice QCD simulations provide results that can be compared with experiments or whether $N_f=2+1+1$ QCD including the charm quark in the sea needs to be simulated. We consider two theories, $N_f=0$ QCD and QCD with $N_f=2$ charm quarks in the sea. The charm sea effects (due to two charm quarks) are estimated comparing the results obtained in these two theories, after matching them and taking the continuum limit. The absence of light quarks allows us to simulate the $N_f=2$ theory at lattice spacings down to $0.023$ fm that are crucial for reliable continuum extrapolations. We find that sea charm quark effects are below $1\%$ for the decay constants of charmonium. Our results show that decoupling of charm works well up to energies of about $500$ MeV. We also compute the derivatives of the decay constants and meson masses with respect to the charm mass. For these quantities we again do not see a significant dynamical charm quark effect, albeit with a lower precision. For mesons made of a charm quark and a heavy antiquark, whose mass is twice that of the charm quark, sea effects are only about $0.1\%$ in the ratio of vector to pseudoscalar masses.


Introduction
In recent years there has been a renewed interest in spectral calculations of charmonium states, which are bound systems made of a charm quark and a charm antiquark¯. The motivation is due to the experimental discovery of many unexpected states, whose properties are still poorly understood [1]. Different hypotheses have been put forward, that see them identified as hybrid mesons, tetra-quarks or some other hitherto unknown form of matter. Other systems of particular interest are the mesons, which are heavy mesons composed of a bottom quark (antiquark) and a charm antiquark (quark). Their peculiarities rely on the fact that they are the only heavy mesons consisting of heavy quarks with different flavors. Therefore, they cannot annihilate into gluons and they are more stable, having decay widths less than a hundred of keV. Only the pseudoscalar has so far been seen at the Tevatron collider at Fermilab [2,3]. The high luminosity of the LHC collider at CERN now allows to measure the spectroscopy and the decay of ★ mesons with much better precision and therefore, comparisons of theoretical predictions with experiments will become more and more important.
Lattice QCD may represent a powerful tool to understand the properties of charmonium and mesons directly starting from the QCD Lagrangian, for recent calculations see [4][5][6][7] (charmonium), [8] (charmonium decay constants) and [9] ( mesons). Depending on the problem at hand, numerical simulations usually include the lightest two, three or four flavors. Moreover, exact isospin symmetry is often assumed, implying that the lightest two quarks (up and down) are mass degenerate. The setups mentioned above are usually referred as f = 2, f = 2 + 1 and f = 2 + 1 + 1 QCD. Even though lattice QCD simulations with f = 2 + 1 + 1 dynamical quarks would be highly desirable to have a better understanding of physical phenomena governed by the strong interaction, the addition of a dynamical charm quark to the sea typically complicates the generations of the gauge configurations on which observables are computed. Also, if lattice spacings are not fine enough, the results may be affected by large cutoff effects, owing to the heavy mass of the charm quark. Given the present difficulties described above, in this work we want to quantify how much quenching the charm quark affects the results obtained through f = 2 + 1 QCD simulations (including mass degenerate up, down and a strange quarks). To this aim, we consider a model of QCD, namely QCD with two degenerate heavy quarks, whose mass is tuned to approximately reproduce the physical value of the charm quark mass. In order to quantify the impact of the dynamical heavy quarks, the results from this model are then compared to the ones found using f = 0 QCD (all quarks are quenched). The advantage of our strategy is that reliable continuum extrapolations can be performed, because light quarks are absent and we can use gauge ensembles with extremely fine lattice spacings ( 0.02 fm) and relatively small volumes. In Refs. [10][11][12] we followed this approach to evaluate the impact of a dynamical charm quark on the masses of charmonium states, e.g. on the pseudoscalar and vector meson masses and . From these studies it turned out that such effects are very small and for the hyperfine splitting ( − )/ they are of around 2%.
Here, we want to extend our previous studies and investigate the charm quark sea effects on the hyperfine splitting of heavy mesons (bound states made of a heavy quark and a charm antiquark, or vice versa) and on the decay constants of charmonium. In particular, the focus is on the decay constants and of pseudoscalar and vector mesons, for which we follow the definitions of Ref. [13,14] = 0| 1 2 0 (0)| ( ì = 0) , (1.1) where the bilinears containing quarks 1 and 2 are and | and | represent the ground state of a pseudoscalar and vector meson polarized in the direction of the vector current, 1 respectively. The meson spatial momentum is ì = 0 and |0 is the QCD vacuum state. 2 In the charmonium system, the lowest-lying states lie below the¯threshold, resulting in relatively narrow widths due to the absence of OZI allowed [16] strong decays. This means that radiative transitions, i.e. transitions from an initial state to a final state via the emission of a photon, can have significant experimentally accessible branching ratios. Therefore, the lattice calculation of the decay constants addressed here provides valuable theoretical insight for experiment on the fully nonperturbative level, such that we consider them as natural and representative observables to quantify charm sea quark effects in hadron physics beyond the mass spectrum. The manuscript is organized as follows. In Section 2 we describe how decoupling works for observables which depend explicitely on the charm quark. In Section 3 we present our ensembles of gauge configurations with f = 2 charm quarks in the sea and f = 0 (pure gauge). The observables, from which we extract the decay constants and the heavy meson masses, are explained in Section 4, together with the tuning of the charm quark mass. The latter involves derivatives with respect to the heavy quark mass and we define suitable mass dependence functions in Section 4. Section 5 describes the distance preconditioning for the computation of the heavy quark propagator. The results of our calculations are presented in Section 6 and summarized in the conclusions (Section 7). Appendix A tabulates the values of the decay constants and their derivatives on our ensembles. Appendix B shows a derivation of the formula for the mass derivative in the effective theory which accounts for the dependence of the scale of the effective theory on the mass of the heavy quarks.

Decoupling in the presence of heavy valence quarks
The decoupling of heavy quarks from low-energy physics is well understood [17][18][19][20] and is known to work very well already with heavy quarks as light as the charm quark [21][22][23]. 1  where (. . . , . . .) denotes the standard inner product for spinors, can be approximated by an effective theory, which to leading order is given by QCD without the heavy quarks When the couplings of the effective theory are correctly "matched", they inherit a dependence on the heavy quark mass of the fundamental theory, and the effective theory reproduces results of the fundamental one, as long as the observables do not contain heavy quarks and are low-energy quantities. The differences between the two theories are of order ((Λ/ ) 2 ), (( light / ) 2 ) and (( / ) 2 )), where is the energy scale associated with the observable at hand, e.g. a momentum transfer in a process involving the light quarks, and ,light are the running masses of light and charm quarks. These differences could in principle be suppressed further by adding additional terms to the effective action, e.g.
It is common practice in lattice QCD, to simulate one of the effective theories instead of full QCD. In particular, and quarks are never included in the simulations while the quark sometimes is included and sometimes not. Even without quarks in the action, charm physics is studied in the "partially quenched" approximation. Before we study how well this approximation works numerically, we would like to argue that it is to some extent covered by the theory of decoupling, even though the observables now do depend on the charm quarks. We can add two terms to that amount to a factor of 1 in the path integral While is a Grassmann-valued fermion field, is a complex-valued "pseudo-fermion" field with the same space-time, color and spin components. After integration over the fermion and pseudo-fermion fields, the fermionic determinant is exactly cancelled by the determinant from the pseudo-fermion Gaussian integral. 3  The decoupling mechanism can now be applied to the fields in , leading to a leading-order effective theory After matching of the parameters (gauge coupling, light quark masses and the mass of the "partially quenched" charm quark), one With the mass of the quark being the same as the one that was integrated out, one can of course expect large (1) corrections, but in some cases the energy scale of the observable may still be low, even if it contains heavy quarks. For example, the binding energies in charmonia are not very high, compared to the charm quark mass. In the same respect, another case are their decay constants which turn out to be not very large.

Numerical setup
Our calculations are performed on the gauge ensembles with parameters summarized in Table 1. The f = 0 ensembles are generated using the standard Wilson plaquette gauge action [25], whilst for the f = 2 case a clover-improved [26] doublet of twisted mass Wilson fermions [27,28] is added.
Since we aim at performing continuum extrapolations using very fine lattice spacings, open boundary conditions in the time direction are applied [29] to avoid the well known problems related to the deficient sampling of topological sectors. The spatial dimensions are kept periodic. Moreover, for the production of our ensembles we benefit from the knowledge of the critical mass [30,31] and the axial current and pseudoscalar density renormalization factors [32][33][34] and [30,35]. For further details about the simulations we refer to our previous works [12,22] and references therein.

Computation of meson decay constants with open boundary conditions
Meson decay constants are related to matrix elements between the vacuum and the meson state at rest, as anticipated in Eqs. (1.1) and (1.2). In the framework of the model considered here, we denote the doublet of the heavy degenerate charm quarks in the physical basis as = (˜1,˜2) and its counterpart in the twisted basis as = ( 1 , 2 ) . At maximal twist, the relation between physical and twisted basis is Therefore, the decay constants of pseudoscalar ( ) and vector ( ) mesons in twisted mass QCD are given by  Simulation parameters of our ensembles. The columns show the lattice sizes, the gauge coupling = 6/ 2 0 , the critical hopping parameter, the twisted mass parameter , the pseudoscalar mass in 0 units, the hadronic scale 0 / 2 defined in [24], the lattice spacing in fm from the scale 0 [12] and the total statistics in molecular dynamics units. For the f = 0 ensembles, we report the values of where is the twisted mass parameter. To derive Eq. (4.2), PCVC relations in the twisted basis have been used, as described in Refs. [36,37]. Since both our quarks, 1 and 2 , have the mass of a charm quark, , and , are comparable to the decay constants and masses of the charmonia states and / . The main differences are related to the fact that in our calculations we neglect the contribution of disconnected diagrams and we are assuming a model of QCD without light quarks and with two charm quarks instead of one.
The twisted mass formulation of QCD provides a particularly convenient setup for the calculation of the pseudoscalar decay constant , because it simplifies the renormalization properties of our operators. In fact, the renormalization factors and obey the relation = 1. Thus, the lattice calculation of does not need any renormalization factors, as already discussed in Refs. [27,36]. As concerns the twisted mass expression for , the matrix element in Eq. (4.3) must be multiplied by the renormalization factor of the axial current, which is known from Refs. [32][33][34] for the ensembles considered here.
To extract the decay constants and given in Eqs. (4.2) and (4.3), we first need to know the values of the masses and (already computed in Ref. [12]) and of the matrix elements 0|¯1 5 2 | and 0|¯1 5 2 | . All these quantities can be determined from the lattice calculation of zero-momentum correlation functions of the form where ( 0 , y) and ( 0 , x) denote the coordinates at which a particle state is created (source) and annihilated (sink), respectively. Integrating over the fermion fields, one obtains where . . . is the expectation value over the fermion fields and −1 1 , −1 2 are the fermion propagators of the quark 1 and 2 , respectively. The trace in Eq. (4.6) can be efficiently estimated stochastically. In particular, we use time-dilution with 16 (1) noise sources per time-slice, which amounts to 16 inversions per 0 value and Dirac structure Γ. In principle, the projection to zero-momentum can be realized performing a sum either over ì or ì. However, here we prefer to keep the two sums, as indicated in Eqs. (4.4) and (4.5), because that typically leads to highly improved signals for the meson correlators.
The pseudoscalar and vector meson masses are extracted from the exponential decay of the correlators (4.4) and (4.5), respectively. For this purpose, we first compute the effective mass and then perform a weighted average in the plateau region as where 0 / = , are the initial and final time-slices of the plateau, chosen such that the results do not depend on source position 0 and excited state contributions are sufficiently suppressed. The weights are given by the inverse squared errors of the corresponding effective masses.
In principle, even the matrix elements 0|¯1 5 2 | and 0|¯1 5 2 | can be extracted through an exponential fit to the meson correlators (4.4) and (4.5). However, since we use open boundary conditions in the time direction, such an approach may lead to unreliable results, because of boundary effects. For this reason, we proceed along the method described in Refs. [38,39], whose advantage is to remove the unwanted boundary effects by forming a suitable ratio of two-point correlation functions.
Following the procedure of Ref. [39], when 0 is far enough from the boundaries (0 0 ) and 0 − 0 is sufficiently large, it is possible to show that the correlators and have a leading asymptotic behavior according to 4 √︁ 2 where , are dimensionless factors depending on the source position 0 and related to the matrix elements between the vacuum and the boundary state. From the relations (4.9)-(4.14) one reads off that the relevant amplitudes can be determined through the following ratios of correlators . , where boundary effects and excited state contributions can be neglected. Thus, to improve the estimate of the matrix elements, we take the plateau averages where and are the start and the end of the plateau. In Figure 1 we report the measurement of the effective quantity in Eq. (4.15) obtained on the f = 2 ensemble , with parameters reported in Table 1. As can be seen, the use of lattices with large temporal extent and very fine lattice spacings allows us to take the plateau average¯for a large range of temporal slices. Similar conclusions also hold for the other ensembles and the ratios . 4 The usual relativistic normalization of states in a finite volume is used, i.e. 0|0 = 1, | = 2 3 and | = 2 3 . Finally, combining the formulae above, the pseudoscalar and vector decay constants and in lattice units can be determined as The crucial point is to obtain a reliable estimate of the correlators ( − 0 , 0 ) and ( − 0 , 0 ) in Eqs. (4.15) and (4.16), as when dealing with heavy quarks the relative precision of the solution of the Dirac equation deteriorates at large distances. A detailed discussion of this issue will be given in Section 5. To improve the quality of the signal, we average over forward and backward correlator. Therefore, we consider the symmetrized correlators¯( As it was demonstrated in Appendix A of Ref. [12], taking the average with the timereflected correlators prohibites the mixing with states of opposite parity. 5 The effective masses for heavy mesons made of a heavy quark ℎ ( ℎ = 2 ) and a charm quark on the W ensemble are displayed, together with the plateau average and its error band. We show the pseudoscalar (violet circles) and vector (green squares) channels.

Heavy mesons
To study the charm loop effects at even higher energies, we measure the masses of heavy mesons made of a heavy quark, ℎ, with mass ℎ = { , 2 }. In particular, we focus on the ground state of the pseudoscalar and vector channels. The masses are extracted according to Eqs. (4.7), (4.8). The main ingredient of the calculation is the propagator −1 ℎ , which is evaluated using the twisted mass parameter ℎ = or ℎ = 2 , depending on the value of the mass ℎ we want to impose. Let us recall that is the twisted mass parameter of the simulations tuned to reproduce the charm quark mass.
From now on, we refer to the masses of these heavy mesons using the notation where denotes the vector or pseudoscalar state and the heavy quark, ℎ, can have a mass equal to or 2 . The masses are extracted through the plateau average (4.8) of the effective mass (4.7). For the ℎ = case, we refer to our previous work [12]. An example of this procedure for ℎ = 2 is shown in Figure 2 for the f = 2 ensemble .

Tuning of the twisted mass parameter
To match f = 0 and f = 2 QCD, we use the low-energy observable ℎ = 1/ √ 0 . For such observables decoupling applies [17,40] and we can assume In order to compare the two theories, the value of the renormalization group invariant (RGI) mass of the charm quark needs to be fixed and, following the procedure of our previous work [12], we choose such that √ 0 is approximately equal to the physical value for the meson and is the same in f = 0 and f = 2 QCD, for each lattice spacing employed (see Table 1).
To set we proceed as described below. On our finest f = 2 lattice ( = 6.0), the RGI mass of the charm quark is fixed by the relation where we use the value = 1510 MeV of Ref. [41] (which agrees with Ref. [42], see also Ref. [43]) and the two-flavor Lambda parameter Λ MS = 310(20) MeV known from [30]. Since on our f = 2 ensembles the hopping parameter is set to its critical value, the renormalized physical quark mass is simply given by = −1 and the condition (4.23) is equivalent to fix the twisted mass parameter through In the equation above, the value of the pseudoscalar renormalization factor at the renormalization scale −1 1 in the Schrödinger Functional scheme ( 1 ) = 0.5184(33), valid for 5.2 ≤ ≤ 6, (4.25) and the relation between the running mass ( 1 ) and the RGI mass are available from Refs. [30,44]. As a value of the Λ parameter of the two-flavor theory in 1 units we take [45] Λ MS 1 = 0.649 (45), (4.27) whilst the ratio 1 / is known from [21]. Plugging the values of Eqs. (4.25)-(4.27) into Eq. (4.24), the charm quark mass can be fixed at a precision of around 10%. Such accuracy is acceptable for us if we devise a method to keep the relative mass differences between the different ensembles under better control. In practice, instead of using Eq. (4.24), we tune the twisted mass parameter to a point ★ such that the renormalized quantity √ 0 is equal to its value on our finest f = 2 ensemble W at = 6.0 (see Table 1) and satisfies √ 0 This value is obtained on the ensemble W by setting the twisted mass parameter through Eq. (4.24). The tuning point Eq. (4.28) is independent of the overall scale Λ MS (known with a 7% accuracy and which limits the overall accuracy).
In the f = 2 theory, we apply a correction to all observables, which is based on the computation of twisted mass derivatives [12]. First, we determine the target tuning point ★ through the Taylor expansion where and √ 0 are the values from the simulations obtained using Eq. (4.24), see Table 1. Subsequently all quantities, denoted by below, are corrected by For a generic primary observable , its twisted mass derivative is simply given by where is the lattice QCD action. However, most quantities we are interested in are non-linear functions of various primary observables, for which we apply the chain rule (4.32) After integrating over the fermion fields, Eq. (4.31) becomes where · gauge are expectation values taken on the ensemble of gauge fields and eff is the effective lattice QCD gauge action which includes the effect of the fermion determinants. The observable˜contains in general the fermion propagator −1 . The sea quark effects originate from eff . They affect the mass derivative of an observable through the covariance of eff and˜.
In the simplest case of f = 0 QCD, the action does not depend on the quark masses and the twisted mass parameter only enters the inversions of the Dirac operator. The mass derivative of an observable in the f = 0 theory therefore only receives a contribution from the last term on the right hand side of Eq. (4.33). Thus, to reproduce the tuning point ★ , Eq. (4.28), for our f = 0 ensembles, we measure the quantities of interest at three different values of the twisted mass parameter such that ★ is eventually obtained through a linear interpolation of the measurements.

Mass dependence functions
The mass dependence of a multiplicatively renormalizable quantity is captured by the renormalized combination ≡ . (4.34) Interesting quantities are for instance the vector or pseudoscalar meson masses as studied in [12], or their decay constants, as determined in this work. In the twisted mass formulation at maximal twist, the RGI quark mass is connected multiplicatively to the twisted mass parameter , which implies that can be replaced by in Eq. (4.34). When computed in the effective theory, the mass dependence stems from an explicit mass dependence of the quantity in question, but also from the mass dependence that the scale inherits through the matching condition between the effective and full theory. In this context it is useful to also define˜≡ Λ=const. . (4.35) In the fundamental theory, Λ is constant and˜= , but in the effective theoryc aptures only the explicit mass dependence of the quantity in question and vanishes e.g. for purely gluonic quantities. It has been computed in [12] for the masses 6 and found to deviate significantly from the corresponding two-flavor values. Most of this deviation is explained by the neglected mass dependence associated with the scale. If the mass parameter in the f = 0 theory is fixed by demanding that has the same value as in the f = 2 theory, the full mass dependence is given by In this equation, denotes the universal mass scaling function introduced in [46], which encodes the complete mass dependence of purely gluonic quantities. Its definition is where 2,0 is defined in Eq. (B.3) and ≡ (2) . A derivation of Eq. (4.36) is relegated to Appendix B.

Distance preconditioning for the Dirac operator
The computation of and in Eqs. and ( ) denote the solution vector and the source field (the latter being non-zero only on a single time-slice 0 ), respectively. Whereas the increase of computational demand for the numerical inversion of the Dirac operator to solve this linear system towards small quark masses is a consequence of the spectral properties of ( , ), a difficulty of entirely different nature arises in the region of heavy quark masses.
To understand the origin of this potentially serious problem, we recall that the exponential decay with time of two-point functions at zero spatial momentum is governed by the ground state energy in its spectral decomposition, which equivalently reflects in the exponential time dependence of the entering quark propagators proportional to exp{− ( 0 − 0 )}, where is the quark mass in the corresponding Dirac operator. However, the exponential decay of the quark propagator is particularly pronounced in the case of heavy quarks, = ℎ , when the time distance between source and sink grows. This entails that an accurate determination of heavy meson correlators over the relevant time span is very difficult to achieve, because, depending on the size of the heavy quark mass in lattice units, the relative precision of the solution of the Dirac equation could actually start to deteriorate at already moderately large time separations 0 − 0 . Moreover, this problem is amplified for correlation functions in the vector channel, in which the ground state is even heavier.
In practice, for heavy quarks, the stopping criterion commonly imposed on the relative global residuum may therefore no longer appear as a reliable quantitative measure to indicate the convergence of the iterative routine used to numerically solve the Dirac equation to desired accuracy. In fact, for large time separations 0 − 0 , contributions to the norm on the l.h.s. suffer from a severe exponential suppression ∝ exp{− ℎ ( 0 − 0 )} and thus become negligible. Even a "brute force" approach of reducing the global residual gl further to catch this may then not be feasible anymore for lattices with very large time extents, regardless of the related increase in computational cost.
To overcome this unwanted exponential suppression problem induced by the heavy quark propagator, we employ the distance preconditioning (DP) technique for the Dirac operator, originally proposed in Ref. [47], in its implementation outlined in [48]. The basic idea is to rewrite the original linear system in such a way that the solution of the new system at time-slices 0 far away from the source gets exponentially enhanced and thereby compensates the rapid decay of the (heavy) quark propagator. While in [47] the associated re-definition of the propagator and source field was introduced on the level of the covariant derivatives in the lattice Dirac operator, the DP implementation of Ref. [48] directly works with the matrix-vector equation to be solved, since this has proven to be most straightforward for, e.g., the SAP-GCR solver routine [49] in use. Schematically written, DP then amounts to replace where in our setup the preconditioning matrix is defined as Here, 0 and ( ) 0 refer to the time-slices of the source and sink insertion, and apart from this, is unity in spatial coordinates as well as in spin and color spaces. The solution of the original system is then simply obtained by = −1 from the solution of the preconditioned one. In Eq. (5.3), 0 is a tunable parameter that needs to be adjusted for each gauge field configuration ensemble in dependence of the respective temporal lattice extent and heavy quark mass in lattice units.
Rather than the global residual, Eq. (5.1), a more suitable measure for the numerical quality of the solution over the whole time range is now provided by inspecting the local residuum, at the time separation of interest ( − 2 0 in our case). Since loc is sensitive to the numerical accuracy of the solution on each time-slice, imposing the stopping criterion of the routine involved to solve the preconditioned system ensures to extract heavy meson correlation functions to sufficient precision also for large time separations. Note that the exponential factor 0 ( 0 − 0 ) in the construction of the preconditioned linear system, Eq. (5.2) above, effectively turns the heavy quark propagator into a propagator corresponding to a "ficticious" light quark. Therefore, the price to pay when adopting this method is an increase of the number of solver iterations such that in practical applications the growth in computational costs has to be reasonably balanced with the gain in accuracy of the heavy meson correlator in question by a proper tuning of 0 . As can be seen from monitoring loc for a representative ensemble and various choices of 0 in Figure 3, the solution becomes more and more accurate as 0 increases. However, in this example we observe that for the solution to have a local residual < 10 −4 for all time-slices (which is enough given the desired statistical accuracy of the meson correlator) the number of solver iterations increases of around a factor 3 compared to the case 0 = 0 (no distance preconditioning applied). Such a small value of the residual is crucial to  Table 1) . We show the relative local residual of the solution as a function of the sink position 0 with respect to the source set at 0 / = 16. On the right, the number of iterations for the solver to converge is reported for the values of 0 plotted in the Figure. extract the meson correlator (and as a consequence the meson decay constants) reliably. Ideally one would choose 0 = , but to keep the computational effort as small as possible, we explored different source positions (with 0 > ) to ensure the ground state dominance over a large range of time-slices and a reasonable number of iterations for the solver to converge. Both and have been determined for all the ensembles under study and the parameters used in the solver setup are collected in Table 2.

Decay constants of charmonium
Combining the measurements of the ratios and , defined in Eqs. (4.15) and (4.16), with the charmonium masses that we previously obtained in Ref. [12], through Eqs. (4.18) and (4.19), we determine the dimensionless quantities Similarly to what we observed for the meson masses in Ref. [12], we find that linear fits in 2 do not work well in the range of lattice spacings 0.02 fm 0.07 fm, if one aims at results with sub-percent precision. Therefore, we compare the continuum limits of pseudoscalar and vector decay constants in f = 0 and f = 2 QCD, performing extrapolations to zero lattice spacing only in the range 0.02 fm 0.05 fm. The results for the f = 0 (empty markers) and f = 2 (full markers) ensembles are shown in Figure 4. The continuum limit values are slightly displaced for clarity.
Comparing the values obtained in the continuum of the two theories, one can observe that dynamical charm effects on the decay constants are barely resolvable, despite the great accuracy of our continuum extrapolations. In Table 3 we summarize our findings. The relative difference of [MeV]   Note that the decay constant of the meson has not been determined experimentally, whilst for the vector meson / the value / = 407(4) MeV is obtained from the partial decay width of / into an electron-positron pair, see [8]. / can be compared to our lattice results for in Eqs. (6.1) and (6.2). The discrepancy is probably due to effects of light sea quarks, disconnected contributions and electromagnetism, which are neglected in this work, and also to the unphysical values of the meson masses. When the scale dependence is taken into account, the errors increase due to a dependence on (2) .
The mass dependence functions for meson masses and decay constants have been computed as well. Figure 5 shows how (2) compares to both (0) and˜( 0) for various quantities. In (2) and˜( 0) , the derivatives are computed using Eqs. (4.32) and (4.31). The missing pieces to evaluate Eq. (4.36) are the continuum value (2) = 0.6996(81) from [12], and the value ( ) = 0.1276(2) determined as in [21] for the f = 2 → f = 0 case,  The mass dependence functions for various quantities in the continuum limit, at the charm quark mass. The results for the meson masses rely on data published in [12]. Note that in the first row (2) = (0) because of the matching condition (B.2).
using four-and five-loop perturbation theory results for the decoupling relations [19,20], the QCD -function and the QCD anomalous dimension [50][51][52][53][54]. Table 4 summarizes the findings on the mass dependence. Note that agreement is only observed when the mass dependence is consistently taken into account. The charm sea quarks effects in the mass derivatives of an observable originate from the covariance of the mass derivative of the effective action with the observables, see Eq. (4.33). This covariance is absent in the f = 0 theory and explains the much larger statistical errors of the f = 2 theory. The derivatives computed on each ensemble are used to shift the decay constants to the tuning point, see Eq. (4.30). Still, they are physical quantities on their own.

Heavy meson masses
In this section we present our results for heavy mesons made up of a heavy quark ℎ ( ℎ = 2 ) and a charm quark , as described in Section 4.2. To guarantee that the cutoff effects in the continuum extrapolations are under control, we only consider the ensembles of Table 1 for which the meson masses satisfy the relation < 1. In Figure 6 we compare the continuum limits of the mass ratio , / , in f = 0 and f = 2 QCD. The numerical results of the continuum extrapolations are summarized in the first row of  Results for , / , in the continuum limit for both the f = 0 and the f = 2 theory. The values for the ℎ = case are taken from Ref. [12]. ratio , / , is found to be 0.092(50)% (around 1.8 ), which is of similar size to the one obtained in our previous work [12] for ℎ = . In Fig. 7  to the f = 2 theory and the red squares to the pure gauge theory, f = 0. We take the Λ parameters in units of the scale √ 0 [24] from [21,30] for f = 2 QCD and from [55] for f = 0 QCD. The RGI quark mass values ℎ = at our tuning point Eq. (4.28) have been determined in [12] for both theories. The data of Figure 7 are listed in Table 6.
The limit ℎ → ∞ is expected to be described by Heavy Quark Effective Theory (HQET). For f = 0 QCD we show in Figure 7 a linear function (dashed line) going through the points at ℎ = ∞, where the vector and the pseudoscalar are degenerate by virtue of the heavy quarks spin symmetry [56][57][58], and ℎ = 2 . Contact with HQET in the spirit of [59], including the HQET-QCD matching factor spin that accounts for the correct mass dependence of the ℎ → ∞ asymptotics, would require larger heavy quark masses than we have simulated in this work.

Conclusions
In this work we present a study of the effects of sea charm quarks in the charmonium decay constants. It is a follow-up of Ref. [12] where we did a similar study for the charmonium spectrum and the RGI charm quark mass. We compute the charm-quark sea effects through which approximately corresponds to a physical charm quark. We find that the effects of two charm sea quarks are below 1% for the decay constant of the pseudoscalar ( ) and of the vector ( / ) ground state particles. We also extracted the derivatives of the decay constants and meson masses with respect to the charm quark mass, because they are required to shift our simulation results to the tuning point Eq. (4.28). The quantities = ( = , , , ) are renormalized and easily computable in the twisted mass formulation at maximal twist, because the twisted mass parameter is multiplicatively renormalizable. In the continuum limit we do not see significant effects of the dynamical charm quark on these quantities.
A computation of the charm sea effects for " " mesons made of a charm quark and an antiquark of twice larger mass reveals that they are only about 1‰ for the ratio of the vector to the pseudoscalar mass.
Our results are relevant for the scale setting [38] of the f = 2 + 1 simulations by the CLS consortium [60,61]. There, a combination of the pseudoscalar decay constants of the pion and kaon is used. They are lower in energy than the decay constants of charmonium. Hence our results indicate that charm sea effects are below the precision of the scale determination in [38]. Finally, the analysis of the " " mesons also demonstrates that charm sea effects are very small for splittings of heavy mesons. For instance, the Υ − Υ splitting has been used in [62] for scale setting.

A Table of decay constants
All the results obtained on our ensembles are listed in Table 7.

B Mass dependence of fermionic observables
Eq. (4.36) describes how quantities with an explicit valence quark mass dependence depend on the quark mass, if decoupling is in place. The parameters of the leading order effective theory, here f = 0 QCD, are fixed by demanding non-perturbatively √︃ .
(0) Λ=const . The latter captures the valence quark mass dependence, but ignores the mass dependence of the scale in the effective theory.