Simulation of an ensemble of $N_f=2+1+1$ twisted mass clover-improved fermions at physical quark masses

We present a general strategy aimed at generating $N_f=2+1+1$ configurations with quarks at their physical mass using maximally twisted mass fermions to ensure automatic $O(a)$ improvement, in the presence of a clover term tuned to reduce the charged to neutral pion mass difference. The target system, for the moment, is a lattice of size $64^3 \times 128$ with a lattice spacing $a\sim 0.08$ fm. We show preliminary results on the pion and kaon mass and decay constants.


Introduction
During the last decade, lattice QCD simulations with light quarks having physical masses on large enough volumes have become feasible owing to improvements in lattice actions and numerical algorithms as well as due to steadily increasing computational power. In this work, we illustrate a strategy aimed at generating N f = 2 + 1 + 1 configurations with quarks at the physical mass using maximally twisted fermions to ensure automatic O(a) improvement. This bonus does not come without a price, as the twisted regularization breaks isospin leading to an undesired (negative) neutral pion to charged mass splitting. It has been shown in [1] that a suitably tuned clover term can ease the problem and allow for stable simulations at affordable lattice spacings and volumes (sec. 2). The tuning procedure is quite complicated as a number of parameters (bare quark masses, critical mass, and the clover term coefficient, c S W ) have to be simultaneously set at appropriate values. In sec. 3 we discuss a viable strategy that allows with relative small effort and good signal-to-noise ratio to desired tuning in a way that weakly depends on the light quark sector. Going to large lattice volumes, as needed to simulate pions with close-to-physical mass, leads to an increased density of very small Dirac operator eigenvalues, which is expected to stabilise the MD behaviour of HMC-like algorithms but also requires the use of very well preconditioned and efficient linear solvers. We will discuss the practical aspects of our Monte Carlo simulation strategy and related issues in sec. 4. Finally in sec. 5, we will present results on the pion and kaon masses and decay constants using the N f = 2 + 1 + 1 ensemble.

