Light hadrons from lattice QCD with light (u,d), strange and charm dynamical quarks

We present results of lattice QCD simulations with mass-degenerate up and down and mass-split strange and charm (N_f = 2+1+1) dynamical quarks using Wilson twisted mass fermions at maximal twist. The tuning of the strange and charm quark masses is performed at two values of the lattice spacing a~0.078 fm and a~0.086 fm with lattice sizes ranging from L~1.9 fm to L~2.8 fm. We measure with high statistical precision the light pseudoscalar mass m_PS and decay constant f_PS in a range 270<m_PS<510 MeV and determine the low energy parameters f_0, l_3 and l_4 of SU(2) chiral perturbation theory. We use the two values of the lattice spacing, several lattice sizes as well as different values of the light, strange and charm quark masses to explore the systematic effects. A first study of discretisation effects in light-quark observables and a comparison to N_f=2 results are performed.


Introduction and Main Results
The beginning of this century has assisted to radical improvements in theory, algorithms and supercomputer technology, leading to a far increased ability to solve non-perturbative aspects of gauge field theories in a lattice regularised framework.Following this path of improving the lattice setup, in this paper, we are reporting about our experiences and results when considering in addition to the u, d light dynamical flavours also the effects of the strange and charm sea quarks.By including a dynamical charm, we are now able to directly study its contribution to physical observables and to quantify the so far uncontrolled systematic effect present in lattice QCD simulations where the charm flavour in the sea is absent.
A number of different lattice fermion formulations are being used by several lattice groups, see refs.[1,2] for recent reviews.Here, we adopt a particular type of Wilson fermions, known as the Wilson twisted mass formulation of lattice QCD (tmLQCD), introduced in [3,4].This approach is by now well established, with many physical results obtained with two light degenerate twisted mass flavours (N f = 2) by our European Twisted Mass (ETM) Collaboration, see refs.[5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22].For a review see ref. [23].In the tmLQCD formulation a twisted mass term is added to the standard, unimproved Wilson-Dirac operator, and the formulation becomes especially interesting when the theory is tuned to maximal twist [4].The major advantage of the lattice theory tuned to maximal twist is the automatic O(a) improvement of physical observables, independently of the specific type of operator considered, implying that no additional, operator specific improvement coefficients need to be computed.Other advantages worth to mention are that the twisted mass term acts as an infrared regulator of the theory and that mixing patterns in the renormalisation procedure are expected to be simplified.
Detailed studies of the continuum-limit scaling in the quenched approximation [24][25][26][27] and with two dynamical quarks [7,10,17,28] have demonstrated that, after an appropriate tuning procedure to maximal twist, lattice artefacts not only follow the expected O(a 2 ) scaling behaviour [4], but also that the remaining O(a 2 ) effects are small, in agreement with the conclusions drawn in ref. [29].
The only exception seen so far is the neutral pseudoscalar mass, which shows significant O(a 2 ) effects.This arises from the explicit breaking of both parity and isospin symmetry, which are however restored in the continuum limit with a rate of O(a 2 ) as shown in [4] and numerically confirmed in refs.[17,30].Moreover, a recent analysis suggests that isospin breaking effects strongly affect only a limited set of observables, namely the neutral pion mass and kinematically related quantities [31,32].
In this paper we report on simulations with twisted mass dynamical up, down, strange and charm quarks.We realise this by adding a heavy mass-split doublet (c, s) to the light degenerate mass doublet (u, d), referring to this setup as N f = 2 + 1 + 1 simulations.This formulation was introduced in [33,34] and first explored in [35].As for the mass-degenerate case, the use of lattice action symmetries allows to prove the automatic O(a) improvement of physical observables in the non-degenerate case [33,34].First accounts of our work were presented at recent conferences [36,37].Recently, results with N f = 2 + 1 + 1 staggered fermions have been reported in [38][39][40].The inclusion of the strange and charm degrees of freedom allows for a most complete description of light hadron physics and eventually opens the way to explore effects of a dynamical charm in genuinely strong interaction processes and in weak matrix elements.
Here, we concentrate on results in the light-quark sector using the charged pseudoscalar mass m PS and decay constant f PS as basic observables involving up and down valence quarks only.In fig. 1 we show the dependence of (a) m 2 PS /2B 0 µ l and (b) f PS as a function of the mass parameter 2B 0 µ l , together with a fit to SU(2) chiral perturbation theory (χPT) at the smallest value of the lattice spacing of a ≈ 0.078 fm and lattice gauge coupling β = 1.95.We summarise the fit results for the low energy constants in table 1.These are the main results of this paper.PS /(2B 0 µ l ) and (b) the pseudoscalar decay constant f PS as a function of 2B 0 µ l fitted to SU(2) chiral perturbation theory, see table 1.The scale is set by the value of 2B 0 µ l at which the ratio f [L=∞] PS

