Scale setting for $N_f=3+1$ QCD

We present the scale setting for a new set of gauge configurations generated with $N_f=3+1$ Wilson quarks with a non-perturbatively determined clover coefficient in a massive O($a$) improvement scheme. The three light quarks are degenerate, with the sum of their masses being equal to its value in nature and the charm quark has its physical mass. We use open boundary conditions in time direction to avoid the problem of topological freezing at small lattice spacings and twisted-mass reweighting for improved stability of the simulations. The decoupling of charm at low energy allows us to set the scale by measuring the value of the low-energy quantity $t_0^\star/a^2$, which is the flow scale $t_0$ at our mass point, and comparing it to an $N_f=2+1$ result in physical units. We present the details of the algorithmic setup and tuning procedure and give the bare parameters of ensembles with two lattice spacings a=0.054 fm and a=0.043 fm. We discuss finite volume effects and lattice artifacts and present physical results for the charmonium spectrum. In particular the hyperfine splitting between the $\eta_c$ and $J/\psi$ mesons agrees very well with its physical value.


Introduction
In [1] a lattice action was proposed to simulate Quantum Chromodynamics (QCD) with N f = 3 + 1 quarks on the lattice. The renormalization and improvement conditions are imposed at finite quark masses. The masses of the up, down and strange quarks are degenerate and their sum is approximately equal to its value in nature. This not only reduces the number of mass parameters, but also facilitates the simulations not having very light up and down quarks. The mass of the charm quark is close to its physical value. Quantities like the charmonium spectrum are fairly insensitive to our approximation of degenerate light quark masses. The action in [1] was designed to simplify the O(a) (a is the lattice spacing) improvement pattern for Wilson quarks [2] in the presence of a dynamical charm quark. The improvement coefficients are determined following the Symanzik improvement program. In a massive scheme for Wilson fermions the Sheikholeslami-Wohlert coefficient c sw (g 2 0 , aM ) [3] in the action depends, besides on the gauge coupling g 2 0 = 6/β also on the quark masses (aM being the bare quark mass matrix). The same is true for the improvement coefficients of operators. The main result of [1] was a non-perturbative determination of c sw using a finite volume scheme. In this work we use the action for large volume simulations to test its scaling properties and provide first physical results for the charmonium spectrum. The hyperfine splitting m J/ψ − m ηc between the ground-state vector and pseudoscalar mesons made from two charm quarks is notoriously sensitive to lattice artifacts stemming from the discretization of the Dirac operator [4][5][6], and therefore a good measure of the quality of the action.
One result of our work is the scale setting for N f = 3 + 1 QCD simulated with the action of [1]. Scale setting provides a relation between the bare coupling of the simulations and the lattice spacing in fm. The scale t 0 (mass dimension -2) is the flow scale t 0 [7]  In our setup the pion and kaon as well as Dand D s -mesons are degenerate. In N f = 3 eq. (1.2) plays no role, in our N f = 3 + 1 simulations discussed here, it is important, that the value of φ 5 corresponds to a charm quark mass, about the same as the one found in nature. We measure t * 0 /a 2 at a mass point defined by eq. (1.1) and eq. (1.2) and assume that 8t 0 = 0.413 (5) (2)fm to obtain the lattice spacing in fm. This particular value was determined in [8,9] in N f = 2 + 1 QCD by determining the ratio of t −1/2 0 , at the unphysical mass point, with a linear combination of pion and kaon decay constants at the physical mass point, in the continuum limit. Relying on a three flavor result is justified by the fact, that when dealing with low-energy quantities like t 0 , our N f = 3 + 1 theory can be well described by an effective theory which to leading order is QCD with just N f = 3 light flavors. A detailed study of non-perturbative decoupling of the charm quark was performed in a model, QCD with N f = 2 mass-degenerate heavy "charm" and zero light quarks [10][11][12][13]. The effects of the decoupled charm quark at low energies are incorporated in the matching of the gauge coupling and the (less relevant) quark masses. Power corrections to decoupling are of the order (E/m charm ) 2 and (Λ/m charm ) 2 and stem from higher dimensional terms in the effective Lagrangian. An important result of the decoupling study was that dimensionless low energy quantities are insensitive to the presence of a dynamical charm quark at our level of accuracy. Indeed, only effects at the permille level were found. We discuss the simulation setup and tuning procedure in section 2 and introduce the observables we measure in section 3. In section 4 we explain the scale setting of our ensembles with two lattice spacings and discuss systematic errors in section 5. These ensembles are not only important for a precise determination of the strong coupling constant, but will be very useful for many other physics applications, e.g. for high precision charm physics. As a first result we therefore present charmonium spectra in section 6 and extract charmonium masses and hyperfine splitting in good agreement with PDG [14] values. We end with concluding remarks and a short outlook in section 7.

