Simulating twisted mass fermions at physical light, strange and charm quark masses

We present the QCD simulation of the first gauge ensemble of two degenerate light quarks, a strange and a charm quark with all quark masses tuned to their physical values within the twisted mass fermion formulation. Results for the pseudoscalar masses and decay constants confirm that the produced ensemble is indeed at the physical parameters of the theory. This conclusion is corroborated by a complementary analysis in the baryon sector. We examine cutoff and isospin breaking effects and demonstrate that they are suppressed through the presence of a clover term in the action.


Introduction
Simulations of Quantum Chromodynamics directly with physical quark masses, large enough volume and small enough lattice spacing have become feasible due to significant algorithmic improvements and availability of substantial computational resources. In fact, state-of-the-art simulations using different discretization schemes are currently being carried out worldwide.
Within the twisted mass formulation [1,2,3], the European Twisted Mass Collaboration (ETMC) has carried out simulations directly at the physical value of the pion mass [4,5] with N f = 2 mass-degenerate up and down quarks at a lattice spacing of a = 0.0913 (2)fm. This is a remarkable result, since explicit isospin breaking effects associated with twisted mass fermions can make physical point simulations at too coarse values of the lattice spacing very difficult. Being able to reach the physical pion mass for the case of N f = 2 flavours was therefore of great importance and many physical quantities have already been computed on the so generated gluon field configurations. Examples are meson properties [5,4,6,7,8,9], the structure of hadrons [4,10,11,12,13,14] and the anomalous magnetic moment of the muon [5].
The success of these N f = 2 flavour simulations strongly suggests to extend the calculation by adding the strange and charm quarks as dynamical degrees of freedom, a situation we will refer to as N f = 2 + 1 + 1 simulation. Adding a quark doublet is a natural step for twisted mass fermions. However, it is known that the presence of a heavy quark doublet in the sea gives rise to larger discretization effects than having only the light up and down quarks.
This paper reports on our successful, but demanding tuning effort to reach a physical situation with the first two quark generations tuned to their physical values in the twisted mass representation. We will present first results for low-lying meson masses and decay constants as well as baryon masses. In addition, we describe a comprehensive determinations of the lattice spacing from the meson and baryon sectors as well as from gradient flow observables. Furthermore, we discuss isospin breaking effects of twisted mass fermions in the neutral and charged pion and in the ∆ sector. Demonstrating the successful generation of an N f = 2 + 1 + 1 ensemble of maximally twisted mass fermions at physical quark masses is the essential result of this paper that lays the ground for a future very rich research program within the twisted mass formulation with an eventually large impact for ongoing and planned experiments.
The outline of the paper is as follows: In section 2 we introduce the employed twisted mass action and discuss details of the parameters used in the Hybrid Monte Carlo simulation. In section 3 we discuss our tuning procedure to reach physical light, strange and charm quark masses, which includes tuning for O(a) improvement and a discussion on the isospin splitting. In section 4 we present mesonic quantities for our ensemble, including a determination of the lattice spacing via the pion decay constant and heavy quark observables. In section 5 we discuss nucleon properties including the determination of the lattice spacing via the nucleon mass. In addition, we discuss possible isospin splitting in the ∆-baryon sector. In section 6 we summarize the different determination of lattice spacing via gluonic, mesonic and baryonic observables and conclude.

Action
We employ the twisted mass fermion formulation, within which observables are automatically O(a) improved when working at maximal twist [2,15]. This formulation has proven very advantageous: It allows one to perform safe, infrared regulated simulations and simplified renormalization in some cases. There is no need for improvement on the operator level due to automatic O(a)-improvement and cut-off effects turn out to be relatively small except for the special case of the neutral (unitary) pion mass.  The action of twisted mass fermions is given by where we choose the Iwasaki improved gauge action for S g [16] which reads with the bare inverse gauge coupling β = 6/g 2 0 , b 1 = −0.331 and b 0 = 1 − 8b 1 . In the case of the light up and down quark doublet, the action takes the form Here, χ = (u, d) t represents the light quark doublet, µ is the twisted and m the (untwisted) Wilson quark mass. The Pauli matrix τ 3 acts in flavour space and D W is the massless Wilson-Dirac operator. Note that the Wilson quark mass m and the clover term i For the heavy quark action, with mass non-degenerate strange (s) and charm (c) quarks, we construct a quark doublet χ h = (s, c) t for which the action reads [15] The important addition compare to eq. (3) is the term µ δ τ 1 with τ 1 again acting in flavour space. The Wilson quark masses in eqs. (3,4) are related to the hopping parameter κ as m = 1/2κ − 4. By tuning the light Wilson bare quark mass m to its critical value m crit the maximally twisted fermion action is obtained for which all physical observables are automatically O(a)-improved [2,15]. Setting m h = m = m crit this property takes over to the heavy quark mass action such that only one bare mass parameter has to be tuned to its critical value which is a great simplification for practical simulations. However quadratic lattice artefacts can be sizable but by introducing a clover term they can be suppressed, e.g. in case of the neutral pion mass as shown in [18,19,5,20]. Here, the clover parameter is set by using an estimate from 1-loop [21] tadpole boosted perturbation theory given by with P the plaquette expectation value. For our target parameter set, shown in tab. 1, the plaquette expectation value is given by P = 0.554301 (6), which is consistent with setting c SW = 1.69.

