Source-linked AI summary
Dynamical Twisted Mass Fermions with Light Quarks: Simulation and Analysis Details
Ph. Boucaud, P. Dimopoulos, F. Farchioni, R. Frezzotti, V. Gimenez, G. Herdoiza, K. Jansen, V. Lubicz, C. Michael, G. Münster, D. Palao, G. C. Rossi, L. Scorzato, A. Shindler, S. Simula, T. Sudmann, C. Urbach, U. Wenger, ETM Collaboration
TL;DR
The paper addresses the need for detailed simulation and analysis procedures behind two-flavour twisted mass lattice QCD at maximal twist. It develops tuning, correlator, error, algorithm, scale, disconnected-contribution, and chiral-analysis methods, reporting precise mesonic quantities and acceptable control of disconnected contributions. The analysis remains limited in assessing cutoff-effect coefficients because data at only one lattice spacing are used.
Problem
The preceding study omitted details needed to assess simulations, maximal-twist tuning, correlator computations, error estimates, autocorrelations, scale setting, disconnected contributions, and chiral descriptions.
Method
The paper supplies theoretical and practical maximal-twist tuning criteria and analyzes correlators, disconnected contributions, autocorrelations, scale setting, and chiral perturbation theory descriptions.
Results
The methods yield precise mesonic quantities and allow quark-disconnected contributions to be computed to acceptable accuracy on unquenched gauge configurations.
Takeaways & Limitations
The study provides detailed procedures for conducting and analyzing Nf = 2 mass-degenerate Wilson twisted mass simulations at maximal twist.
Takeaways & Limitations
With data at only one lattice spacing, the paper cannot evaluate scaling behaviour or determine the coefficients multiplying a2 cutoff terms for the observables studied.
Abstract
from arXiv · showhide
In a recent paper [hep-lat/0701012] we presented precise lattice QCD results of our European Twisted Mass Collaboration (ETMC). They were obtained by employing two mass-degenerate flavours of twisted mass fermions at maximal twist. In the present paper we give details on our simulations and the computation of physical observables. In particular, we discuss the problem of tuning to maximal twist, the techniques we have used to compute correlators and error estimates. In addition, we provide more information on the algorithm used, the autocorrelation times and scale determination, the evaluation of disconnected contributions and the description of our data by means of chiral perturbation theory formulae.
1 Twisted mass fermions
The section introduces two-flavour Wilson twisted mass lattice QCD and explains how tuning the untwisted mass enables controlled discretization effects. It also details practical maximal-twist tuning and its numerical precision requirements.
- 1 Twisted mass fermions: The simulations use two mass-degenerate Wilson twisted mass quarks with the tree-level Symanzik improved gauge action.The Wilson parameter is set to r = 1, while the gauge action uses b1 = −1/12 and b0 = 1 − 8b1.
- 1 Twisted mass fermions: Maximal twist is defined by tuning the bare untwisted mass m0 to its critical value mcrit, equivalently through the hopping parameter κ.The relation is κ = 1/(8 + 2am0).
- 1.2 Maximal twist and residual O(a2) artifacts: Tuning to maximal twist can provide automatic O(a) improvement, leaving parity-even expectation values with O(a2) discretization errors in the continuum limit.The practical criterion uses a vanishing untwisted PCAC mass, imposed at the lowest simulated twisted mass when appropriate.
- 1.2 Maximal twist and residual O(a2) artifacts: Residual cutoff effects can grow as the pion mass decreases, so the strategy is to tune ξπ(µq) to zero or reduce it to O(am^2_QCD) by adjusting κcrit.Setting mPCAC = 0 at µq,min sufficiently suppresses ξπ when µq < ΛQCD.
- 1.3 Numerical precision for tuning to maximal twist: For lattice spacings near 0.1 fm, |ε/µq| ≲ 0.1 is considered acceptable for determining the critical mass, ideally decreasing with a.This condition need only hold at µq,min, which should remain roughly fixed in physical units as a decreases.
- 1.3 Numerical precision for tuning to maximal twist: Because the analysis uses data at only one lattice spacing, it cannot determine the coefficients of residual a2 cutoff effects or directly evaluate scaling behaviour.Preliminary collaboration results are reported as indicating small residual effects consistent with O(a) improvement.
2 Computations in the charged meson sector
Charged-meson correlators are computed in the twisted basis using stochastic time-slice sources, linked-source extensions, and simultaneous correlated fits to multiple correlators. These methods reduce inversions while enabling extraction of energies and decay-related quantities.
- Operator interpretation: At maximal twist, twisted-basis operators require translation to physical-basis fields because parity and isospin are not exact lattice symmetries.Operators can interpolate states with opposite parity or different isospin, producing O(a) contaminations in some correlators.
- Charged correlators: The charged pion correlator contains only connected diagrams and can be computed from one-flavour propagators using a propagator identity.The identity relates the d-quark propagator to the u-quark propagator, so a common source suffices for the charged channel.
- Correlator fits: Simultaneous fits to N × N correlator matrices with M states determine energies and couplings, from which quantities such as afπ and amPCAC are evaluated.Correlated fits use the covariance structure to assess whether the fitted χ2 is acceptable.
- Stochastic sources: Stochastic time-slice sources exploit gauge configurations efficiently while maintaining a manageable noise-to-signal ratio, with one noise sample per configuration sufficient in the described setup.The one-end trick is formulated for zero-momentum pseudoscalar sources, while linked sources extend it to general Dirac structures.
- Stochastic sources: Linked fuzzed sources extend the one-end trick to nonzero momentum and spatially nonlocal mesonic operators, at the cost of additional inversions.The method combines propagators from fuzzed and unfuzzed linked sources.
3 Computations in the neutral meson sector
Neutral-meson calculations must include disconnected contributions and account for parity and isospin violations induced by twisted-mass lattice effects. The analysis uses connected and disconnected correlators built from stochastic sources and operator matrices.
- Lattice-symmetry effects: O(a2) cutoff effects can violate parity and isospin, causing neutral and charged mesons to have different masses and requiring disconnected diagrams for neutral isovector mesons.These violations also allow correlators to receive contributions from states with different formal parity or isospin.
- Neutral correlators: Neutral correlators combine connected contributions with disconnected loop contributions, unlike the charged pion correlator.The disconnected contribution is evaluated from u-quark loops at both source and sink, with variance-reduction methods described separately.
- Operator interpretation: Neutral-meson operators are associated with continuum states only in the a → 0 limit, ignoring O(a) contaminations from states of different parity and isospin.The operator correspondence includes both isotriplet and isosinglet channels.
- Computational setup: Neutral calculations use local and fuzzed sources to form 6 × 6 or 2 × 2 correlator matrices, with linked stochastic sources for u- and d-quark propagators.The source construction extends the charged-meson setup to the additional neutral-channel contractions.
- Momentum dependence: The study distinguishes vector-meson correlators with momentum parallel and perpendicular to the spin, while omitting other additional correlators at anisotropic momentum.The comparison can probe mixing of rho mesons with their ππ decay products.
4 Simulation algorithm and error analysis
The simulations use mass-preconditioned HMC with multiple time scales, while statistical uncertainties are assessed through autocorrelation-aware methods. The lightest ensemble exhibits especially long PCAC-mass autocorrelations, motivating careful blocking and cross-checks.
- Simulation ensembles: All B1–B5 ensembles use β = 3.9, κ = 0.160856, and 24^3 × 48 lattices, with approximately 5000 equilibrated trajectories per twisted-mass value.The lightest mass ensemble was extended from about 5000 to approximately 10000 trajectories by combining replicas.
- Simulation algorithm: The gauge configurations are generated with mass-preconditioned HMC and multiple time-scale integration using SW, 2MN, or 2MNp schemes.Trajectories have length 1/2 and use two pseudo-fermion fields.
- Simulation algorithm: At the lightest quark mass, 5000 trajectories cost about 17 rack days on BlueGene/L, with approximately 18% code efficiency and about 115 Tflop per trajectory.These figures characterize the computational cost of the production run.
- Error analysis: Autocorrelation times are estimated with both data blocking and the Γ-method, while jackknife, bootstrap, and covariance-aware fits account for statistical dependence.The Γ-method also determines autocorrelation times for non-primary fermionic observables.
- Autocorrelation times: The PCAC-mass estimator has substantially longer autocorrelation times than pseudoscalar mass and decay-constant observables, attributed to the twisted-mass phase structure.For the lightest ensemble, the integrated autocorrelation is 32(9) trajectories.
5 The scale from the static potential
The paper determines the hadronic scale r0/a from the static potential using variationally optimized Wilson loops, correlated fits, and interpolation in the quark separation. It assesses excited-state, interpolation, mass-dependence, and autocorrelation effects to quantify the reliability of the result.
- Potential extraction: The static potential is extracted from Wilson loops using improved temporal links, spatial smearing, and a variational method to enhance ground-state overlap.Five smearing levels produce a 5 × 5 correlation matrix, whose largest-eigenvalue projection targets the ground state.
- Scale determination: The scale r0/a is obtained by fitting the interpolated potential, with final uncertainties estimated by jackknife and bootstrap procedures using binning factor 4.Table 10 reports measurement counts, χ2 per degree of freedom, and r0/a results.
- Interpolation and fitting: The interpolation uses a two-step procedure: determine V(r) separately at each distance, then fit a potential ansatz while including spatial and temporal cross-correlations.Potential uncertainties are estimated with a non-parametric bootstrap.
- Potential extraction: The analysis selects effective-mass plateaus by balancing excited-state contamination at small t against noise at large t, using correlated χ2 tests.The selected fit ranges are illustrated for r/a = 4 across ensembles B1–B5.
- Systematic checks: Residual autocorrelations may cause the quoted r0/a errors to be somewhat underestimated because fits become unreliable beyond bin size 4 before binning errors stabilize.The interpolation ambiguity is typically covered within 1–2 standard deviations when using r/a = 4–7.
- Systematic checks: The mass dependence of r0/a is well described by ansatz (I), while possible spontaneous-chiral-symmetry-breaking effects near 0.5 fm are negligible within statistical errors.Ansatz (II) cannot be completely ruled out, but ansatz (I) has better χ2/d.o.f. and is supported by fixed-distance potential fits.
6 Some selected results
The paper details charged and neutral pseudoscalar analyses, including correlator fitting, disconnected contributions, decay constants, renormalization, and maximal-twist checks. Results include compatible operator determinations, controlled excited-state and two-pion effects, and a substantial neutral–charged mass difference alongside consistent decay constants.
- Charged pseudoscalar channel: The charged pseudoscalar analysis extracts masses and decay constants from correlators, with errors estimated using the procedures described in sections 2.1 and 4.1.Results are reported in Table 12.
- Charged pseudoscalar channel: Three interpolating operators give compatible effective masses from t/a ≈10, supporting ground-state dominance for t/a > 9.The operators are local-local, local-fuzzed, and fuzzed-fuzzed.
- Charged pseudoscalar channel: A reliable unconstrained determination of the first excited state was not possible, but fixing it to 3 times the ground-state mass produced an acceptable fit.
- Charged pseudoscalar channel: Two-pion contamination effects are hardly detectable for the relevant mPS and t values despite the small statistical errors.The expected O(a2) contamination is negligible relative to the three-pion contribution in this regime.
- Neutral pseudoscalar channel: The neutral pseudoscalar analysis uses connected and disconnected correlators, a 4 × 4 correlator matrix, and blocked bootstrap errors based on 80-trajectory Monte Carlo segments.Ratios in Figure 6 display the disconnected contribution.
- Neutral pseudoscalar channel: Nonzero-momentum neutral-meson results agree with Table 13, while at aµq = 0.004 an energy of 0.309(27) corresponds to a mass 0.164+47The nonzero-momentum method avoids vacuum subtraction and provides a crosscheck.
- Neutral pseudoscalar channel: The coefficients are estimated as c = −5.0(1.2) at µq = 0.004 and c = −6.7(2.8) at µq = 0.0085, consistent within errors and with the first-order phase-transition scenario.
- Decay constants and renormalization: Neutral and charged decay constants are compatible within 1 to 1.5 standard deviations, unlike the roughly 50 MeV neutral–charged pseudoscalar mass difference.The neutral decay-constant comparison uses ZA = 0.76(2).
7 Chiral Perturbation Theory analysis of fPS and mPS
The paper applies chiral perturbation theory to pseudoscalar masses and decay constants while accounting for finite-size, autocorrelation, and cross-correlation effects. The analysis finds finite-size effects relevant mainly at the lightest mass and emphasizes that unknown NNLO terms can dominate the systematic uncertainty.
- Analysis strategy: The analysis describes fitting pseudoscalar masses and decay constants with one-loop χPT, while testing whether two-loop corrections affect the selected dataset.The two-loop investigation found the dataset insensitive to higher-loop corrections.
- Finite-size effects: Finite-size corrections are included because they significantly affect amPS and especially afPS at the lowest and next-to-lowest µq values.
- Analysis strategy: Continuum χPT is used to describe dependence on both finite spatial size L and bare quark mass µq despite a large neutral-pion O(a2) artifact.Lattice χPT and Symanzik analyses are cited as the justification.
- Analysis strategy: The χPT fit depends on four unknown parameters, B0, F, Λ3, and Λ4, determined from the data.
- Error analysis: Autocorrelations are handled by blocking, with 32-measurement blocks corresponding to more than 60 Monte Carlo trajectories and producing safely uncorrelated data.The blocked data are used for the covariance matrix, χ2, and fit-parameter errors.
- Error analysis: Suppressing off-diagonal covariance terms changes error bars only at the percent level, indicating negligible cross-correlations in the full dataset.This is consistent with fits using separate subgroups for masses and decay constants.
- Finite-size effects: The smallest mPSL exceeds 3, and finite-size effects are relevant only at aµq = 0.004; at larger masses they remain below statistical errors.The mPS deviations are barely larger than statistical errors of about 0.5%.
- Finite-size effects: Finite-size effects are below a few percent but remain comparable to statistical errors and can therefore affect the extracted low-energy constants.Higher-order finite-size formulae improve the calculation.
8 Summary
The paper documents practical methods for maximally twisted Nf = 2 simulations, correlator calculations, error analysis, scale setting, pseudoscalar observables, and χPT fits. These methods support precise measurements and provide a technical reference for ongoing lattice QCD work.
- Maximal-twist tuning: The authors justify tuning mPCAC/µq ≤0.1, with an error ∆mPCAC/µq ≤0.1, to achieve maximal twist.They argue that this precision leaves physical quantities with controlled O(a2) lattice artifacts.
- Charged correlators: Fuzzed stochastic time-slice sources combined with the one-end trick and random source locations significantly reduce noise in charged-meson two-point correlators.The demonstrated noise reduction applies at least to two-point correlators in the meson sector.
- Neutral correlators: Stochastic volume sources and variance-reduction methods enable quark-disconnected contributions to be computed with acceptable accuracy on unquenched configurations.The methods target the intrinsically noisy disconnected components of neutral-meson correlators.
- Algorithms and errors: Small enough autocorrelation times, together with Γ- and binning-method error analyses, permit trustworthy uncertainty estimates for physical observables.The paper describes both the Monte Carlo algorithm features and the procedures used to analyze correlated data.
- Scale determination: Better than 1% accuracy is reached for r0 in the chiral limit, while its quark-mass dependence is consistent with being quadratic in µq.The force parameter is used to check scaling toward the continuum limit.
- Physical observables: Effective-mass plots show stable Euclidean-time plateaux, supporting precise charged and neutral pseudoscalar masses and related quantities.The collected observables include the untwisted PCAC quark mass and the renormalization constant ZV.
- Chiral analysis: The χPT analysis estimates B0, F, Λ3, and Λ4 while examining higher-order stability and finite-size effects.The paper explains how errors on the fitted low-energy constants are obtained.
- Role of the paper: The methods are presented as a technical reference for the collaboration’s ongoing research using maximally twisted Wilson fermions.
Appendices
The appendices provide operator formulae in twisted and physical quark bases, explain their relation through the twist angle, and summarize their renormalization properties. Maximal twist simplifies most of these expressions.
- Operator representations: The appendices list bare quark bilinear operators relevant to the paper in the twisted quark basis and physical basis.The twisted basis is defined by the fermionic action used in the paper.
- Twist angle: The twist angle satisfies tan ω = µq/(m0 −mcrit), with mcrit determined by the maximal-twist tuning procedure.
- Basis transformation: The displayed operator expressions follow from the relation between twisted-basis χ and physical-basis ψ quark fields.
- Renormalization: All listed bare operators renormalize multiplicatively except P′3 and S′0, which mix additively with the identity.The divergences are cubic for P′3 and quadratic for S′0, with the latter vanishing as µq →0.
- Maximal twist: At maximal twist, ω = ±π/2, substantial simplifications occur in the operator and renormalization formulae.
B Evaluation of disconnected loops
The appendix develops stochastic estimators and variance-reduction techniques for quark-disconnected loops, whose intrinsic noise requires evaluation across time slices and many configurations. The combined methods can reduce stochastic noise below gauge-field fluctuations.
- Motivation: Disconnected correlator components are intrinsically noisier than connected components, so loops must be evaluated accurately across time slices and many gauge configurations.The target is to make stochastic error negligible relative to the intrinsic gauge noise.
- Stochastic estimators: Stochastic sources generally have whole-lattice support, and solving Mφ = ξ provides the fields used in loop estimators.Here M denotes the lattice Dirac matrix for a given flavour.
- Estimator construction: The loop estimator sums over colour, spin, and spacetime indices, while each disconnected loop uses independent stochastic samples and a single time-slice restriction.Independent samples avoid unwanted biases between the two loops arising from Wick contractions.
- Hopping-parameter method: The hopping-parameter method evaluates the first four terms of the hopping-parameter expansion of P XM−1 exactly on each gauge configuration.It reduces stochastic noise without much additional computational effort.
- Twisted-mass reduction: A twisted-mass variance-reduction identity applies to ΣX(1/Mu −1/Md), where an explicit aµq factor reduces fluctuations and the one-end trick avoids further inversions.
- Numerical effectiveness: At β = 3.9 and µq = 0.004, the equation (67) method gives an error 6 times smaller than a conventional stochastic volume source.The resulting stochastic contribution is negligible compared with intrinsic gauge noise.
- Numerical effectiveness: The applicable variance-reduction method makes stochastic noise smaller than intrinsic gauge-field noise in neutral-meson correlators.
C Γ-method and data-blocking
The appendix describes the Γ-method and data-blocking procedure used to estimate statistical errors for physical observables.
- Error estimation: The Γ-method and data-blocking procedure are used to estimate statistical errors of physical observables.
C.1 Γ-method
The Γ-method estimates errors for ensemble averages and nonlinear secondary observables by accounting for autocorrelations and propagating fluctuations through a linearized observable.
- C.1 Γ-method: The Γ-method estimates the standard deviation of an ensemble average from autocorrelation information.It uses estimators of the autocorrelation function and integrated autocorrelation time.
- C.1 Γ-method: The autocorrelation function estimator is formed from products of deviations of measurements from the true observable value.The measurement index labels individual samples, while the expectation brackets denote theoretical averaging.
- C.1 Γ-method: For nonlinear observables F = f(A), the method linearizes f around the true primary observables to estimate the variance of F.Hadron masses obtained from correlator fits are given as a typical secondary-observable example.
- C.1 Γ-method: The method defines a derived primary quantity as a linear combination of primary observables to propagate errors to nonlinear functions.This construction supplies the variance estimate for the secondary observable.
- C.1 Γ-method: The finite-statistics bias in the Taylor expansion is O(N^-1) and can be neglected when the number of measurements is sufficiently large.Replacing the true value A with the ensemble estimate introduces another bias of comparable order.
C.2 Binning method
The binning method accounts for autocorrelations by grouping consecutive measurements into blocks and estimating errors from the blocked data.
- C.2 Binning method: Data-blocking, or binning, uses a block size B analogous to the Γ-method window W when accounting for autocorrelations.For sufficiently large B, the integrated autocorrelation time can be estimated from blocked measurements.
- C.2 Binning method: The blocked-observable error σ̄F(B) is estimated with jackknife after measurements are grouped into blocks of size B.The resulting error estimate is used to infer autocorrelation effects.
C.3 Error on the error: Γ-method vs data-blocking √
The Γ-method and binning balance statistical and systematic errors differently when choosing a window or block size for autocorrelation analysis.
- C.3 Error on the error: Γ-method vs data-blocking √: The Γ-method balances a statistical error increasing with W/N against a systematic bias decreasing exponentially with W.The systematic bias is modeled as δsyst(σ̄a) ∼ 1/2 exp(−W/τ).
- C.3 Error on the error: Γ-method vs data-blocking √: An optimal Γ-method window Wopt can be selected from the onset of an error plateau or by minimizing the total error.The window bounds the autocorrelation sum.
- C.3 Error on the error: Γ-method vs data-blocking √: For binning, the statistical error increases with B/N while the systematic bias decreases with B^-1.Thus, block-size selection also requires balancing competing error contributions.
C.4 Further remarks
The authors compare Γ-method and binning-based error analyses, identify a limitation of binning for strongly autocorrelated PCAC masses, and describe smeared Wilson-loop construction for the static potential.
- C.4 Further remarks: In 15 analysed cases, the authors compare two Γ-method windowing criteria to test the robustness of their error estimates.One criterion follows an approximately optimal algorithm, while the other stops when the autocorrelation estimate becomes negative from fluctuations.
- C.4 Further remarks: For nonlinear observables, binning is combined with bootstrap or jackknife, with jackknife bin sizes chosen so the error stabilizes.The bootstrap uses B = 4, 8, 16, 32 trajectory units; jackknife targets Bopt/τint ≈ 10 or larger.
- C.4 Further remarks: Binning and Γ-method estimates are generally comparable, but binning underestimates the PCAC quark-mass error when significant autocorrelations are present.The authors attribute this to the Γ-method’s more favorable dependence on the number of measurements for autocorrelation errors.
- C.4 Further remarks: The analysis therefore uses the Γ-method for plaquette, amPCAC, and other quoted fermionic quantities, while binning gives similar results in most cases.The PCAC quark mass is the stated exception where binning was not reliable.
- C.4 Further remarks: Static-potential measurements use temporally smeared links and spatial APE smearing to improve Wilson-loop signals and reduce excited-state contamination.The spatial smearing preserves the stated gauge and discrete symmetries after SU(3) projection.
- C.4 Further remarks: The static quark-antiquark correlator matrix is built from smeared spatial strings and improved temporal links, using five string operators to form a 5 × 5 matrix.The construction corresponds technically to spatially smeared and temporally improved Wilson loops.