Simulations
The ensembles were generated using the open-source (GPL v2) package openQCD version 1.6 [15] in plain C with MPI parallelization. This software has been successfully used in various largescale projects. The measurements of the hadronic correlation functions were carried out with the open-source (GPL v2) program "mesons" [16], which is based on various openQCD modules, in particular the openQCD inverters, and hence shares its architecture (C+MPI) and good scalability.
The simulations are performed on lattices of size N t ×N 3 s . Gauge fields are periodic in spatial directions and absent on temporal links sticking out of the lattice, i.e. we have open boundaries in temporal direction imposed on time slice 0 and N t − 1, hence the physical time extent is T = (N t − 1)a with the lattice spacing a. For a general introduction to lattice QCD simulations with open boundary conditions and twisted-mass reweighting we refer to [17,18].

Action
The full action reads S = S g + S f , with a gauge action given by the Lüscher-Weisz action [19,20] with tree-level coefficients c 0 = 5/3 and c 1 = −1/12, where the summation is over all oriented plaquettes p and rectangles r (double-plaquettes) on the lattice. U p/r is the product of four/six SU(3) gauge fields U µ (x) around the elementary plaquette p or rectangle r. The free parameter of the gauge action is the bare coupling g 2 0 ≡ 6/β and the weights w(p) = w(r) = 1 except for spatial plaquettes and rectangles at T = 0 and T = (N t −1)a where w(p) = w(r) = 1/2, i.e. the coefficient of the gluonic boundary improvement term is set to its tree level value c G = 1. We do the same with the coefficient of the fermionic boundary improvement term c F = 1. The fermion action then reads where ∇ µ and ∇ * µ are the covariant forward and backward derivatives, respectively, and the Sheikholeslami-Wohlert action S SW in eq. (2.2) is needed for O(a) improvement [3].
Simulating a physical charm quark results in lattice artifacts due to large cutoff effects of order proportional to a large value of the lattice charm quark mass (am c ≈ 0.5). In a mass independent scheme the mass dependent O(a) lattice artifacts are cancelled by improvement terms where a mass independent coefficient multiplies aM [22,23]. Some of these coefficients are only known to one loop in perturbation theory which may be acceptable for the light quarks but not if the quark mass is large. Therefore a massive renormalization scheme with close to realistic charm mass m c and Symanzik improvement for Wilson quarks was proposed [1]. This results in massdependent renormalization parameters for coupling and quark masses as well as a mass-dependent coefficient in the clover action (F µν is the standard discretization of the field strength tensor [22]), which we determine using the non-perturbative fit formula obtained in [1] c sw (g 2 0 , aM ) = While in a full massive scheme the renormalization and improvement factors would depend on all quark masses, in practice a hybrid approach is more feasible, where the dependence on the heavy quarks (just the charm quark in our case) is absorbed into the factors, but light quarks are treated like in a mass-independent scheme. This allows to use the same value of c sw along chiral trajectories with constant charm quark mass, while at the same time avoiding large O(am c ) artifacts. In fact, eq. (2.4) can be used whenever up, down and strange quarks are light, and the charm quark has close to its physical mass. In this work we describe simulations at the SU (3) flavor symmetric point, where M = diag(m l , m l , m l , m c ) 1 .