Algorithm
For the generation of the gauge field configurations we use as a basis the Hybrid Monte Carlo (HMC) algorithm [22,23] as described in Ref. [24,25]. For the light quark sector Hasenbusch mass preconditioning [26,27] is applied. In particular, we employ four determinant ratios with mass shifts ρ = {0.0; 0.0003; 0.0012; 0.01; 0.1}. The heavy quark determinant is treated by a rational approximation [28,29] with 10 terms tuned such that the (eigenvalue) interval [0.000065, 4.7] is covered. For the molecular dynamics integration we use a nested second order minimal norm integrator. This results in 12 integration steps for the smallest mass term in the light and heavy quark sector and 192 steps for the gluonic sector [20]. We use the software package tmLQCD [25] which incorporates the multi-grid algorithm DDalphaAMG for the inversion of the Dirac matrix [30]. The force calculation in the light quark sector is accelerated by a 3-level multi-grid approach optimized for twisted mass fermions [31]. Moreover, we extended the DDalphaAMG method for the mass non-degenerate twisted mass operator. The multi-grid solver used in the rational approximation [32,33] is particularly helpful for the lowest terms of the rational approximation, as well as for the rational approximation corrections in the acceptance steps, where it yields a speed up of two over the standard multi-mass shifted conjugate gradient(MMS-CG) solver. We checked the size of reversibility violation of this setup yielding a standard deviation < 0.01 for δ∆H and |1 − ∆H | < 0.02 fulfilling the criteria discussed in [34]. Here, δ∆H is the difference of the Hamiltonian at integration time t = 0 and the Hamiltonian of the reversed integrated field variables after one trajectory is performed.