Setup
The action used in our simulation is given by with S g the Iwasaki improved gauge action [2], β the bare coupling constant and χ ℓ and χ h pseudofermion fields. The light quark-doublet is made of the up-and down-quarks, described in the isospin symmetric limit by a degenerated twisted mass operator with twisted mass parameter µ ℓ . The heavy quark doublet is made by the strange and charm quarks and are described by the mass non-degenerated twisted mass operator given by with bare mass parameter µ σ and µ δ . The Wilson Dirac operator D W (κ, c sw ) depends on the bare mass parameter κ = 1 2(4+m) . A Sheikholeslami-Wohlert-term [3] is added. For degenerated masses one gets that det{D † (µ ℓ )D(µ ℓ )} = det{D 1+1 (µ ℓ , 0)}.
In the twisted mass discretization of lattice QCD all the correlators that are expectation values of parity even multilocal operator products are by symmetry free of linear lattice artefacts if the Partially Conserved Axial Current (PCAC) mass vanishes, which is equivalent to tuning the bare mass parameter κ(m) to its critical value [4,5]. In this case, the renormalized light quark mass is at maximal twist and related to m ℓ = µ ℓ /Z P with Z P the pseudoscalar renormalization constant. This fine tuning of κ(m) → κ crit (m) will be discussed in the next section.
The twisted mass Wilson-Dirac operator corresponding to the light quark sector is protected from zero eigenmodes for a finite twisted mass term µ ℓ > 0. In fact D † (µ ℓ )D(µ ℓ ) = D † W D W + µ 2 ℓ with min{λ} ≥ µ ℓ . From a numerical point of view this holds the promise to be able to reach small quark masses with µ ℓ > 0 including the physical value.
However, the lattice regularization breaks some continuum symmetries that are recovered after extrapolating to zero lattice spacing. Besides the breaking of the chiral symmetry due to the Wilson term, a finite twisted mass term breaks isospin symmetry by introducing a mass-splitting in the pion triplet.
The neutral to charge pion mass splitting can be estimated in chiral perturbation theory. To first order one finds [6] (m 2 where c 2 depends on low energy constants and it is found to be negative for twisted mass fermions [7]. This means that the neutral pion mass may vanish for finite light quark masses triggering a phase transition into a non-physical phase induces by lattice artefacts.
In Fig. 1 we show the isospin mass splitting between the charged and neutral pion. For N f = 2 + 1 + 1 twisted mass fermion ensembles with c S W = 0 the mass splitting between the neutral and charged pions is sizable for the three different lattice spacings examined. It follows that for lattice spacing a 0.05 fm the neutral pion mass would vanish for physical light quark masses making simulation at the physical point impossible. However, it was shown that by adding a clover term the isospin breaking effects can be suppressed [1]. This was demonstrated in the case of the N f = 2 twisted mass fermion ensembles with c S W = 1.57551 and a lattice spacing of a = 0.0938 fm. In particular, for a charged pion mass of m π = 130 MeV the pion mass difference between the neutral and charge pions is smaller than |m 2 π 0 − m 2 π ± | 6 MeV and simulation of twisted mass fermions at physical light quarks are possible for a < 0.1 fm.  Figure 1. Left: The squared neutral pion mass is plotted against the squared charged pion mass in units of w 0 for N f = 2 + 1 + 1 ETMC ensembles with no clover term i.e. taking c S W = 0. We use a generalized fit excluding ensembles with m π > 400 MeV and show results for three lattice spacings and extrapolating to a = 0. The line corresponds to the physical value of the charged pion mass. We used the definition of [8] for w 0 . Right: The neutral pion mass for the N f = 2 ensemble with a clover term.

Tuning towards physical quark masses
We start by discussing the tuning procedure towards physical quark mass, using the N f = 2 + 1 + 1 twisted mass fermion discretization at a target lattice spacing of a ∼ 0.8 fm. The inclusion in the action of a clover term suppresses isospin violation effects. Using an estimate from 1-loop tadpole boosted perturbation theory given by c S W 1 + 0.113g 2 plaq [9] with g 2 plaq = g 2 0 /P and P the plaquette expectation value. This yields c S W = 1.69. The tuning procedure consists of two major parts: i) tuning towards the critical mass and ii) tuning the heavy quarks towards their physical values. The two steps are strongly correlated. To tune to critical mass we use the vanishing of the PCAC mass, which depends on the light valence quarks and to tune the strange and charm quark masses we use dimensionless rations depending only on the strange and charm valence quarks.
with Z A the axial renormalization constant roughly given by ∼ 0.8. For the tuning of the heavy quark masses, we introduce an additional step by setting the mass parameters µ s and µ c of the Osterwalder Seiler fermions such that The final step is to match the masses of the Kaons involving an Osterwalder Seiler and a nondegenerated twisted mass valence strange quark, respectively. In this way the matching of valence and sea quark mass parameters is independent of finite size effects and with a negligibly small dependence on the light quark mass entering only through sea quark effects. The approach was as follows • We start with simulation using lattices of spatial extent N s = 24 at four different mass values κ and m π ∼ 310 MeV. The parameters β and µ σ , µ δ are initially estimated by a set of preliminary studies aimed at investigating their interdependencies [11]. Using a linear fit for the PCAC mass as a function of the bare light quark mass parameter we estimate κ crit as shown in Fig. 2) and recheck the value of c sw and the lattice spacing a by using Wilson flow observables.
• With the estimated critical mass value we simulate a larger lattice size with spatial extent N s = 32 and we check and re-tune the heavy quark masses by calculating the properties of the kaon and D-mesons.
• Based on the re-tuned heavy quark parameters µ σ and µ δ we run three simulations again with spatial lattice size N s = 32 at different κ values and a pion mass of m π ∼ 230 MeV. Again we estimate the critical mass by using a linear fit to the PCAC mass for the generated ensembles as shown in Fig. 2. We then cross-check the estimated value by simulating an ensemble with a pion mass m π ∼ 180MeV.
• Finally we simulate our target lattice size with N s = 64 3 and temporal extent N T = 128 using pion masses close to the physical value. The PCAC masses of the target ensemble are shown Fig. 2 and final tuned parameters are found in Table 1. On the first 1200 MDUs we find that the average PCAC mass satisfies |m PCAC | < 0.00005, which fulfills the bound given in Eq. (4) that ensures sufficient suppression of the linear lattice artefacts.
We use the Hybrid Monte Carlo (HMC) algorithms [12] with Hasenbusch mass-preconditioning [13] for the light quark sector and the rational HMC [14] for the heavy quark sector. The software package tmLQCD [15] with the multi-grid solver DDalphaAMG [16] is employed. For the inversions of the light and heavy quark sector a three-level DD-αAMG solver [17,18], and even-odd preconditioned conjugate gradient (CG) solver and a multi-mass shifted CG solver are used (see Ref. [19] for more details). The gauge fields are integrated by using a 6-level nested minimal norm scheme of second order. We show the distribution of the maximal force terms for the light and heavy quark sectors in Fig. 3. The rational approximation is done for the interval [0.000014 1] by involving 10 terms, while a correction term is included in the acceptance step. Using 12 integration steps on the outer level and 384 on the inner level, the integration with length τ = 1 yields a stable simulation, where the maximal energy violation is given by max(δH) = 5.67 with an acceptance ratio of 77%. In the case of the plaquette it can be seen that the use of the Hasenbusch mass-preconditioning yields a quite stable run. As shown in Fig. 4 the history of the plaquette shows some larger oscillations for MDU< 800. This provides an estimate for the integrated autocorrelation time τ int . We find τ int = 31 (13). Similar oscillations are found in the history of the PCAC mass, also shown in Fig. 2. They are seen to be anticorrelated with the plaquette oscillations. Nevertheless, we find that other observables, such as the topological charge, are not affected, and shows frequent oszillations. As shown in Fig. 4 the algorithm can efficiently sample the topological charge with an autocorrelation time of τ int = 8(3).