Algorithms
The generation of these configurations proceeds according to a variant of the Hybrid Monte-Carlo (HMC) algorithm [24]. The classical equations of motion are solved numerically for trajectories of length τ = 2 in all simulations, leading to Metropolis proposals which are accepted with an acceptance rate P acc . In order to reduce the computational cost and obtain a high acceptance rate, we split the action and corresponding forces as explained in detail in this section.
Since the Wilson Dirac operator is not protected against eigenvalues below the quark mass, we use the second version of twisted mass reweighting, suggested in ref. [25], for the asymmetric even-odd preconditioned [28] Dirac operator The u/d quark doublet 2 is then simulated with a weight proportional to where we introduced an infrared cutoff by the twisted mass µ 0 . To further reduce the fluctuations in the forces, we factorize the last term in eq. (2.6) according to [29][30][31] , with aµ i given by {0.0005, 0.005, 0.05, 0.5} for all our ensembles. The individual factors can be represented by pseudo-fermions such that the combination of twisted-mass reweighting and factorization leads to a pseudo-fermion action for the light quark doublet with N + 2 = 6 terms where φ i are six pseudo-fermion fields with support on the even sites of the lattice. The u/d quark doublet comes with a reweighting factor which is estimated stochastically. The strange and charm quarks are simulated with the RHMC algorithm [26,27], decomposing into reweighting factors (to be estimated stochastically again) and Zolotarev optimal rational approximations of (D † qDq ) −1/2 [32] in the spectral ranges [r a , r b ] q ofD † qDq with a number of poles N q p , optimized during the tuning process, which uniquely determine the parameters A q , ω q and ν q . The Zolotarev rational function eq. (2.11), in the case of the strange quark, is further broken into several factors, introducing separate pseudo-fermion fields for the five smallest ω s,k , ν s,k , whereas the determinant of the remaining factors (first seven poles) is expressed as a single pseudo-fermion integral. In practice, this product is rewritten via partial fraction decomposition into a sum The effective action for the strange quark with N s p = 12 poles reads The charm quark contribution to the action is not further factorized such that we have 13 pseudo-fermion fields and 14 actions in total, with forces summarized in table 9. The molecular dynamics equations are integrated using multi-level higher order (2nd and 4th order OMF) integrators [33]. Linear systems involving the Dirac operators are dealt with using low-mode-deflation [34] and SAP-preconditioned [35] Krylov space solvers based on local coherence [36,37], the tuning of corresponding parameters is discussed in the next section.

Parameter tuning
For the algorithmic parameters, we started with the setup of CLS's H400 simulation, cf. [38], to which we added the charm quark. The new contribution to the action was not further factorized and the corresponding forces were integrated on the intermediate level of our three level integrator. The gauge fields are integrated on the innermost level with the fourth order integrator suggested by Omelyan, Mryglod, and Folk (OMF) [33], most of the fermion forces are on the intermediate level (4th order OMF) and only the most expensive doublet factor (the one with φ 0 in eq. (2.15)) and four largest poles of the strange rational function are on the outermost (coarsest) level using a second order OMF. For most fermion forces we use the locally deflated solver [34,36,37] with a deflation subspace block size 4 4 and 24 deflation modes per block. We set the relative residua to 10 −12 in the acceptance-rejection step and to 10 −11 in the force computation. The parameters of the rational functions, i.e. number of poles and spectral ranges, are adjusted during thermalization by calculating the reweighting factors W s,c and the true spectral range of |γ 5D | for a subset of the generated gauge field configurations. The simulation and algorithmic parameters are summarized in table 9.

Wilson-flow observables
The Wilson flow can be used to define quantities with a finite continuum limit in lattice QCD, starting out with a smoothing flow equation [7,39] with smeared gauge field V µ (x, t) at flow time t and the corresponding Wilson gauge action Using a symmetrized clover definition of the field strength tensorĜ µν (x, t) constructed from the smooth fields V µ (x, t), the time slice action density E(x 0 , t) and the global topological charge Q(t) can be defined as Since significant boundary effects in E(x 0 , t) and Q(x 0 , t) were observed in [38], we take a plateau average in a range of x 0 values away from the boundaries (given in table 3) to define the vacuum expectation value of the energy E(t) ant topological charge Q(t). Correlators constructed from gauge fields at t > 0 are automatically renormalized [40], which allows to define low-energy scales t 0 [7], t c and w 2 0 [41] as the flow times t at which respectively, the partial derivative in the case of w 0 is evaluated numerically.

Meson correlators and masses
We want to extract meson masses from the zero momentum correlation functions Integrating out the fermions leaves us with a single, connected diagram of the form withΓ ≡ γ 0 Γ † γ 0 . The trace acts in color and spin space and the propagators S i are given by the inverse Dirac operator An efficient way to carry out the calculation is to use stochastic sources, we use 32 noise vectors per time-slice, which amounts to 32 inversions per y 0 value and Dirac structure. In Table 1 we list the meson interpolators for the various states investigated and their particle names which are the closest relatives in nature. An improved signal and exact symmetries are achieved by defining the averages The meson masses are extracted from the exponential decay of these correlators at x 0 a via weighted plateau averages of the effective masses . (3.14) and weights w given by their inverse squared errors. Table 1 also lists the plateau ranges for the different ensembles, chosen such that excited state contributions are negligible.