Tuning of the light quark sector
As shown in Ref. [35,36] a most suitable and theoretically sound condition for the desired automatic O(a) improvement for twisted mass fermions is achieved by demanding a vanishing of the partially conserved axial current (PCAC) quark mass  Table 2: Summary of the parameters of the ensembles used for the tuning and final runs: L is the lattice spatial size with the time direction taken to be 2L, aµ is the twisted mass parameter of the mass degenerate light quarks, κ is the hopping parameter (common to all flavours), N th are the number of thermalized trajectories in molecular dynamics units (MDU), aµ σ and aµ δ are the bare twisted mass parameter of the mass non-degenerate fermion action used for the heavy quark sector. The ensembles cB211.072.64.r1 and cB211.072.64.r2 represent the targeted large volume runs at the physical point.
with A a µ the axial vector current and P a the pseudoscalar current. In the twisted basis and for light, mass degenerate quarks, the axial and pseudoscalar currents can be calculated via using τ + = (τ 1 + iτ 2 )/2 where τ i are the Pauli matrices. The tuning procedure to maximal twist requires a value of the hopping parameter κ = κ crit where m PCAC (κ crit ) = 0. Note that the corresponding definition of the critical mass am crit = 1/(2κ crit )−4 is a function of aµ , aµ σ , aµ δ . Thus, even if the 1/a divergence in m crit is independent from µ , µ σ and µ δ , determining am crit at the µ σ and µ δ values of interest is important in order to keep lattice artifacts small which are introduced by the heavy quark doublet [37,38]. Instead the dependence of am crit on µ reflects much milder discretization errors. In practice, we allow for some tolerance to this strict condition and following Ref. [39] we impose that within errors. In eq. (7) Z A is the renormalization constant of the axial current. Fulfilling the condition eq. (7) is numerically consistent with O(a)-improvement of physical observables, where it entails only an error of order O((Z A ·m PCAC /µ ) 2 ). Hence for < 0.1 follows for the targeted lattice spacing the error is comparable to other O([aΛ QCD ] 2 ) discretization errors. This allows an O(a 2 ) scaling of physical observables towards the continuum limit. In order to tune to κ crit , we have generated several ensembles with fixed volumes of size 24 3 · 48 and 32 3 · 64, as listed in Table 2. For a fixed twisted mass parameter of the up and down doublet, we scan over several values of the hopping parameter κ, see Table 2. After fixing κ crit in this manner we proceed by tuning the light and heavy twisted mass parameters to realize physical pion, kaon and D-meson masses and decay constants. This procedure, which is described in more detail below will provide the input parameters for the target large volume simulations, denoted as the ensembles cB211.072.64.r1 and cB211.072.64.r2 in Table 2.
Initially, we had attempted to start our N f = 2+1+1 simulations at a smaller value of β = 1.726 that would correspond to the lattice spacing of our N f = 2 ensemble with a ∼ 0.095 fm [5]. However, it turned out that tuning to maximal twist for a physical value of the pion mass for this β-value was not feasible. Nevertheless, our simulations at β = 1.726 for pion masses in the range between 170 MeV and 350 MeV allowed us to develop a tuning strategy to realize the situation of maximal twist and also to reach the physical kaon and D-meson masses. This tuning strategy was then used at the finer lattice spacing as discussed in the present paper. The occurrence of instabilities of the simulations at β = 1.726 when approaching the physical pion mass is, in fact, not unexpected. With twisted mass fermions, going to sufficiently small values of the light twisted mass parameter at a fixed lattice spacing one either enters the Aoki [40] or the Sharpe-Singleton [41] regime, see for an recent overview [42]. For the Sharpe-Singleton case, which is realized in our unquenched simulations, a sizable O(a 2 ) negative shift of the neutral pion mass occurs.
Let us consider the region close to maximal twist, where |ω − π/2| 1 or, equivalently, m = m 0 − m crit µ . Here the pion mass splitting can be related the PCAC quark mass by [40,43] am PCAC ∼ Zam m 2 π m 2 where Z = Z m Z P /Z A is a combination of the untwisted quark mass (Z m ), the pseudoscalar (Z P ) and the axial (Z A ) renormalization factors and am 0 denotes the bare quark mass. The charged pion mass is denoted throughout this paper by m π , while the neutral pion is given by m π (0) . The twisted mass angle ω can be defined via the gap equation, see [40,43]. From eq. (8) it is clear that the tuning necessary to satisfied eq. (7) becomes very hard for a large pion mass difference m 2 π − m π (0) 2 0. In the Sharpe-Singleton scenario a first order phase transition is predicted from chiral perturbation theory. In simulations on finite lattices this leads to large fluctuations and jumps of physical observables [44,45,46,47,48,49], driving the simulations to become unstable. This makes it very hard to tune successfully to maximal twist. In our simulations at β = 1.726 we observed a strong dependence of the PCAC quark mass on the bare mass parameter m 0 , which made it difficult to tune to the critical hopping parameter for a pion mass below 170 MeV. Although at β = 1.726 we did not investigate in detail which of the lattice ChPT scenario is realized, the fact that at β = 1.778 we find (see Section 3.3) a neutral pion mass ∼ 20% smaller than the charged one suggests that a Singleton-Sharpe lattice scenario occurs in the scaling region with our chosen action (see Sec. 2).
In order to avoid the aforementioned difficulties, we therefore decided to choose a finer value of the lattice spacing that would facilitate tuning to critical mass at the physical point. We found that a value of β = 1.778, corresponding to a ≈ 0.08 fm, allows us to tune to maximal twist successfully. In the following, we consider therefore a lattice volume of size 64 3 · 128, which is sufficiently large to suppress finite size effects but at the same time can be simulated with reasonable computational resources, given the algorithmic improvements that were discussed in section 2.1. For the tuning process of κ, which is a function of the light, strange and charm quark mass parameters, we use the 24 3 · 48 and 32 3 · 64 lattices, see section 3.2. for more details. The dependence of the PCAC quark mass on κ at fixed light twisted mass parameter is shown in Fig. 1. Note that it can be assumed that eq. (8) is valid here for the range −0.4141 am −0.4135 i.e. |a(m − m crit )| < 0.0003. Using simple linear fits for the L = 32 ensembles, we determine a critical value of κ, κ crit = 0.1394265. We then employ this κ-value for our large volume ensembles.
For the simulations on the 64 3 · 128 lattices we first thermalize one configuration using 500 trajectories. We then use this configuration as a starting point for two replicas, each having a final statistics of about 1500 MDUs. In Fig. 1 we depict the Monte Carlo history of the PCAC quark mass for these two replicas, where we show, for better visibility, one history plotted by reversed history. The PCAC quark mass fluctuates around zero and does not show particularly large autocorrelation times nor any indication of a first order Sharpe-Singleton transition. Performing the average over the two replica runs, we find m PCAC /µ = 0.03 (2). Thus, the condition of eq. (7) is nicely fulfilled. Note that here we do not include the renormalization factor Z A . However, our first estimate is that Z A ≈ 0.8 and anyhow smaller than one, making the condition even better fulfilled. We therefore conclude that the tuning to maximal twist is achieved for the N f = 2 + 1 + 1 setup. And, as we will demonstrate below, the parameters of the cB211.072.64 runs are chosen such that we indeed simulate at, or very close to the physical values of the pion, the kaon and the D-meson masses.