Observables
In this section, we show preliminary results extracted from two-point functions using our target ensemble of tab. 1. See Ref. [20] for a detailed discussion on how these measurements are performed. In the left panel of Fig. 5 we show the pion mass dependence of the pion decay constant f π in comparison with results obtained using N f = 2 + 1 + 1 twisted mass fermions with c S W = 0 and smallest pion mass of around m π ∼ 210 MeV in units of the gradient flow observables t 2 0 /w 2 0 a 2 . The results are fitted to a quadratic function in m 2 π to guide the eye. Note that finite size effects were not corrected. The measured pion decay constant of our target ensemble agrees with the continuum value, which was determined using continuum extrapolations of the gradient flow observables determined in [8] and a pion decay constant of 130.2 MeV.
In the right panel of Fig. 5 we show the nucleon mass m N in units of t 2 0 /w 2 0 a 2 as a function of the charged pion mass squared m 2 π . We fit the results using N f = 2 + 1 + 1 twisted mass fermion ensembles without a clover term with the help of leading order chiral perturbation theory [21,22]. The agreement between the chirally extrapolated value with the value obtained using our current physical N f = 2 + 1 + 1 ensemble indicates small cut-off effects for the nucleon mass in agreement with previous observations [23].
In the heavy quark sector, we use our tuned values of the strange and charm quarks to evaluate the decay constant of the kaon and D-mesons. In Fig. 6 we show the ratio of the pseudoscalar mass over the pseudoscalar decay constant as a function of the pion mass and compare it to the value obtained using the physical N f = 2 ensemble with a clover term [1]. The ratio involving the D-meson agrees with the expectation while the ratio by the Kaon is above the physical value. A possible reason for this behavior could be the tuning conditions which primarily affects the charm quark mass possibly yielding large cut-off effects.  Figure 6. Left: The ratio of the decay constant of the D-meson over its mass for the N f = 2 + 1 + 1 ensemble (red star) and for the N f = 2 twisted mass clover ensemble with a pion mass of m π = 130 MeV (magenta star). The green symbols correspond to the continuum points, using continuum extrapolated gradient flow observables t 0 and w 0 for n f = 2 + 1 in [8] (square) and for the N f = 2 in [24] (diamond). Right: The ratio of the decay constant of the kaon over its mass as a function of the pion mass in units of t 2 0 /w 2 0 with the same symbols as those used for the D-meson.

Conclusion
In this work we illustrate how it is posible to set up stable simulations for N f = 2 + 1 + 1 twisted mass fermions with a clover term to maximal twist and approximately physical values of the light, strange and charm quark masses. Introducing a clover term help to supress isospin breaking effects in the pion triplet and makes the simulation within the twisted mass fermion approach at the physical point feasible. Our tuning conditions for fixing the heavy quark masses and the light bare quark mass to maximal twisted are presented. Moreover, we show that, by using Hasenbusch-mass preconditioning combined with the RHMC, simulations are stable and that the force terms are under control. Using the N f = 2 + 1 + 1 ensemble we show results on the pion decay constant, the nucleon mass and the ratios of the pseudoscalar mass over the decay constant [20]. We find a nucleon to pion mass ratio m N /m π = 6.87(7) and kaon and D-meson mass ratio m K /m π = 3.58 (2) and m D /m π = 13.2(9). Additional we listed pseudoscalar meson masses and their ratio with their decay constant in tab. 2. The lattice spacing as determined from the nucleon mass turns out to be a ∼ 0.082(2) fm.