Reweighting, mass derivatives and tuning corrections
With reweighting, QCD expectation values are obtained from the simulated ones via with corresponding weights W = W ud W s W c given by the product of reweighting factors eq. (2.8) and eq. (2.10) estimated stochastically. We are interested in functions of primary observables and their quark mass derivatives df /dm. Therefore we first reweight our primary observables according to eq. (3.15) and calculate their reweighted mass derivatives via [9] d  where ∂S ∂mq = xq (x)q(x) gives rise to x tr [S(x, x)] after integration over the quark fields in the third term of eq. (3.17) and also in the second term, if no other quark fields are present in the observable O i . This trace is evaluated stochastically using 32 U(1) random sources per configuration. When O i contains fermionic fields, the second term in eq. (3.17) produces additional Wick contractions that depend on O i . For the case of the meson correlators these are worked out and measured. For derived observables f , we then use the chain rule The partial derivatives can be calculated analytically or estimated numerically where a suitable h is obtained from the fluctuations of O i , and we do not have explicit derivatives ∂f ∂m for the derived observables considered here. In order to correct for mis-tunings of φ 4 eq. (1.1) and φ 5 eq. (1.2), we solve the following system of linear equations to get the mass shifts ∆m l and ∆m c for light and charm quarks and use these to correct arbitrary observables f accordingly (3.20)

Scale setting
For starting parameters we choose a bare coupling β = 3.24, light quark masses given by κ u,d,s = 0.134484 and a charm quark mass by κ c = 0.12, obtained from a rough interpolation of the simulation parameters in [1]. We thermalize on spatially smaller lattices and subsequently double the spatial dimensions, using eight copies of the smaller lattices as the new initial configurations. Unfortunately, the smaller lattices are not helpful for parameter tuning, since we are dealing with large finite volume effects, as discussed in section 5.2. We measured the flow observables and meson masses on a more-or-less thermalized subset of configurations of the final lattice size to evaluate φ 4 and φ 5 . With the starting parameters the tuning point was missed by quite a bit, therefore, we simulate several other mass points and calculate the derivatives of φ 4 and φ 5 with respect to bare quark masses, in order to gradually reach the correct tuning point. This is non-trivial, since the mass dependence of t 0 and the meson masses in φ 4 and φ 5 go in opposite directions, see table 10 and 11 for results on the ensembles A2 and B. With the final parameters κ u,d,s = 0.13440733 and κ c = 0.12784, at β = 3.24 we produced two high statistics ensembles A1 and A2 with two different lattice sizes, and a smaller ensemble A0 on an even smaller volume to study finite size effects. Assuming decoupling of the charm quark, i.e., , our value of t 0 /a 2 ≈ 7.4 corresponds to a lattice spacing a ≈ 0.054 fm, where we used the N f = 3 result 8t 0 = 0.413 (5)(2)fm from [8,9]. This makes our smaller volume of ensemble A1 (L/a = 32) L ≈ 1.73 fm with m π L = 3.5, which is a bit small, but finite size effects seem to be under control, as the comparison with L/a = 48 shows, see further section 5.2. In a next step we tuned an ensemble B at a finer lattice spacing on a 144×48 3 lattice with a bare coupling β = 3.43 using the same procedure as described above, yielding final parameters κ u,d,s = 0.13599 and κ c = 0.13088, and a lattice spacing a ≈ 0.043. This gives a physical lattice extent L ≈ 2.06fm and m π L = 4.3.   Simulation parameters are summarized in Table 2, flow measurement and topological charge results are presented in Table 3 and corresponding masses and tuning results are shown in Table 4. The second row values for ensembles A2 and B are shifted to the correct tuning points φ 4 and φ 5 , using the strategy described in section 3.3. The corresponding quark mass shifts a∆ m l and a∆ mc are given in Table 5 together with the final lattice spacings.

Lattice artifacts
A major challenge in the simulation of QCD with heavy quarks, is the control of lattice artifacts. In order to asses their size for our ensembles, and to find out whether they behave like the expected ens. N ms    leading scaling violations, we look at two slightly different definitions of t 0 in eq. (3.4), using either Wilson's plaquette discretization of the action density E(t) in eq. (3.2), or a more symmetric definition based on the symmetric "clover" discretization of the field strength tensor, see [7]. Both definitions of t 0 have to agree in the continuum limit, but may differ at finite lattice spacing, hence, the ratio of t clov 0 and t plaq 0 has to be one up to cutoff effects. The two values of t 0 are highly correlated and their ratio can be evaluated very accurately. Figure 1 shows one minus the ratios for ensembles A2 at β = 3.24 and B at β = 3.43 (shifted to the correct tuning point using eq. (3.20)), supposed to vanish in the continuum limit at a rate proportional to a 2 . Due to the extremely high precision of this particular observable, the data cannot be described by a pure a 2 function shown in the left plot of Figure 1, as that would lead to an unacceptably large χ 2 /dof = 19.2. However, in absolute terms, the deviation from pure a 2 scaling is tiny, only about 1.7 permille on the coarser lattice for this particular observable, see the right plot of Figure 1. Given that gradient flow observables in general and the plaquette definition of the action density in particular are known to have large lattice artifacts (cf. [7]) , we are convinced that for the majority of observables continuum extrapolations linear in a 2 will suffice.