Tuning of the heavy quark sector
In tuning the mass parameters of the heavy quark sector we exploit the fact that the value of the critical hopping parameter, as determined in the light quark sector, can be employed also for the heavy quark action while preserving automatic O(a)improvement of all physical observables [3,50]. Nevertheless, tuning the heavy twisted mass parameters to reproduce the physical values of the strange and charm quark masses is a non-trivial task, owing to the O(a 2 ) flavour violation [15] inherent to the heavy sector fermion action in eq. (4). In order to tackle the problem, it is convenient to employ in an intermediate step the so-called Osterwalder Seiler (OS) fermions [51] in the valence which avoids these mixing effects. The OS-fermions can be used in a well defined mixed action setup as valence fermions at maximal twist with the same critical mass, m crit , as determined in the unitary setup [3]. The flavour diagonal action, denoted as Osterwalder Seiler fermion action, is given by with χ f a single-flavour fermion field. The renormalized valence masses µ OS,ren c,s /Z P can be matched to the corresponding renormalized quark masses via with Z P and Z S denoting the non-singlet pseudoscalar and scalar Wilson fermion quark bilinear renormalization constants. Then correlation functions using OS or unitary valence quarks are equivalent in the continuum. Moreover they still yield O(a) improved physical observables. The general idea to tune the heavy quark twisted mass parameters is to start with an educated guess in the unitary setup and to tune the OS charm and strange valence masses by imposing two suitably chosen physical renormalization conditions. The so determined parameters of the OS action, i.e. aµ OS s for the strange quark and aµ OS c for the charm quark, can then be translated to new heavy quark twisted mass parameters (aµ σ and aµ δ of eq. (4)) via eq. (10), in the unitary setup and, together with a slight retuning of κ crit , a new unitary simulation can be performed. With a convenient choice of the physical renormalization conditions, here C 1 and C 2 (see below), this parameter tuning procedure can be carried out on a non-large lattice (in the present case, 32 3 · 64) and at larger than physical up/down quark mass.
In this work, we follow the above described strategy. As physical conditions we choose where m Ds is the D s -meson mass and f Ds the D s -meson decay constant. The condition C 2 has a strong sensitivity to the charm quark mass while C 1 fixes the strange-to-charm mass ratio. They show only small residual light quark mass dependence arising from sea quark effects. We expect these conditions to be essentially free from finite-size effects due to the heavy D s -meson mass. This setup leads indeed to an only small error for the final parameter choices. Details on our measurements of meson masses and decay constants for twisted mass fermions are given in the Appendix A.
As a first step, we work on gauge ensembles produced with µ around three times larger than the physical up-down average quark mass and with educated guess values of µ σ , µ δ and m 0 . We choose the OS quark masses µ OS c and µ OS s such that condition C 1 is fulfilled. We then vary the OS quark masses, while maintaining condition C 1 , over a broad enough range such that also condition C 2 is satisfied within errors.
In a second step, we match the heavy charm and strange twisted mass of the unitary action (4) to the OS fermion quark mass parameters via eq. (10). The value of aµ σ is directly determined from aµ OS s and aµ OS c , while aµ δ is fixed by the ratio Z P /Z S . The latter can be estimated by adjusting aµ δ such that the kaon mass evaluated in the unitary formulation (m tm K ) and its counterpart computed with valence OS fermions (m OS K ) are equal. Although the kaon mass value can be unphysical due to having a too large value of µ and possible finite size effects, the matching condition actually relates only heavy quark action parameters. It fixes the relation of aµ δ to aµ OS s and aµ OS c , or equivalently the ratio Z P /Z S . In that way it is insensitive to both the finite lattice size and the actual value of µ up to O(a 2 ) artifacts. Since the matching steps described so far were implemented only on the valence quark mass parameters of the unitary and OS actions using gauge ensembles with so far different values of the sea quark mass parameters, one still needs to generate new gauge configurations at the so-determined values of aµ σ and aµ δ . Now on these new ensembles it can be re-checked whether the condition C 2 and the matching condition m tm K = m OS K , as well as the maximal twist condition eq. (7) in the light quark sectors, are fulfilled with sufficient accuracy. If this happens not to be the case, the procedure has to be iterated.
More concretely, we start with an initial guess for the heavy quark mass parameters given by aµ δ = 0.1162 and aµ σ = 0.1223, which we deduce from a number of tuning runs on a lattice of size 24 3 × 48 and 32 3 × 64 along the lines of ref. [52]. These parameters are realized for the ensemble Th1.200.32.k2, which is moreover very close to maximal twist. We then employ OS fermions in the valence sector and vary the values of µ OS s and µ OS c -while maintaining condition C 1 -such that condition C 2 is fulfilled. This is illustrated in Fig. 2 As explained above, the values in eq. (12) already determine aµ σ . To determine aµ δ , we first compute the kaon mass in the OS setup at aµ used in the unitary setup and aµ OS s from eq. (12). Having found the OS kaon mass, we go back to the ensemble Th1.200.32.k2 and tune in the unitary heavy quark valence sector µ δ such that we match the OS kaon mass. We then take the so found value of µ σ and µ δ for our simulations on the target large volume lattice. In this process a useful guidance is provided by assuming Z P /Z S = 0.8 known to be a typical value from our previous simulations. As we will discuss later, this assumption for Z P /Z S turns out to be rather close to the values we determine on the cB211.072.64 ensembles. Our final result for the action parameters in the heavy sector of maximally twisted mass fermions then read aµ σ = 0.12469 and aµ δ = 0.13151.
Due to the retuning of the heavy quark masses κ crit has to be re-tuned as well. To this end, several ensembles with volumes of 32 3 × 64 at light twisted mass values of aµ = 0.002 and aµ = 0.00125 were generated to determine the critical hopping parameter for the simulation at aµ = 0.00072 resulting in κ crit = 0.1394265. In this work, it turned out that we only needed one iteration of the above procedure using the Th1.200.32.k2 ensemble. After this first step, the tuning conditions for the heavy quark masses were checked again on the Th2.200.32.k2 ensemble (see table 2) and found to hold to a good accuracy within statistical errors. A similar finding holds also on our target ensemble cB211.072.64 ensembles. If we impose again an exact matching between m OS K and the unitary m tm K on the two cB211.072.64 ensembles we find the ratio of the pseudoscalar to the scalar renormalization constants to be Using this value of Z P /Z S the values of µ σ,δ of eq. (13) are close to the corresponding parameters at the physical point (the cB211.072.64 ensembles) that match our tuning conditions. Indeed the actually employed sea quark mass parameters correspond to a sea strange (charm) quark mass 6% lighter (4% heavier) than those derived a posteriori from imposing the same tuning and matching conditions on the physical point ensembles. It is also very nice to observe that by enforcing these conditions with very high precision one would obtain at the physical point with a ∼ 0.08 fm a kaon mass in isosymmetric QCD less than 1% smaller than its experimental value.