/m
[L=∞] PS assumes its physical value [41] f π /m π = 130.4(2)/135.0(black star).The lattice gauge coupling is β = 1.95 and the twisted light quark mass ranges from aµ l = 0.0025 to 0.0085, see eq. ( 3) for its definition, corresponding to a range of the pseudoscalar mass 270 m PS 490 MeV.The kaon and D meson masses are tuned to their physical value, see table 4. The lightest point (open symbol) has not been included in the chiral fit, see the discussion in section 3.2.
A comparison between data obtained with N f = 2 + 1 + 1 and N f = 2 flavours of quarks -see sections 3.4 and 4, and ref. [17] -reveals a remarkable agreement for the results involving light-quark observables such as the pseudoscalar mass and decay constant or the nucleon mass.This provides a strong indication in favour of the good quality of our data in this new setup.In particular, barring cancellations due to lattice discretisation errors, these results would suggest that the dynamical strange and charm degrees of freedom do not induce large effects in these light-quark observables.In the N f = 2 case, data collected at four values of the lattice spacing have allowed us to properly quantify all systematic errors present in the determination of light-quark observables [17].In this first work with N f = 2 + 1 + 1 flavours, we consider data at two close values of the lattice spacing, while we defer to a forthcoming publication the inclusion of additional ensembles at a significantly lower lattice spacing and a more complete analysis of the systematic effects.
The rest of this paper is organised as follows.In section 2 we describe the gauge action and the twisted mass fermionic action for the light and heavy sectors of the theory.The realisation of O(a) improvement at maximal twist is also presented.In section 3 we define the simulation parameters, describe the tuning to maximal twist as well as the tuning of the strange and charm 0.0782(6) Table 1 Results of the fits to SU(2) χPT for the ensemble at β = 1.95.Predicted quantities are: the low energy constants l3,4 , the charged pseudoscalar decay constant in the chiral limit f 0 , the mass ratio 2B 0 µ l /m 2 PS at the physical point and the pion scalar radius r 2 NLO s .The first quoted error is from the chiral fit at β = 1.95, the second error is the systematic uncertainty that conservatively accommodates the best fitted central values of the three fits reported in table 9, section 4. The small error on the quoted lattice spacing comes exclusively from the fit at β = 1.95.The scale is set by fixing the ratio f 4(2)/135.0to its physical value [41].The chirally extrapolated Sommer scale r χ 0 is determined separately and not included in the χPT fits.For a comparison with the N f = 2 ETMC results, see [17].
quark masses and the relevance of discretisation effects.Section 4 includes a discussion of the fits to SU(2) χPT also for data on a slightly coarser lattice, a ≈ 0.086 fm, and provides a first account of systematic uncertainties.Our conclusions and future prospects are summarised in section 5.

Lattice Action
The complete lattice action can be written as where S g is the pure gauge action, in our case the so-called Iwasaki action [42,43], S l is the twisted mass Wilson action for the light doublet [3,4] and S h the one for the heavy doublet [33,34].
The choice of the gauge action is motivated by the non trivial phase structure of Wilson-type fermions at finite values of the lattice spacing.The phase structure of the theory has been extensively studied analytically, by means of chiral perturbation theory [44][45][46][47][48][49][50], and numerically [51][52][53][54][55][56].These studies provided evidence for a first order phase transition close to the chiral point for coarse lattices.This implies that simulations at non-vanishing lattice spacing cannot be performed with pseudoscalar masses below a minimal critical value.
The strength of the phase transition has been found [53,56] to be highly sensitive to the value of the parameter b 1 in the gauge action in eq. ( 2).Moreover, in [35] it was observed that its strength grows when increasing the number of flavours in the sea from N f = 2 to N f = 2 + 1 + 1, at otherwise fixed physical situation.Numerical studies with our N f = 2 + 1 + 1 setup have shown that the Iwasaki gauge action, with b 1 = −0.331,provides a smoother dependence of phase transition sensitive quantities on the bare quark mass than the tree-level-improved Symanzik [57,58] gauge action, with b 1 = −1/12, chosen for our N f = 2 simulations.
Another way to weaken the strength of the phase transition is to modify the covariant derivative in the fermion action by smearing the gauge fields.While the main results of this work do not use smearing of the gauge fields, we report in section 3.7 on our experience when applying a stout smearing [59] procedure, see also [60].

Action for the Light Doublet
The lattice action for the mass degenerate light doublet (u, d) in the so called twisted basis reads [3,4] where m 0,l is the untwisted bare quark mass, µ l is the bare twisted light quark mass, τ 3 is the third Pauli matrix acting in flavour space and is the massless Wilson-Dirac operator.∇ µ and ∇ * µ are the forward and backward gauge covariant difference operators, respectively.Twisted mass light fermions are said to be at maximal twist if the bare untwisted mass m 0,l is tuned to its critical value, m crit , the situation we shall reproduce in our simulations.The quark doublet χ l = (χ u , χ d ) in the twisted basis is related by a chiral rotation to the quark doublet in the physical basis where the twisting angle ω l takes the value |ω l | → π 2 as |m 0,l − m crit | → 0. We shall use the twisted basis throughout this paper.