Finite volume effects
Next, we study finite volume effects of am π following [42,43], who propose an analytic scaling formula from chiral perturbation theory (χPT) in the "p-expansion" (5.1) We take the pion mass and decay constant in ξ π at the SU(3) flavor symmetrical point determined on the finest lattice in [9], m π (∞) our result from ensemble A2, N f = 3 and m(n) are integer multiplicities listed in [42] for n = 1 . . . 20. In the left plot of Fig. 2 we show the analytic χPT prediction together with our results on ensembles A0, A1 and A2 at the same lattice spacing. In the range of pion masses and volumes considered the agreement between the one-loop analytical prediction and our lattice data is poor, especially for Lm π < 3.5. Finite volume effects seem to be under control for m π L > 4, as can be seen by direct comparison of ensembles A1 and A2, see e.g. table 3-8, most values agree within errors. Another source of finite volume effects stems from the finite temporal extent with the open boundaries. Translational invariance is broken and only a portion of the temporal lattice sufficiently far away from the boundaries can be used in the averages. The right plot of Fig. 2 shows these boundary effects for some of the flow observables and the topological charge squared.

Decoupling of charm
How well decoupling of a dynamical charm quark holds has been studied in a sequence of papers [10][11][12][13] Based on the outcome of these studies we are optimistic that using the t * 0 value from the N f = 3 simulations introduces negligible errors for N f = 3 + 1 QCD with a physical charm quark. The studies were carried out in a model, namely N f = 2 QCD with two degenerate charm quarks and no light quarks. Here we are in the position to investigate decoupling in a more realistic setting with three degenerate light quarks. We look at dimensionless quantities from ratios of flow observables t 0 /t c and t 0 /w 2 0 and attempt continuum limit extrapolations from the corresponding values from ensembles A2 and B, listed in table 8. In figure 3 we plot the ratios and their extrapolations and compare with corresponding extrapolations on N f = 3 CLS ensembles [38]. The ratios were formed from flow measurements with clover (black) and plaquette (red) definitions of the action density E(t). The continuum extrapolations are performed by constrained linear fits which take the correlations into account. The deviations between N f = 3 + 1 and N f = 3 in the continuum are far below 1% and of the order of magnitude found in the model study [11]. We notice that for a reliable continuum extrapolation of the N f = 3 data the combination of the two different definitions of the action density and fine lattice spacings 0.05fm are necessary. While the complicated a-dependence of observables relying on the 0 (green circles) and squared topological charge plaquette definition of the action density might seem intimidating, we would like to point out that this is probably a "worst case observable" which is avoided in practice, for good reasons. O(a 2 ) improved variants [44] of both the action density and the flow equations should be employed, when possible.

Autocorrelations
When approaching small lattice spacings, the HMC algorithm is known to suffer from critical slowing down. Autocorrelations grow at least ∝ a −2 , but have been observed to grow much faster for quantities like the topological charge [17]. One method to circumvent the freezing of the topological charge, is to use open boundary conditions in one of the lattice directions [17], and that is what we are doing here. As can be seen in Fig. 4 the topological charge moves freely, as do other "slow" quantities like the flowed action density. Other quantities exhibit the expected slowing down ∝ a −2 , which however is still manageable on our finest lattice. E.g. we find integrated autocorrelation times τ int,Q 2 ≈ 8 ± 3 [4 MDU] for the topological charge squared and τ int,t 0 ≈ 35 ± 21 [4 MDU] for t 0 . These values are large, but small enough compared to our total statistics. This is corroborated by fits to the tail of the autocorrelation functions which allow us to estimate the exponential autocorrelation time τ exp (corresponding to the slowest mode in the Markov chain) listed in Table 2. For our error analysis we use the Γ-method [46] adding a tail to the autocorrelation function to account for τ exp [47]. For most observables the errors are only marginally different with respect to the standard analysis of [46].