O(a 2 ) isospin-breaking lattice artifacts in the pion sector
An important aspect when working with twisted mass fermions at maximal twist is to keep the size of isospin violations small. This isospin breaking manifests itself by the fact that the neutral pion mass becomes lighter than the one of the charged pion. In leading order (LO) of chiral perturbation theory this effect is described by with the twisted mass angle given by ω = atan(µ /Z A m PCAC ) and c 2 a low energy constant characterizing the strength of O(a 2 )-effects of twisted mass fermions. As shown in Ref. [5,18], using a clover term the value of the low energy constant c 2 decreases. Indeed, employing a clover term, simulations at physical quark masses become possible as demonstrated in Ref. [5]. It turns out that c 2 < 0 for twisted mass fermions [53] leading to the the Sharpe-Singleton scenario [41]. In order to calculate the neutral pion mass one needs to compute disconnected twopoint functions that are notoriously noisy. To suppress the noise in the computation of the two-point functions we use a combination of exact deflation, projecting out the 200 lowest lying eigenvalues, and 6144 stochastic volume sources corresponding to an eight-distance hierarchical probing [54,55]. The disconnected correlator needed is given by where the ensemble and time average of the vacuum contribution is subtracted from the disconnected operator. Note that we used global volume noise sources to extract the disconnected contribution, however methods which do not subtract the vacuum expectation value explicitly could be more effective as pointed out in [5,6,56]. We have found that the disconnected contribution dominates the correlator for time distances t/a > 10, as can be seen in Fig. 3. However we include the connected contribution in the plateau average, leading to a neutral pion mass given by am π (0) = 0.044 (9) .
Note, that for the connected contribution small statistics of around 250 measurements is used, which results in a relatively large statistical error. The charged pion mass is straight forward to compute and we find for the charged pion mass am π = 0.05658 (6). This gives an isospin splitting in the pion mass of 22(16)% and the low energy constant c 2 of eq. (15) reads 4c 2 a 2 = −0.0013(8) assuming ω = π/2. Thus, introducing a clover term for N f = 2 + 1 + 1 twisted mass fermions suppresses isospin breaking effects effectively, i.e. by a factor of 6 compared to an N f = 2 + 1 + 1 ensembles with twisted mass fermions without a clover term and a pion mass of 260 MeV at a similar lattice spacing of a = 0.078(1) fm [53,57], where it was found that the mass splitting is given by (am π (0) ) 2 − (am π ) 2 = −0.0077(4) . The suppression of the pion isospin breaking effects, thanks to the use of the clover term, is the underlying reason why we can perform our simulations at the physical point with N f = 2 + 1 + 1 flavours of quarks.

Pseudoscalar meson sector
In order to check, whether we are indeed at (or close to) the targeted physical situation, we studied the charged pion, the kaon and the D-meson masses and decay constants. These observables are rather straightforward to compute with good accuracy. A detailed description of the calculation of these quantities with twisted mass fermions can be found in appendix A .