Action for the Heavy Doublet
We introduce a dynamical strange quark by adding a twisted heavy masssplit doublet χ h = (χ c , χ s ), thus also introducing a dynamical charm in our framework.As shown in [34], a real quark determinant can in this case be obtained if the mass splitting is taken to be orthogonal in isospin space to the twist direction.We thus choose the construction [33,34] where m 0,h is the untwisted bare quark mass for the heavy doublet, µ σ the bare twisted mass -the twist is this time along the τ 1 direction -and µ δ the mass splitting along the τ 3 direction.
The bare mass parameters µ σ and µ δ of the non-degenerate heavy doublet are related to the physical renormalised strange and charm quark masses via [33] ( where Z P and Z S are the renormalisation constants of the pseudoscalar and scalar quark densities, respectively, computed in the massless standard Wilson theory.
A chiral rotation analogous to the one in the light sector transforms the heavy quark doublet from the twisted to the physical basis where the twisting angle ω h takes the value |ω h | → π 2 as |m 0,h − m crit | → 0.

O(a) improvement at maximal twist
One of the main advantages of Wilson twisted mass fermions is that by tuning the untwisted bare quark mass to its critical value, automatic O(a) improvement of physical observables can be achieved.
Tuning the complete N f = 2+1+1 action to maximal twist can in principle be performed by independently choosing the bare masses of the light and heavy sectors am 0,l and am 0,h , resulting, however, in a quite demanding procedure.
On the other hand, properties of the Wilson twisted mass formulation allow for a rather economical, while accurate alternative [4,34,35], where the choice am 0,l = am 0,h ≡ 1/2κ − 4 is made, and the hopping parameter κ has been introduced.
Tuning to maximal twist, i.e. κ = κ crit , is then achieved by choosing a parity odd operator O and determine am crit (equivalently κ crit ) such that O has vanishing expectation value.One appropriate quantity is the PCAC light quark mass [29,52,53] m PCAC = x ∂ 0 A a 0,l (x, t)P a l (0) 2 x P a l (x, t)P a l (0) where and we demand m PCAC = 0.For the quenched [25] and the N f = 2 case [17], this method has been found to be successful in providing the expected O(a) improvement and effectively reducing residual O(a 2 ) discretisation effects in the region of small quark masses [29].
The numerical precision required for the tuning of m PCAC to zero has been discussed in [8].Contrary to the N f = 2 case [5,8], where this tuning was performed once at the minimal value of the twisted light mass considered in the simulations, we now perform the tuning at each value of the twisted light quark mass µ l and the heavy-doublet quark mass parameters µ σ and µ δ .This obviously leaves more freedom in the choice of light quark masses for future computations.
Although theoretical arguments tell us that O(a) improvement is at work in our setup, a dedicated continuum scaling study is always required to accurately quantify the actual magnitude of O(a2 ) effects.In section 3.4 we provide a first indication that such effects are indeed small, at least for the here considered light meson sector; currently ongoing computations at a significantly smaller lattice spacing will allow for a continuum limit scaling analysis in this setup.
3 Simulation Details

Simulation Ensembles
We performed simulations at two values of the lattice gauge coupling β = 1.90 and 1.95, corresponding to values of the lattice spacing a ≈ 0.086 fm and a ≈ 0.078 fm, respectively.The parameters of each ensemble are reported in table 2. The charged pion mass m PS ranges from 270 MeV to 510 MeV.Simulated volumes correspond to values of m PS L ranging from 3.0 to 5.8, where the smaller volumes served to estimate finite volume effects, see table 3. Physical spatial volumes range from (1.9 fm) 3 to (2.8 fm) 3 .
As already mentioned, the tuning to κ crit was performed independently for each value of the mass parameters aµ l , aµ σ and aµ δ .The mass parameters of the heavy doublet aµ σ and aµ δ reported in table 2 are related to the strange and charm quark masses.In particular, they are fixed by requiring the simulated kaon and D meson masses to approximately take their physical values, as discussed in section 3.3.The simulation algorithm used to generate the ensembles includes in the light sector, a Hybrid Monte Carlo algorithm with multiple time scales and mass preconditioning, described in ref. [61], while in the strange-charm sector a polynomial hybrid Monte Carlo (PHMC) algorithm [62][63][64]; the implementation of ref. [65] is publicly available.
The positivity of the determinant of the Dirac operator is a property of the mass-degenerate Wilson twisted mass action, which does not necessarily hold in the non degenerate case for generic values of the mass parameters µ σ and µ δ . 1 The positivity is monitored by measuring the smallest eigenvalue λ h,min of Q † h Q h , where Q h = γ 5 τ 3 D h and D h is the Wilson Dirac operator of the nondegenerate twisted mass action in eq. ( 5).We observe that λ h,min is roughly proportional to the renormalised strange quark mass squared.Since we choose the mass parameters µ σ and µ δ such that the strange quark takes its physical value, a spectral gap in the distribution of Q † h Q h is observed, implying that the determinant of D h does not change sign during the simulation.While this is sufficient for the purpose of this study, we shall provide a detailed discussion of this issue in a forthcoming publication.2 Summary of the N f = 2 + 1 + 1 ensembles generated by ETMC at two values of the lattice coupling β = 1.90 and β = 1.95.From left to right, we quote the ensemble name, the value of inverse coupling β, the estimate of the critical value κ crit , the light twisted mass aµ l , the heavy doublet mass parameters aµ σ and aµ δ and the volume in units of the lattice spacing.Our notation for the ensemble names corresponds to X.µ l .L, with X referring to the value of β used.The run A100.24s is used to control the tuning of the strange and charm quark masses.3 For each ensemble, from left to right the values of m PCAC /µ l , m PS L, the integrated autocorrelation time of the plaquette, m PS and m PCAC in units of the trajectory length.Every ensemble contains 5000 thermalised trajectories of length τ = 1, except A40.24 which contains 8000 trajectories.
To generate correlators we use stochastic sources and improve the signal-tonoise ratio by using the "one-end trick", following the techniques also employed in our N f = 2 simulations [8].We have constructed all meson correlators with Fig. 2. The ratio m PCAC /µ l for the ensembles at β = 1.90 and 1.95 at the largest simulated volumes and as a function of 2B 0 µ l .For both ensembles the ratio m PCAC /µ l satisfies the 10% level criterion, except for the lightest point at β = 1.90 and β = 1.95 (open symbols), also affected by larger statistical errors.We assume Z A = 1, while the actual value Z A 1 can only improve all tuning conditions.local (L), fuzzed (F) and Gaussian smeared (S) sources and sinks.The use of smeared or fuzzed sources has stronger impact on the extraction of the kaon and D meson masses; results for the latter are reported in section 3.3, while a companion paper [66] discusses the adopted strategy for the less straightforward determination of these masses in the unitary N f = 2 + 1 + 1 Wilson twisted mass formalism.

Tuning to Maximal Twist
To guarantee O(a) improvement of all physical observables while also avoiding residual O(a2 ) effects with decreasing pion mass, the numerical precision of the tuning to maximal twist -quantified by the deviation from zero of m PCAC -has to satisfy |Z A m PCAC /µ l | µ l , µσ, µ δ aΛ QCD [5,8,17].The left-hand side contains the renormalised ratio of the untwisted mass over the twisted lightquark mass.A similar condition should be fulfilled by the error on this ratio.For the current lattice spacings, aΛ QCD ≈ 0.1, while the values of the axial current renormalisation factor Z A have not yet been determined.Nevertheless, since Z A enters as an O(1) multiplicative prefactor, and it is expected to be Z A 1 for our ensembles 2 , we adopt the conservative choice Z A = 1 in verifying the tuning condition.
Satisfying this constraint clearly requires a good statistical accuracy in the determination of the PCAC mass.The values of m PCAC /µ l reported in table 3 and shown in fig. 2 are well satisfying the tuning condition to maximal twist, with the exception of the lightest mass point at β = 1.90 and β = 1.95.We notice that the autocorrelation time of m PCAC reported in table 3 grows with decreasing values of the light quark mass µ l , thus rendering the tuning more costly for the two lightest points.For the ensemble B25.32, we are currently performing a new simulation aiming at a more accurate tuning to κ crit .We are also testing a reweighting procedure [36] in κ on the same ensemble, in view of applying it to the other not optimally tuned ensemble A30.32, and to future simulations.In what follows, we use the lightest mass points for consistency checks, and we exclude them from the final χPT fits.We also remind the reader that the small deviations from zero of am PCAC will only affect the O(a 2 ) lattice discretisation errors of physical observables [8].

Tuning of the Strange and Charm Quark Masses
The mass parameters µ σ and µ δ in the heavy doublet of the action in eq. ( 5) can in principle be adjusted so as to match the renormalised strange and charm quark masses by use of eq. ( 6).In practise, in this work, we fix the values of µ σ and µ δ by requiring that the simulated kaon mass m K and D meson mass m D approximately take their physical values.
A detailed description of the determination of the kaon and D meson masses is separately given in [66], while figures 3(a) and 3(b) show the resulting dependence of (2m 2 K − m 2 P S ) and m D upon the light pseudoscalar mass squared for both ensembles, and compared with the physical point.Table 4 summarises their numerical values, while the corresponding values for aµ σ and aµ δ are given in table 2. Observe also that, in order to be able to properly tune the strange and charm quark masses to their physical values, aµ σ must be chosen larger than aµ δ , since (see eq. ( 6)) the ratio Z P /Z S is significantly smaller than one [66].
While the kaon and D meson masses at β = 1.95 are sufficiently well tuned to their physical values, the ensembles at β = 1.90 with aµ δ = 0.190 carry a heavier kaon mass.The latter is instead visibly closer to its physical value for aµ δ = 0.197, as can be inferred from figure 3(a).We are currently performing simulations with aµ δ = 0.197 for other light quark masses.Moreover, another set of values of µ σ and µ δ are currently being used at β = 1.90 to generate ensembles with a slightly lower D meson mass and a third value of the kaon mass, in order to properly interpolate the lattice data to the physical strange quark mass.4 For each ensemble, the values of the kaon mass and the D meson mass as determined in [66].

Discretisation Effects in Light-quark Observables
In this section we explore discretisation effects in the analysed light-quark observables.To this aim we also make use of the determination of the chi- 1.5 1.6 1.7 q q q q q q q 0 5 10 15 20 25 30 β q 1.90 q 1.95 q 3.90 (m PS m N ) 2 (mPS fPS) 2 0 2 4 6 8 10 q q q q q q q 0.00 0.05 0.10 0.15 β q 1.90 q 1.95 q 3.90 rally extrapolated r 0 value for our data samples, as discussed in the following section 3.5.
In figures 4(a) and 4(b) we study the sensitivity of the charged pion mass and decay constant to possible discretisation effects, by comparing the N f = 2 + 1 + 1 data at β = 1.90 and β = 1.95 and the results obtained in twisted mass simulations with two dynamical flavours [17].The alignment of all data points at different values of β is in itself an indication of small discretisation effects.The comparison and good agreement with the N f = 2 data seems also to suggest no significant dependence upon the inclusion of dynamical strange and charm quarks for these light observables, at least at the present level of accuracy and provided that no cancellations occur due to lattice discretisation effects.However, only a more complete study at significantly different lattice spacings will allow to draw conclusions.
In the same spirit, we show in figure 5 an analogous ratio plot where the nucleon mass data points are included.The alignment of all data and the good extrapolation to the physical point is again evident.We defer to future publications the analysis of the baryon spectrum and the study of discretisation effects in strange-and charm-quark observables.

The Sommer Scale r 0
The Sommer scale r 0 [67] is a purely gluonic quantity extracted from the static inter-quark potential.Since the knowledge of its physical value remains rather imprecise, we use the chirally extrapolated lattice data for r 0 /a only as an effective way to compare results from different values of the lattice spacing.In this work, the lattice scale is extracted by performing χPT inspired fits to the very precise data for af PS and am PS , and by using the physical values of m π and f π as inputs.
Figures 6(a) and 6(b) display the data for r 0 /a at both values of the lattice coupling β = 1.90 and 1.95, and as a function of the bare lattice mass squared.The data are reasonably well described by a quadratic dependence, as also previously found for our N f = 2 ensembles.For a more detailed discussion of the possible functional forms and their theoretical interpretation see [37].To extrapolate to the chiral limit, we have performed fits using the largest available volume at each value of the pseudoscalar mass.The chirally extrapolated values for our N f = 2 + 1 + 1 ensembles are r χ 0 /a = 5.231(38) at β = 1.90 and r χ 0 /a = 5.710(41) at β = 1.95,where the lightest points of both ensembles have been excluded from the extrapolation, consistently with the fact that they do not satisfy our most stringent tuning condition to maximal twist.
In order to meaningfully compare the dependence upon the light quark mass at the two different lattice couplings β = 1.90 and 1.95, we estimated the slope of the functional form r 0 /r χ 0 = 1 + c r (r χ 0 m PS ) 4 , where the explicit lattice spacing dependence has been removed.We observe a mild dependence on the light quark mass and similar slopes c r [β = 1.90] = −0.0379(37) and c r [β = 1.95] = −0.0234(69).It is also worth noticing that the dependence upon the light quark mass of the N f = 2 + 1 + 1 data and that observed in the N f = 2 case [37] are not significantly different.

Effects of Isospin Breaking
A most delicate aspect of the twisted mass formulation is the breaking of the isospin symmetry.Clear evidence for this breaking has been found in the N f = 2 simulations by ETMC when comparing the neutral with the charged pion masses.Indeed, while the discretisation effects in the charged pion were observed to be very small, significant O(a 2 ) corrections appear when studying the scaling to the continuum limit of the neutral pion [17].Notice, however, that similar effects have not been observed in other quantities that are in principle sensitive to isospin breaking but not trivially related to the neutral pion mass.These observations are supported by theoretical considerations detailed in [31,32].
In the N f = 2 + 1 + 1 case, it turns out that the isospin breaking effect in the mass difference of charged and neutral pion masses is larger than for N f = 2 at fixed physical situation3 , as can be inferred from table 5. On the other hand, the same theoretical considerations as in [32] do apply to the case of N f = 2 + 1 + 1 flavours, and it is expected that the same class of physical observables as for N f = 2 will not be significantly affected by isospin  5 Measurements of the masses of the charged and the neutral pion.We compare runs at β = 1.95 and N f = 2 runs [17] with comparable lattice spacing and similar charged pion masses in physical units.All masses are reported in units of the chirally extrapolated r 0 for the same ensemble, see table 9, and r χ 0 /a = 5.316 (49) for N f = 2.We also report on the approximate value of c, giving the slope of the a 2 dependence of the pion mass splitting.
breaking corrections.Having said that, a careful measure of this effect for each observable or class of observables is anyway mandatory.The increase of the pion mass splitting with increasing the number of flavours in the sea is in line with the observation [35] of a stronger first order phase transition when moving from N f = 2 to N f = 2 + 1 + 1, as discussed in section 2.1.Indeed, the endpoint of the phase transition [44,45] corresponds to the critical value of the light twisted mass µ l,c where the neutral pion mass vanishes.The mass difference can be described by , where the coefficient c is related to µ l,c [44,45] and it is therefore a measure of the strength of the first order phase transition.Hence, a larger value of c means that simulations are to be performed at smaller values of the lattice spacing to reach, say, the physical point.Table 5 reports  We are currently performing simulations at a significantly different and lower lattice spacing than the present ensembles.They will allow to determine the slope c for N f = 2+1+1 more accurately and to better quantify the conditions to approach the physical point.

Stout Smeared Runs
In addition to our main simulation ensembles, we also performed runs with stout smeared gauge fields in the lattice fermionic action.The stout smearing as introduced in [59] was designed to have a smearing procedure which is analytic in the unsmeared link variables and hence well suited for HMC-type updating algorithms.In an earlier work with N f = 2 quark flavours [60] we showed that using smeared gauge fields in the fermion operator is reducing the strength of the phase transition in twisted quark mass simulations and  7 The masses in lattice units for the ensembles with one level of stout smearing.therefore allows to reach smaller quark masses at a given lattice spacing.
The definition of the stout smeared links can be found in [59], and for the parameter ρ connecting thin to fat gauge links we choose ρ = 0.15.In principle, such smearing can be iterated several times, with the price of rendering the fermion action delocalised over a larger lattice region.We made a conservative choice to maintain the action well localised and performed a single smearing step.As shown in [60], this kind of smearing does not substantially change the lattice spacing, and for the sake of comparison we thus kept the same value of β as in one of the non stout-smeared runs.On the other hand, the hopping parameter has to be tuned again, since the additive renormalisation of the quark mass is expected to be smaller.The parameters of our runs are given in Table 6.These runs have been done with the two-step polynomial Hybrid Monte Carlo (TS-PHMC) update algorithm [68].Results for the hadron masses are collected in Table 7, where the quoted errors include an estimate of the systematic error induced by variations of the fitting range.The method of estimating and combining statistical and systematic errors for the case of the kaon and D meson masses is described in [66].
As the values of m PCAC /µ l in table 7 show, the hopping parameters are well tuned to maximal twist.The masses in the run with smallest light twisted mass aµ l = 0.0040 (ensemble A st 40.24) satisfy r 0 m PS = 0.668 (10), r 0 m K = 1.315 (13) and r 0 m D = 4.25 (29).This means that the pion is lighter than in the corresponding run without stout smearing (see table 8) and the kaon and D meson masses are closer to their physical value.The smaller pion mass should be interpreted as due to a quark mass renormalisation factor closer to one.For the same reason the tuned twisted masses in the heavy doublet aµ σ = 0.170, aµ δ = 0.185 are smaller than in the runs without stout smearing.It is also interesting to compare the mass splitting of the charged and neutral pion between runs with and without stout smearing.For the ensemble A st 60.24 we obtain a neutral pion mass r χ 0 m 0 PS = 0.409 (34) and a charged pion mass r χ 0 m ± PS = 0.7861 (56), in units of the chirally extrapolated value r χ 0 /a = 5.280 (25), providing an estimate of the slope c = −12.6(0.8).Notice that the mass dependence of r 0 /a in table 6 is reduced as compared to the runs with no stout smearing, and a quadratic dependence on the bare quark mass has been used for the extrapolation to the chiral limit, consistently with the analysis of section 3.5.For the corresponding ensemble A60.24 without stout smearing, using data in tables 8 and 9, we obtain instead r χ 0 m 0 PS = 0.560(37), r χ 0 m ± PS = 0.9036 (71), and a slope c = −13.8(1.2),slightly but not significantly different from the stout-smeared case.
The runs with stout-smeared gauge links show somewhat better characteristics than the ones without stout smearing, but the improvements are not dramatic, at least with one level of stout smearing.More iterations would further accelerate the approach to lighter masses and are expected to further reduce the charged to neutral pion splitting.However, it is a delicate matter to establish how physical observables other than the spectrum will be affected.Based on these considerations and given the present pool of data, the final results in this study are obtained with non stout-smeared simulations.
4 Results: f PS , m PS and Chiral Fits We concentrate in this section on the analysis of the simplest and phenomenologically relevant observables involving up and down valence quarks.These are the light charged pseudoscalar decay constant f PS and the light charged pseudoscalar mass m PS .
The present simulations with dynamical strange and charm quarks, sitting at, or varying around, their nature given masses, should allow for a good measure of the impact of strange and charm dynamics on the low energy sector of QCD and the electroweak matrix elements.As a first step, one can determine the low energy constants of chiral perturbation theory (χPT).The values of af PS and am PS for our ensembles at β = 1.95 and β = 1.90 are summarised in table 8.In contrast to standard Wilson fermions, an exact lattice Ward identity for maximally twisted mass fermions allows for extracting the charged pseudoscalar decay constant f PS from the relation without need to specify any renormalisation factor, since Z P = 1/Z µ [3].We have performed fits to NLO SU (2) 8 Lattice measurements of the charged pseudoscalar mass am PS , the charged pseudoscalar decay constant af PS and the Sommer scale in lattice units r 0 /a for our two ensembles at β = 1.90 (A set) and β = 1.95 (B set).The value of the light twisted mass aµ l and the spatial length L/a are also shown.Quoted errors are given as (statistical)(systematic), with the estimate of the systematic error coming from the uncertainty related to the fitting range.
We thus simultaneously fit our data for the pseudoscalar mass and decay constant to the following formulae, where the contributions F , D and T parametrising finite size corrections, discretisation effects and NNLO χPT effects, respectively, will be discussed below: with the pseudoscalar mass squared at tree level defined as χ µ ≡ 2 B 0 µ l and the chiral expansion parameter by ξ ≡ χ µ / (4πf 0 ) 2 .The low energy constants l 3 and l 4 receive renormalization corrections according to li = l i + ln [Λ 2 /χ µ ], with Λ the reference scale.During the fitting procedure, where all quantities are defined in lattice units, we set the reference scale to a single lattice spacing to let its constant logarithmic contribution vanish.Once the scale of the simulation has been set, the low energy constants are rescaled to the scale of the physical pion mass to recover the physical values l3 and l4 .
Systematic errors can arise from several sources: finite volume effects, neglecting of higher orders in χPT and finite lattice spacing effects.These different corrections are accounted for explicitly in eq.(11).Finite volume corrections are described by the rescaling factors denoted by F m 2 PS and F f PS , computed in the continuum theory.Notice that the discretisation effects present in the neutral pion mass, see section 3.6, generate peculiar finite volume corrections which have been recently analysed in ref. [69].We shall comment on them later.We investigated the effectiveness of one loop continuum χPT finite volume corrections, as first computed in [70], which do not introduce any additional low energy constants.However, the resummed expressions derived by Colangelo, Dürr and Haefeli (CDH) in [71] describe the finite volume effects in our simulations better, be it at the expense of the introduction of two new free parameters, and are thus adopted for this analysis.To O(ξ 2 ), these corrections read with geometric contributions defined as The K i are the modified Bessel functions and the low energy constants l 1 and l 2 again receive renormalisation corrections.Equations ( 12) and ( 13) use the shorthand notation λ n = √ nm PS L. The ρ n in eq. ( 12) are a set of multiplicities, counting the number of ways n 2 can be distributed over three spatial directions 4 .Because the finite volume corrections in the case of the volumes used in the chiral fits are fairly small to begin with and subsequent terms quickly decrease, the sums over n can be truncated rather aggressively without real loss of precision.It is therefore unnecessary, in practise, to go beyond the lowest contributions.The parameters l 1 and l 2 , which are in fact low energy constants appearing at NLO in χPT, cannot be determined well from the small finite volume corrections alone.Priors are therefore introduced as additional contributions to the χ 2 , weighting the deviation of the parameters from their phenomenological values by the uncertainties in the latter.The values used as priors are -0.4(6) for l1 and 4.3(1) for l2 [71], as reported in table 9. We used the largest available volumes for each ensemble, in the χPT fits.For those points, the difference between the finite volume and the infinite volume values estimated via CDH formulae for f PS and m 2 PS are within 1%, except for the runs B85.24 and A60.24 (see table 2 and table 8), where they are about 1.5% for both quantities.
Because of the automatic O(a) improvement of the twisted mass action at maximal twist, the leading order discretisation artefacts in the chiral formulae of (11) are at least of O(a 2 ), and O(a 2 µ) for m 2 PS .The mass and decay constant of the charged pion have been studied up to NLO [44,45,50] in the context of twisted mass chiral perturbation theory (tmχPT).The regime of quark masses and lattice spacings at which we have performed the simulations is such that µ l aΛ 2 QCD .In the associated power counting, at maximal twist, the NLO tmχPT expressions for the charged pion mass and decay constant preserve their continuum form.The inclusion of the terms proportional to D m 2 PS ,f PS , parametrising the lattice artifacts in eq. ( 11), represents an effective way of including sub-leading discretisation effects appearing at NNLO.The finite lattice spacing artefacts can of course not be determined using only data from a single lattice spacing.In addition, including these terms when analysing data with an insufficient range in a, may lead to mixing of these degrees of freedom with continuum parameters and thereby destabilise the fits.Hence, these terms were neglected for the separate fits, but included to arrive at a qualitative estimate of these systematic effects in a combined fit to the data at both lattice spacings.
Finite size effects on our data at finite lattice spacing can be analysed in the context of twisted mass chiral perturbation theory as recently proposed in ref. [69]. 5However, our present limited set of data with only a small number of different volumes all of them at a single value of the lattice spacing, is not sufficient to apply such an analysis.We plan, however, to perform dedicated runs on different volumes to confront our data to the finite size effect formulae of ref. [69] and to estimate in particular the size of the pion mass splitting in this alternative way.
Two new parameters k m and k f enter these corrections.Again, a limited range of input pion masses may lead to poorly constrained values of these newly introduced parameters, some degree of mixing among different orders and fit instabilities.To retain predictive power and stability, additional priors are given for k m and k f , both priors set to 0(1), analogously to what is done for l 1 and l 2 in the CDH finite volume corrections.
To set the scale at each lattice spacing, we determine aµ phys , the value of aµ l at which the ratio m 2 PS (L = ∞)/f PS (L = ∞) assumes its physical value.We can then use the value of f PS , or equivalently m PS , to calculate the lattice spacing a in fm from the corresponding physical value.We also perform a chiral fit combining the two different lattice spacings.With only two different values of β, that are in fact fairly close to each other, a proper continuum limit analysis cannot be performed.Instead, we treat this combined fit as a check on the presence of lattice artefacts and the overall consistency of the data.Without a scaling variable, such as the Sommer scale r 0 , the data from different lattice spacings cannot be directly combined.Rather, the ratios of lattice spacings and light quark mass renormalisation constants (Z µ = 1/Z P ), as well as the renormalised B 0 parameter are left free in the fit.
In order to estimate the statistical errors affecting our fitted parameters, we generate at each of the µ l values 1000 bootstrap samples for m PS and f PS extracted from the bare correlators, organised by blocks.For each sample, and combining all masses, we fit m 2 PS and f PS simultaneously as a function of µ l .The parameter set from each of these fits is then a separate bootstrap sample for the purposes of determining the error on our fit results.By resampling f PS and m PS on a per-configuration basis, correlations between these quantities are taken into account.
Our final results for the separate and combined fits are summarised in table 9.The χPT fit ansätze provide a satisfactory description of the lattice data, with a χ 2 /d.o.f = 5.68/3 1.9 at β = 1.95, χ 2 /d.o.f = 4.31/5 0.9 at β = 1.90, and 16.9/11 1.5 for the combined fit.We also predict the scalar radius of the pion at next to leading order The numerical values in table 9 for the combined fit show a very good agreement with the results from the separate fits, and with errors at the percent level throughout.The fits for f PS and m PS at β = 1.95 are displayed in fig- -0.07820(59) 0.07775(39) -Table 9 Results of the fits to SU(2) χPT for the ensembles at β = 1.95 and β = 1.90, separate and combined.The largest available volumes are used for each ensemble.Predicted quantities are: the low energy constants l3,4 (while l1,2 are introduced with priors), the charged pseudoscalar decay constant in the chiral limit f 0 , the mass ratio 2B 0 µ l /m 2 PS at the physical point and the pion scalar radius r 2 NLO s .The scale is set by fixing the ratio f 4(2)/135.0to its physical value [41].The chirally extrapolated Sommer parameter r χ 0 is determined separately and not included in the chiral fits.For a comparison with the N f = 2 ETMC results, see [17].The data presented here do not allow yet for a complete account of the systematic effects, but we extract estimates of their magnitude by extending the fits with additional terms as written down in eq.(11).Checks were done for χPT NNLO terms and O(a 2 ) corrections separately.Including NNLO corrections does not lower the total χ 2 of the fit, while we do observe a shift of several standard deviations for the lower order parameters already present in the NLO fit.Using these shifted values to obtain the implied NLO approximation produces fits with much larger values of χ 2 .We conclude that the current data lack the precision and range in quark masses to constrain NNLO effects, the added degrees of freedom mix with NLO effects and destabilise the fit instead.In practise, we conclude that the systematic error from the truncation of χPT is unobservable at the current level of precision.Inclusion of O(a 2 ) corrections leads to similar observations, as the difference between the lattice spacings and the statistical accuracy of the data is too small to result in a stable fit.The fit mixes D f PS and D m 2 PS on the one hand and f 0 , B 0 and the rescaling in the lattice spacing and the quark mass on the other.PS /2B 0 µ l and (b) the pseudoscalar decay constant f PS as a function of 2B 0 µ l , for the ensemble at β = 1.90, fitted to SU(2) chiral perturbation theory, eq. ( 11).The scale is set by aµ phys , the value of aµ l at which the ratio f [L=∞] PS

/m
[L=∞] PS assumes its physical value [41] f π /m π = 130.4(2)/135.0(black star).The light twisted masses used in the fit range from aµ l = 0.004 to 0.010.The lightest point (open symbol) lies outside our most conservative tuning criterion to maximal twist, and is not included in the fit.
The chirally extrapolated Sommer scale r χ 0 has been determined separately, using a fit of r 0 /a with quadratic dependence on the bare light quark mass, as shown in figures 6(a) and 6(b), and using the lattice spacing determined by the chiral fits.As also reported in table 9, the obtained values are r χ 0 = 0.4491 (43) fm at β = 1.90 and r χ 0 = 0.4465(48) fm at β = 1.95,where only statistical errors are quoted.For consistency, we also verified that a combined chiral fit with the inclusion of r 0 /a, as data points and additional fit parameter, gives results anyway in agreement with the strategy adopted here.
For our final estimates of the low energy constants l3,4 and the chiral value of the pseudoscalar decay constant f 0 we use the predictions from the β = 1.95 ensemble based on two important observations.First, the strange quark mass in this ensemble is better tuned to the physical value.Secondly a reduced isospin breaking is observed at this finer lattice spacing.The results for the β = 1.90 ensemble and the combined fits serve instead as an estimation of systematic uncertainties.As a result of the current N f = 2 + 1 + 1 simulations we thus quote l3 = 3.70(7)(26) l4 = 4.67(3)(10) , and f 0 = 121.14(8)(19) MeV, where the first error comes from the chiral fit at β = 1.95, while the second quoted error conservatively accommodates the central values from the β = 1.90 and combined fits as a systematic uncertainty.The predictions for l3 and l4 are in good agreement and with our two-flavour predictions [17] and with other recent lattice determinations [2,72].

Conclusions and Outlook
In this paper we have presented the first results of lattice QCD simulations with mass-degenerate up, down and mass-split strange and charm dynamical quarks using Wilson twisted mass fermions at maximal twist.This study constitutes a first step in our effort to describe low energy strong dynamics and electroweak matrix elements by fully taking into account the effects of a strange and a charm quark.
We have considered ensembles at slightly different lattice spacings simulated with Iwasaki gauge action at β = 1.95 with a ≈ 0.078 fm and β = 1.90 with a ≈ 0.086 fm.The charged pseudoscalar masses range from 270 to 510 MeV and we performed fits to SU(2) chiral perturbation theory with all data at a value of m PS L 4. This analysis provides a prediction for the low energy constants l3 = 3.70(7)(26) and l4 = 4.67(3) (10), for the charged pseudoscalar decay constant in the chiral limit f 0 = 121.14(8)(19) MeV and for the scalar radius at next-to-leading order r 2 NLO s = 0.724(5)(23) fm 2 .A companion paper [66] describes the less straightforward determination of the kaon and D-meson masses for the same ensembles.
We have compared our results in the light meson sector with those obtained for N f = 2 flavours of maximally twisted mass fermions, ref. [17].There, an extrapolation to the continuum limit, a study of finite size effects and checks against higher order χPT have been performed, leading to a controlled determination of systematic errors.The comparison we have carried through does not show any significant difference between N f = 2 and N f = 2 + 1 + 1 flavours, at least at the present level of accuracy.These results would suggest that effects of the strange and charm quarks are suppressed for these light observables, as it should be expected.The same comparison has also been used for a first investigation of lattice discretisation errors.As figures 4(a) and 4(b) show, the N f = 2 + 1 + 1 data are completely consistent with the corresponding ones obtained for N f = 2, where the discretisation effects have turned out to be very small.Thus, it can be expected that also for the case of N f = 2 + 1 + 1 flavours the lattice spacing effects will be small, at least for the light meson sector considered here.Notice however that, at the present level of accuracy, there is still the possibility that cancellations occur between physical contributions due to dynamical strange and charm quarks and lattice discretisation effects.A more accurate study at a significantly lower lattice spacing will allow to draw conclusions.
One aspect of the twisted mass formulation is the breaking of isospin symmetry.Its effect is likely to be most pronounced in the lightest sector, where lattice discretisation effects at O(a 2 ), affecting the neutral pseudoscalar mass only, generate a mass splitting between the charged and the neutral pseudoscalar mesons.While this mass splitting for N f = 2 + 1 + 1 flavours has been found here to be larger than in the N f = 2 simulations at fixed physical situation, we do not find further effects in other quantities computed so far.This observation is supported by theoretical arguments [31,32] and consistent with our experience in the N f = 2 flavour case.
We consider the present results to be encouraging to proceed with the N f = 2 + 1 + 1 flavour research programme of ETMC.In particular, we want to perform the non-perturbative renormalisation with dedicated runs for N f = 4 mass-degenerate flavours, an activity which we have started already.Furthermore, we want to compute the quark mass dependence of many physical quantities towards the physical point where the pion assumes its experimentally measured value.We are currently performing simulations at a significantly different and lower lattice spacing than the present ensembles.Both strategies, smaller quark masses and smaller lattice spacings, will allow us to estimate systematic effects on a quantitative level and to obtain in this way accurate physical results in our N f = 2 + 1 + 1 flavour simulations with statistical and systematical errors fully under control.

Fig. 1 .
Fig. 1.(a) The charged pseudoscalar mass ratio m 2PS /(2B 0 µ l ) and (b) the pseudoscalar decay constant f PS as a function of 2B 0 µ l fitted to SU(2) chiral perturbation theory, see table 1.The scale is set by the value of 2B 0 µ l at which the ratio f

Fig. 3 .
Fig. 3. (a): 2m 2 K − m 2 PS , and (b): m D , as a function of m 2 PS , for β = 1.95 (blue) and β = 1.90 (orange).The physical point is shown (black star).The kaon and D meson masses appear to be properly tuned at β = 1.95.The ensembles at β = 1.90, µ δ = 0.190 have a larger value of the strange quark mass, while the red point at β = 1.90, aµ δ = 0.197 appears to be well tuned.Data points have been scaled with the lattice spacing a = 0.08585(53) fm for β = 1.90, and a = 0.07820(59) fm for β = 1.95, obtained in this work and where the errors are only statistical.

Fig. 6 .
Fig.6.The Sommer scale r 0 /a as a function of (aµ l ) 2 for (a) β = 1.90 and (b) β = 1.95.The lines represent a linear extrapolation in (aµ l ) 2 to the chiral limit.The lightest point (open symbol) is not included in the fits and we have always used the largest available volume for a given value of the mass.
on the values of m ± PS , m 0 PS and c for some examples taken from the β = 1.95 ensemble and the N f = 2 ensemble with the closest values of the lattice spacing and physical charged pseudoscalar mass.As anticipated, the coefficient c increases in absolute value from N f = 2 to N f = 2 + 1 + 1.

ures 1 (
a) and (b), while in figures 7(a) and (b) we show the analogous fits at β = 1.90.Figures 8(a) and (b) show the results for the fit combining the two β values.

Fig. 7 .
Fig. 7. (a) The charged pseudoscalar mass ratio m 2PS /2B 0 µ l and (b) the pseudoscalar decay constant f PS as a function of 2B 0 µ l , for the ensemble at β = 1.90, fitted to SU(2) chiral perturbation theory, eq.(11).The scale is set by aµ phys , the value of aµ l at which the ratio f

Fig. 8 .
Fig. 8. (a) The charged pseudoscalar mass ratio (m PS /2B 0 µ l ) 2 and (b) the pseudoscalar decay constant f PS as a function of 2B 0 µ l , for the combined ensembles at β = 1.90 and β = 1.95, and fitted to eq. (11).The scale is set as in figure 7 (black star).The light twisted masses used in the fit range from aµ l = 0.0035 to 0.010.The lightest point at β = 1.90 (open orange symbol) and at β = 1.95 (open blue symbol) lie outside our most conservative tuning criterion to maximal twist, and are not included in the fit.
Ensemble m PCAC /µ l m PS L τ int ( P ) τ int (am PS ) τ int (am PCAC ) Parameters of the runs with stout smearing on L/a = 24, T /a = 48 lattices.The number of thermalised trajectories with length τ = 1 is given by N traj. .The label "st" in the ensemble name refers to the use of stout smearing, compared to the non stout-smeared ensemble in table 2.