Chiral extrapolation effects in charmonia
Although our N f = 3 + 1 ensembles are not at the physical mass point, this work also provides a formidable starting point for N f = 2 + 1 + 1 simulations of QCD with improved Wilson quarks.
The quantities φ 4 and φ 5 are chosen such, that the (renormalized) charm quark mass m c and the sum of the degenerate (renormalized) light quark masses m u + m d + m s are very close to their physical values. The physical mass point can be approached along chiral trajectories where m u = m d is decreased and m s is increased, while the trace of the quark mass matrix is kept constant, i.e. the sum of the differences to the physical masses is zero (∆ mu + ∆ m d + ∆ ms = 0). The chiral extrapolations are expected to be very flat for quantities that do not contain light (u,d,s) valence quarks. The derivatives of these quantities with respect to the light quark masses are equal to each other at the mass-symmetric point, so in the expression for the correction between symmetrical and physical mass point, the leading term vanishes, e.g. for the mass of the η c meson

Charmonium spectrum
The new ensembles with the fine lattice spacings are very well suited for a study of charmonia, already at the coarsest lattice, the charmonium masses that we measure, neglecting disconnected contributions at the moment, are very close to their values in nature. In figure 5 we show the meson spectra of our ensembles A2 and B. We get a clear signal up to the J/ψ state and can extract reasonable plateau values for higher lying states summarized in table 6. We find good agreement with PDG [14] data because the states contain only charm valence quarks which in our simulations have their physical mass value. Further, we get a very precise result for the charmonium hyperfine splitting (m J/ψ − m η )/m η with perfect agreement to the experimentally known value 0.038, cited by the Particle Data Group, see table 7.
We also evaluate a number of other dimensionless quantities from ratios of flow observables and meson masses, summarized in table 8. For the latter and the charmonium masses in physical units we attempt continuum limit extrapolations via linear fits of the corresponding values from ensembles A2 and B. These results are also listed in tables 6-8.

Conclusions & outlook
We presented the scale setting and tuning of N f = 3 + 1 QCD using a massive renormalization scheme with a non-perturbatively determined clover coefficient from [1]. We produced two ensembles A1 and A2 with lattice sizes 96 × 32 3 and 128 × 48 3 at a lattice spacing a = 0.054 fm, see also [49], and an ensemble B on a 144 × 48 3 lattice at a finer lattice spacing a = 0.043fm.
As a first physics result, we measure the masses of the charmonium states η c , J/ψ, χ c 0 , χ c 1 and h c , which agree with their PDG values. In particular, we can reproduce the charmonium hyperfine splitting (m J/Ψ − m η )/m η within permille level precision of the experimentally known value 0.038. First continuum limit extrapolations of the measured quantities are attempted.
The next steps are to double the statistics of ensemble B and produce a third ensemble C on a 192 × 64 3 lattice at an even finer lattice spacing a = 0.032fm, in order to make more reliable continuum limit extrapolations. For the future, we plan to measure disconnected contributions to the charmonium masses and extract excited states. A further development of this project is the simulation of N f = 2 + 1 + 1 QCD close to the physical light quark masses. Another longterm   (5) 3.096900 (6) 3.4148 (3) 3.51066 (7) 3.52538(11)    goal is the determination of the strong coupling α S or equivalently the Λ-parameter in the case of four flavors. The currently most precise result comes from lattice simulations [8], see also [50]. One of the remaining uncertainty comes from the use of perturbative decoupling for the matching of the gauge couplings of QCD with four and three flavors. In [12] this uncertainty was estimated to be below 1.5% for the ratios of the Λ parameters which is still below the 3.5% precision of the Λ-parameter of [8]. But new techniques like [51] promise a larger precision for α S . The running of α S in four flavor QCD has been computed in [52]. Such a result, can be combined with the scale setting for N f = 3 + 1 QCD presented in this work to obtain α S fully non-perturbatively in four flavor QCD.
ens.   Table 9: Parameters of the algorithm: We give the forces used for the gauge and fermion fields with their integration levels, the integrators for the different levels and number of steps per trajectory resp. outer level, the hopping, c sw and twisted-mass parameters, the number of poles N p and the ranges used in the RHMC for strange and charm quarks, as well as the total length of the Markov chain and the acceptance rates.   Derivatives with respect to quark masses of pion/kaon, Dand D s -mesons and charmonium states η c , J/ψ, χ c 0 , χ c 1 and h c effective masses and flow observables on ensemble B.