Light Meson sector
The first goal of this section is to determine the value of the lattice spacing within the pion sector. The extracted value will then be compared to the one from a similar investigation in the nucleon sector in section 5. In principle, the lattice spacing could be determined already from our cB211.072.64 target ensembles given in Table 2, having a twisted mass parameter of µ = 0.00072 and yielding a pion mass to decay constant ratio of m π /f π = 1.073 (3), which is rather close to the physical one. However, it is helpful to also use other ensembles, listed in Table 2, which are all tuned to maximal twist, namely Th1.350.24.k2, Th2.200.32.k2, Th2.150.32.k2 in addition to the cB211.072.64 ensembles. By employing chiral perturbation theory (χPT) to describe the quark mass dependence of the pion decay constant and pion mass, we obtain a robust result for the value of the lattice spacing. Since the ensembles that are not at the physical point have partly only a small volume, we include finite volume corrections from chiral perturbation theory to the χPT formulae used [58]. We depict in Fig. 4 the ratio m 2 π /f 2 π and the pion decay constant itself as function of the light bare twisted quark mass.
In Fig. 4 we also show the fits to NLO χPT [60,61,62], which for the ratio m 2 π /f 2 π read m 2 π f 2 π = 16π 2 ξ (1 + P ξ + 5ξ log (ξ )) and for the pion decay constant with the finite volume correction terms F F V E fπ , F F V E mπ [58]. Here ξ = 2B 0 µ /Z P [(4πf 0 ) 2 ] where B 0 and f 0 are low energy constants. From the fits, we determine the values of 2B 0 /Z P = 4.52 (6) and af 0 = 0.0502(3). The fitting constants P, R are related to the NLO low energy constants by We determine the finite volume correction terms by fixing the low energy constants using the results of Ref. [38]. For our target ensemble cB211.072.64 with m π L = 3.62 we find that the finite volume effects yield corrections of less than 0.5% for the pion The ratio m 2 π /f 2 π is plotted against aµ . Right: The pion decay constant af π is plotted against the light twisted mass math parameter aµ . The solid lines are fits to NLO chiral perturbation theory with the error as shaded band, see eq. (19) and eq. (20). The dotted lines are fits for which the chiral logs are neglected. The pion mass and decay constant are corrected for by finite volume correction terms [58] and [59] respectively. mass and less than 0.5% for the pion decay constant. By using the fit functions from χPT and fixing the ratio m 2 π,phys /f 2 π,phys ≡ 1.034 we find for the light twisted mass parameter aµ ,phys = 0.00067(1). We then use this value in eq. (20) to determine the lattice spacing. We get a fπ = 0.07986(15)(35) fm, (22) with the first error the statistical and the second the systematic by using the physical value of the pion decay constant, f π,phys = 130.41 (20) MeV [63]. We follow the procedure adopted in Ref. [39] for determining a systematic error by performing several different fits, adding or neglecting finite volume terms. Such fits employ e.g. the finite volume corrections of Ref. [59] using the calculated low energy constant c 2 of eq. (18) different orders in chiral perturbation theory and including or excluding the ensemble Th2.150.32.k2 due to larger finite size effects. The systematic error is then given by the deviations of these different fits from the central value given in eq. (22). Although we include ensembles like Th2.150.32.k2 or Th1.350.24.k2 which have large finite size effect of up to 8% in the pion decay constant, the systematic uncertainties are suppressed due to the fact that we are using ensembles close to physical quark masses which stabilize the fits. Thus this demonstrates the importance of working at physical quark masses. Moreover this is confirmed by an estimation of the lattice spacing which takes only the pion mass and decay constant from cB211.072.64 into account. Requiring a vanishing pion mass in the chiral limit, the lattice spacing and the physical twisted mass value can be fixed by assuming a linear dependence of µ on a 2 m 2 π and m 2 π /f 2 π . The so determined lattice spacing agrees with eq. (22) and reads a = 0.0801(2) fm.

Heavy meson sector
As discussed in section 3.2, the heavy sea quark parameters used in the simulation are tuned by employing the ensemble Th1.200.32.k2. With these parameters the kaon mass on the cB211.072.64 ensembles is smaller as compared to the OS kaon mass using the parameters of eq. (12). By employing the tuning condition of eq. (11) we therefore  (12)). By using aµ ,phys = 0.000674 the strange to light quark mass ratio reads µ OS s µ l = 0.01892(13) 0.00067(1) = 28.1 (5) .
The kaon and D-meson masses and the respective decay constants as well as the corresponding quantities for the D s -meson are all computed at three different values of µ OS s and µ OS c . We use a linear interpolation of m 2 K , m D and m Ds with respect to the heavy OS quark masses. Using the values for µ OS s and µ OS c of eq. (24) this allows us to determine the masses and decay constants for these mesons. In Fig. 5 we show the decay constants of the kaon and the D-meson and compare them with the results extracted from the N f = 2 clover ensembles [5]. We employ 244 measurements for the cB211.072.64 and 100 for the Th2.200.32.k1 ensemble. The ratios of the kaon and D-meson masses to decay constants for the cB211.072.64 ensembles are found to be (7) and where the former ratio has a central value slightly larger than the physical ratio m phys K /f phys K = 3.162(18) [64], while the latter agrees well within errors with the value m phys D /f phys D = 9.11 (22) [63]. These results indicate that discretization effects for our setup are small in the heavy quark sector. For a more rigorous check, a direct calculation at different values of the lattice spacing will be carried out. Table 3: The masses and the decay constants of the charged pseudoscalar mesons as well as the plaquette P and m PCAC are presented. am π = 0.05658 (6) am K = 0.2014(4) am D = 0.738 (3) m π /f π = 1.0731(30) m K /f K = 3.188 (7) m D /f D = 8.88 (11) am PCAC = 0.189(114)10 −4 P = 0.5543008(60)

Baryon sector
As another test, whether we are in the desired physical condition, we analyzed the nucleon mass which can also provide an independent determination of the lattice spacing, which can be compared to the one found in the meson sector. We measured the nucleon mass on the two cB211.072.64 ensembles by using interpolating fields containing the operator with C = γ 4 γ 2 the charge conjugation matrix. We then constructed the two point correlation function which provides the nucleon mass in the large time limit. We used 50 APE smearing steps with α AP E = 0.5 [65] in combination with 125 Gaussian smearing steps with α gauss = 0.2 [66,67] to enhance the overlap of the used point sources with the lowest state. We extracted the nucleon mass for t 0 by a plateau average over the effective mass aE ef f = log(C p (t + a)/C p (t)) shown in Fig. 6. The plateau average of the nucleon mass, given by am N = 0.3864(9) on the cB211.072.64 ensemble, is in agreement with a two-state fit with am N,2st = 0.3850 (12) as shown in the left panel of Fig. 6.

