Simulation of an ensemble of $N_f=2+1+1$ twisted mass clover-improved fermions at physical quark masses
Jacob Finkenrath Constantia Alexandrou Simone Bacchio Panagiotis Charalambous Petros Dimopoulos Roberto Frezzotti Karl Jansen Bartosz Kostrzewa Giancarlo Rossi Carsten Urbach
hep-lat




DESY 17-165






Computation-based Science and Technology Research Center, The Cyprus Institute, Nicosia,  20, Constantinou Kavafi Str., Nicosia 2121, Cyprus
Department of Physics, University of Cyprus, P.O. Box 20537,  Nicosia 1678, Cyprus
Bergische Universität Wuppertal, Wuppertal, Germany
Dipartimento di Fisica, Università di Roma “Tor Vergata" and INFN,
Sezione di Roma 2, Via della Ricerca Scientifica - 00133 Rome, Italy
Centro Fermi - Museo Storico della Fisica e Centro Studi e Ricerche
Enrico Fermi, Compendio del Viminale, Piazza del Viminiale 1, I-00184, Rome, Italy
NIC, DESY, Zeuthen, Germany
HISKP (Theory), Rheinische Friedrich-Wilhelms-Universität Bonn, Bonn, Germany




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 × 128 with a lattice spacing a∼ 0.08 fm.
We show preliminary results on the pion and kaon mass and decay constants. 



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

    
Jacob Finkenrath1Speaker, j.finkenrath@cyi.ac.cy Constantia Alexandrou1,2 Simone  Bacchio2,3Panagiotis Charalambous1,2

Petros Dimopoulos4,5Roberto Frezzotti4Karl Jansen5Bartosz Kostrzewa6

Giancarlo Rossi4,5Carsten Urbach6

    December 30, 2023
=====================================================================================================================================================================================================================================





§ 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 <cit.> 
that a suitably tuned clover term can ease the problem and allow for stable simulations at 
affordable lattice spacings and volumes (sec. <ref>). 
The tuning procedure is quite 
complicated as a number of parameters (bare quark masses, critical mass, and the clover term coefficient,
c_SW) have to be simultaneously set at appropriate values.
In sec. <ref> 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. <ref>. Finally in sec. <ref>, 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

    S = β· S_g  - χ_ℓ^†{ D^†(μ_ℓ)D(μ_ℓ)}^-1χ_ℓ - χ_h^†{D_1+1^†(μ_σ,μ_δ) D_1+1(μ_σ,μ_δ) }^-1/2χ_h

with S_g the Iwasaki improved gauge action <cit.>, β 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
  
    D_1+1(μ_σ,μ_δ) =  D_W(κ,c_sw) ⊗1 + i μ_σ (γ_5 ⊗τ_3) - μ_δ (1 ⊗τ_1) = [ D_W + i γ_5 μ_σ            -μ_δ;            -μ_δ D_W - i γ_5 μ_σ ]
  
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 <cit.>
is added. For degenerated masses one gets that
{ D^†(μ_ℓ)D(μ_ℓ)} = { 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 <cit.>. 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 <cit.>

    (m_π^0^2 - m_π^±^2) ≈  4 · c_2 · a^2

where c_2 depends on low energy constants and it is found to
be negative for twisted mass fermions <cit.>. 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. <ref> we show the isospin mass splitting between the charged
and neutral pion. For N_f = 2+1+1 twisted mass fermion ensembles with c_SW = 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 <cit.>.
This was demonstrated in the case of  the N_f = 2 twisted mass fermion  ensembles with c_SW = 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_π^0^2 - m_π^±^2 |)≲ 6  MeV
and simulation of twisted mass fermions at physical light quarks are 
possible for a<0.1  fm.
   

  
   
   





§ 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_SW≅ 1+0.113 g^2_plaq <cit.> with 
g^2_plaq = g_0^2/P and P the plaquette expectation value.
This yields c_SW = 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.






For the tuning towards the critical mass we impose <cit.>

    Z_A m_PCAC(κ)/μ_ℓ≲ 0.1

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

    A_1 = μ_c/μ_s≡ 11.8     and     B_2 =  m_D_s/f_D_s≡ 7.9  .

The final step is to match the masses of the Kaons involving an Osterwalder Seiler
and a non-degenerated 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 <cit.>. 
	 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. <ref>) 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. <ref>.
        We then cross-check the estimated value by simulating an ensemble with a pion mass m_π∼ 180 MeV.
  
  * 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. <ref>
and final tuned parameters are found in Table <ref>.
On the first 1200 MDUs we find that the average
PCAC mass satisfies |m_PCAC| <  0.00005, which 
fulfills the bound given in Eq. (<ref>) that ensures sufficient suppression of the linear lattice artefacts.
  



§ STABLE ALGORITHMS


  
We use the Hybrid Monte Carlo (HMC) algorithms <cit.>
with Hasenbusch mass-preconditioning <cit.> for the light quark sector 
and the rational HMC <cit.> for the heavy quark sector. 
The software package tmLQCD <cit.> with the multi-grid solver DDalphaAMG <cit.> is employed.
For the inversions of the light and heavy quark sector 
a three-level DD-αAMG solver <cit.>, and even-odd preconditioned conjugate 
gradient (CG) solver and a multi-mass shifted CG solver are used (see Ref. <cit.> 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. <ref>.
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. <ref> 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. <ref>. 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. <ref> 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. <ref>. See Ref.<cit.> for a detailed discussion on how these measurements are performed.
In the left panel of Fig. <ref> 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_SW=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<cit.> and a
pion decay constant of 130.2  MeV.

In the right panel of Fig. <ref> 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 <cit.>. 
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 <cit.>.

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. <ref>  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 <cit.>.
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.





§ 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 <cit.>.
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. <ref>.
The lattice spacing as determined from the nucleon mass  turns out to be a∼ 0.082(2)  fm.





 §.§ Acknowledgments

We thank the members of the ETM collaboration for helpful and inspiring discussions. Additional thanks goes
to U. Wenger for providing the gradient flow data for the N_f=2+1+1 ensembles with c_SW=0.
J. F. thanks  PRACE Fourth Implementation
Phase (PRACE-4IP) program of the European Commission under grant agreement No 653838 for financial support.
This project has received funding from the  Horizon 2020 research and innovation program of the European Commission
under the Marie Sklodowska-Curie grant agreement No 642069.
Part of this work was supported by the Sino-German CRC110.
The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. for funding the project pr74yo by 
providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre.