Determination of the lattice spacing
As an alternative way to determine the lattice spacing, one can use the nucleon mass. A direct way would be to use the physical ratio from which, by using the (lattice) pion mass determined above, the lattice spacing can be estimated directly by the value of the lattice nucleon mass. Indeed, with the the pion mass am π = 0.05658(6) the nucleon to pion mass ratio 0.3864(9)/0.05658(6) = 6.83(2) is close to its physical value of m phys N /m phys π = 0.9389/0.1348 = 6.965 where we take the average of neutron and proton mass [63] and the pion mass in the isospin symmetric limit [64]. However, as in the case of the meson sector, using more data points at heavier pion masses and χPT to describe their quark mass dependence, a more robust result can be obtained. More concretely, we have employed chiral perturbation theory at O(p 3 ) [68,69] for the nucleon mass dependence on the pion mass, i.e. Similar to [70] we use the nucleon masses of the N f = 2 + 1 + 1 ETMC ensembles without a clover term, determined in [71] to perform the chiral fit of eq. (28). In our analysis we neglect cutoff effects, which appear to be small and not visible within our statistical errors. The same holds true for finite volume effects, see [70]. We fixed f π = 0.1304(2) GeV and g A = 1.2723(23) [63] in eq. (28). The resulting fit to eq. (28) is shown in Fig. 6 (right) and allows to determine the lattice spacing as a m N (β = 1.778) = 0.08087(20)(37) fm .
The first error is statistical while the second error is the deviation between the estimate obtained from eq. (28) and taken the mass from the two-state fit. Note that the statistical error in the nucleon mass is comparable with pion sector due to three orders of magnitude larger number of inversions.

O(a 2 ) isospin splitting in the baryon sector
The finite twisted mass value can result into a mass splitting of hadrons which are symmetric under the isospin symmetry of the light flavor doublet. As pointed out in sec. 3.1 this indeed leads to a sizable effect in the neutral-charged pion mass splitting.
Here, we want to discuss the splitting in the baryon sector in case of the ∆-baryon employing the two cB211.072.64 ensembles. Note that for the used lattice size the lowest decay channel of the Delta baryon, which is a nucleon+pion state with correct parity, is heavier than the Delta baryon itself. Thus, for the simulations performed here, the Delta can be treated as a stable state.
We measured the ∆-baryon correlator by using the following interpolating fields Note that J µ ∆ + and J µ ∆ ++ is symmetric under u → d to J µ ∆ 0 and J µ ∆ − respectively. We neglect the potential mixing of ∆ with the spin-1/2 component which is suppressed [72]. Thus the correlators for the ∆ ++ is given by C ∆ = Tr[C]/3 with C ij = Tr[(1 + γ 4 )/2 J i ∆ ++ (t)J j ∆ ++ (0) ] and gives an average value of am ∆ = 0.5251(72) by using a plateau average over the effective mass. Now we define the splitting in the mass by where we average over the symmetric parts. In Fig. 7 we show the effective relative mass splitting given by δm ∆,ef f /m ∆ + . In addition we plot the relative effective mass m ef f (t) of the ∆ + particle subtracted from its plateau average to illustrate where the plateau of the ∆-baryon starts. We find that the relative splitting in the ∆ mass is δm ∆ /m ∆ + = 0.0098 (65) and hence close to zero within errors. This result is in agreement with [5] where it was found that the isospin splitting of the twisted mass action in the baryon section is suppressed.

Lattice spacing
The lattice spacing can be evaluated by matching lattice observables to their physical counterparts. This has been done, as described in sections 4.1 and 5.1, in the meson sector by employing the pion decay constant and in the baryonic sector using the nucleon mass, respectively. Differences in the values obtained for the lattice spacing as determined using different physical observables can shed light on cut-off effects. We discuss in this section an additional method to determine the lattice spacing, which is provided by the gradient flow scale setting parameters t 0 [73] and w 0 [74]. Following the procedure described in these articles and in particular as applied to the twisted mass setup [5], we extrapolate the gradient flow observables to the chiral limit using a fit ansatz linear in aµ , which corresponds to LO χpt [75]. The resulting curve is shown as Fig. 8. We follow a similar procedure for the extrapolation of In Table 4 we summarize the values of the lattice spacing as determined from the pion mass and decay constant, the nucleon mass and the gradient flow parameters t 0 and w 0 . As it can be noticed, there are small deviations of the lattice spacing between the meson and the baryons sector and in any case they are comparable to the one we have observed in the simulations with N f = 2 flavours of quarks. That indicates that cutoff effects do not increase for our N f = 2 + 1 + 1 flavour setup used here. We would like to stress, that we plan to carry out further simulations at different and, in particular, smaller values of the lattice spacing in future works.

Conclusions
The first successful simulation of maximally twisted mass fermions with N f = 2 + 1 + 1 quark flavours at the physical values of the pion, the kaon and the D-meson masses has been presented. By having a lattice spacing of a = 0.08029(41) fm, we find that the simulations are stable when performed with physical values of the quark mass parameters. In particular, we are able to carry out a demanding but smooth tuning procedure to maximal twist and to find the values of the light, strange and charm bare quark masses, which correspond to the physical ones for the first two quark generations. Table 4: We give the values of the lattice spacing determined by using different physical quantities as inputs, including in the errors the input systematic uncertainties. The final value of the lattice spacing is derived via a weighted average of a fπ and a m N where for the final error a 100% correlated data is assumed [76]. The residual systematic uncertainty on the lattice spacing, which stems from higher order cutoff effects, should be of relative size O(a 2 ) and looks numerically smaller than 2%.
[fm] quantities in lat. units a t 0 0.0811 (14) t 0 /a 2 | µ =0.00072 = 3.246 ( 7) a w 0 0.0823( 8) w 2 0 /a 2 | µ =0.00072 = 4.512 (16) a fπ 0.07986 (38) af π | µ =0.00072 = 0.05272 (10) a m N 0.08087 (44) m N /m π | µ =0.00072 = 6.829 (19) average 0.08029 (41) In our setup, which employs a clover term, the cutoff effects appear to be small. Several observations corroborate this conclusion: as already mentioned above, the simulations themselves are very stable; when fixing the quark mass parameters through the selected physical observables, other physical quantities, as collected in Table 3 come out to be consistent with their physical counterparts; the O(a 2 ) effects originating from the isospin breaking of twisted mass fermions are small and significantly reduced compared to our earlier simulations with N f = 2 + 1 + 1 flavours at non-physical pion masses; deviations of the lattice spacing from the meson sector, the baryon sector and gradient flow observables, as listed in Table 4, are small and of the same size as in our former N f = 2 flavour simulations.
This work focuses on the tuning procedure both to maximal twist and to the physical values of the quark masses. We include the pseudoscalar meson masses and decay constants as well as the nucleon and ∆ masses, in order to demonstrate that we indeed reach the targeted physical setup. We are planning to compute many more quantities in the future connected to hadron structure, scattering phenomena, electroweak observables and heavy quark decay amplitudes. In addition, we have already performed the tuning for a second, finer lattice spacing and we are in the process of generating configurations. The combination of results for various physical quantities from the present lattice spacing of a ≈ 0.08 fm, from the ongoing finer lattice spacing and from an already existing lattice spacing of a ≈ 0.1 fm, which is however not exactly at the physical point, will allow us to explicitly check the size of cut-off effects and eventually take the continuum limit.
We thus conclude that we have given a successful demonstration that simulations of maximally twisted mass fermions with N f = 2 + 1 + 1 quark flavours can be carried out with all quarks of the first two generations tuned to their physical values. This clearly opens the path for the ETM collaboration to perform simulations towards the continuum limit with a rich research programme being relevant for phenomenology and ongoing and planned experiments.

A Mesonic Correlators
In general, the charged 2-point pseudoscalar correlators can be defined by C q,q PS (t) = P ± q,q (t) P ± q,q (0) † using the interpolating field P ± q,q (t) = xχ q (x, t) iγ 5 τ ± χ q (x, t) , with the quark flavors q, q ∈ { , s, c}. For sufficiently large times the charged pseudoscalar correlator is dominated by the lowest energy, such that 2m P S e −m P S t + e −m P S (T −t) (36) and the mass m P S and matrix element G P S = | 1 P S |P ± q,q |0 | can be extracted in a standard way via plateau averages for large time extent. In case of maximal twist the matrix element G P S is directly connected to the pseudoscalar decay constant by [1,2] f π = (µ q + µ q )G P S sinh(m P S )m P S .
Due to the flavor mixing in case of the mass non-degenerate twisted mass operator we adopt a non-unitary setup [3] for the heavy quark doublet, namely the Osterwalder-Seiler fermion regularization [51]. As shown in [3], this mixed action introduces effects which are only of order O(a 2 ) and are hence suppressed for small µ and fine lattice spacings. The OS fermions correspond to the twisted mass discretization in single flavor space, where µ = ±µ q . The sign of µ is always chosen such that the two valence quarks in the interpolating fields eq. (35) have opposite signs.

B Autocorrelation
The autocorrelation of the Hybrid Monte Carlo algorithm slows down critically for very fine lattice spacings with a < 0.05 fm. This can be seen in the freezing of the  topological charge [77]. For our lattice with a ∼ 0.08 fm we found that the topological charge can fluctuate between the different sectors leading to small autocorrelation times of τ int (Q) = 13(5) [MDU]. As pointed out in [78] the energy density at finite flow times develops larger autocorrelation times in the regime with a 0.05 fm. Although we have relative small statistics we calculated integrated autocorrelation time for the plaquette and the gradient flow observables t 0 /a 2 as shown in figure 9 by using the Γ-method [79]. We found a pion mass dependence given by with b = 2.2(5) case of the plaquette while b = 2.0(7) in case of the gradient flow observable t 0 . A possible explanation for the quark mass dependence of the autocorrelation time τ int is a phase transition in case of finite twisted mass term for vanishing neutral pion masses. Although the isospin splitting is suppressed in our case, observables like the gradient flow observables shows an increase with inverse of the squared pion mass. This behavior is also seen in the PCAC mass, where moderate integrated autocorrelation times were found which can be clearly seen in the Monte Carlo history (see right panel of fig. 1). However in other quantities like the pseudoscalar mass, the pseudoscalar decay constant or nucleon observables τ int is very small and a quark mass dependence can be not observed.