Performance of Low Mode Averaging on Twisted-Mass Fermion Ensembles at the physical pion mass point

A direct contraction of propagators built from $S_{\eta}^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}}}$ as in eq. 3 produces purely IR and purely UV contributions, along with several mixed IR–UV terms. In practice, this typically requires dedicated contraction kernels and considerable implementation effort. By contrast, the right-hand side (r.h.s) of eq. 6 allows for a much simpler procedure: compute the full stochastic correlator, subtract its stochastic IR component, and replace it with an improved—ideally exact—estimate. Moreover, even in lattice setups where the two sides of eq. 6 are not identical, the construction in the r.h.s. still provides a valid and more efficient realisation of LMA, with double counting of the IR contribution properly removed.

We note that both approaches, left- and right-hand sides of eq. 6, have been employed in the literature. In fact, the formulation based on the full $S_{\eta}^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}}}$ dates back to the original LMA works [49], while the second strategy was introduced in the context of combining LMA with AMA [22], as it naturally enables independent improvement of the IR contribution via LMA and the UV contribution via TSM. To our knowledge, however, it has not been emphasized that, in certain commonly used computational setups, the two procedures are in fact exactly equivalent¹¹1Actually, in parts of the literature, the full correlation function $C_{\eta}^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}}}$ is sometimes regarded as superior because the mixed terms are assumed to yield additional improvement since their IR part is computed exactly. This may indeed happen in setups where eq. 6 does not hold, e.g. as discussed in Ref. [51] for baryons computed with stochastic sources. For the cases considered here, however, we demonstrate that the two constructions are identical.. Recognising this equivalence or, more generally, adopting the second strategy in place of the first, has several important consequences:

Noise reduction

When eq. 6 holds, it makes explicit that LMA improves only the purely IR contribution $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}(t)$. All UV contributions, as well as any IR–UV mixing, are governed by stochastic noise. This observation is crucial for identifying the class of observables, Euclidean-time regions, and quark-mass regimes in which LMA can deliver significant variance reduction and meaningful performance gains. In particular, noise reduction is expected only where $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}(t)$ dominates the correlator. In general, this is not straightforward to identify, as the extent of this region depends sensitively on the number of deflated modes, as we will demonstrate with our numerical results.

Reduced contraction costs

Evaluating the r.h.s. of eq. 6 is considerably easier and more efficient than computing $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}}}_{\eta}$ through explicitly expanding it in terms of the IR and UV contractions. The correlators $C_{\eta}$ and $C_{\eta}^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$ are expressed in terms of ordinary stochastic contractions differing only in the contracted propagator, while $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$ is the only genuinely new object. Effectively, the $2^{n}$ IR-UV contraction terms—with $n$ being the number of propagators involved in the contraction—are replaced by only three, since all interference contributions are avoided.

Post-production LMA

If the stochastic sources $\eta$ are reproducible—i.e., the same exact stochastic source can be generated in independent runs by e.g. storing the point-source coordinates or the random seeds—then stochastic correlators can be improved a posteriori. Suppose $C_{\eta}$ has already been computed using the full stochastic propagator, and one subsequently wishes to apply LMA. Then $C_{\eta}^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$ and $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$ can be evaluated independently, without recomputing the full propagator, but using only the exact IR information.

Optimised production

The separation implied by eq. 6 naturally suggests an optimised computational strategy, particularly at physical quark masses. The evaluation of $C_{\eta}$ requires Dirac matrix inversions, for which modern multigrid solvers are optimal but exhibit poor strong scaling, and are therefore best run on as few nodes as possible. By contrast, the computation of the IR sector involves a large number of exact eigenvectors; while memory-intensive, it scales efficiently with the number of nodes. The correlators $C_{\eta}^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$ and $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$ can thus be computed in a separate stage on large node counts, decoupled from the full Dirac operator inversions, provided that the stochastic sources $\eta$ are reproducible. On the other hand, in the direct $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}}}_{\eta}$ construction, the UV projection forces inversions and IR modes to be treated simultaneously within a single workflow, often requiring execution on more nodes than optimal due to the increased memory demands.

In conclusion, we advocate implementing LMA as defined by the right-hand side of eq. 6, or better, its generalisation in eq. 79 of the conclusions, and in general ensuring that correlation functions are reproducible so that LMA can be applied a posteriori when needed.

We first briefly summarise the key properties of Wilson twisted-mass fermions that are relevant for discussing their spectrum and deriving the mass dependence of LMA. In the twisted basis, i.e., the one in which the mass parameter $\mu$ couples to a pseudoscalar density, the maximally-twisted Wilson–Dirac operator with a quark mass $\mu$ reads as²²2For simplicity in the following discussion, we set the Wilson parameter to $r=1$; the case $r=-1$ is completely equivalent.

$$D(\mu)=D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}+i\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\,\mu\mbox{\quad with\quad}D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}=\left[\gamma\cdot\tilde{\nabla}+W_{\mathrm{cr}}\right]\,,$$

(7)

where $D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}$ is the massless $\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}$-Hermitian Wilson–Dirac operator, consisting of the hopping term $\gamma{\cdot}\tilde{\nabla}$ and a diagonal mass term $W_{\rm cr}$ tuned to its critical value, i.e. to vanishing residual mass, with an optional clover term included. The associated Hermitian operator

$$Q_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}\equiv D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}=Q_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}^{\dagger}$$

(8)

has eigenpairs $\left(\Lambda_{j},v_{j}\right)$ satisfying

$$Q_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}v_{j}=\Lambda_{j}v_{j}\mbox{\quad with\quad}\Lambda_{j}\in\mathbb{R},.$$

(9)

It follows that the equivalent twisted-mass Dirac operator reads as

$$Q(\mu)\equiv D(\mu)\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}=Q_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}+i\mu$$

(10)

and since the mass term enters multiplied by the identity operator, it shares the same eigenvectors of $Q_{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}$, while its eigenvalues are shifted into the complex plane by $\mu$,

$$Q(\mu)v_{j}=\lambda_{j}(\mu)\,v_{j}\mbox{\quad with\quad}\lambda_{j}(\mu)\equiv\Lambda_{j}+i\mu\,.$$

(11)

These vectors are also eigenmodes of the positive-definite squared operator,

$$Q(\mu)^{\dagger}Q(\mu)\,v_{j}=Q(-\mu)Q(\mu)\,v_{j}=\left(\Lambda_{j}^{2}+\mu^{2}\right)v_{j}\qquad\Longrightarrow\qquad|\lambda_{j}(\mu)|\geq|\mu|{\rm~~for~any~}j\,.$$

(12)

The above relation highlights a central property of the twisted-mass operator, namely, the parameter $\mu$ provides an infrared cutoff that prevents the appearance of near-zero modes [43]. As a consequence, at non-zero mass, the squared twisted operator is always positive definite, with condition number $\kappa\geq\lambda_{\rm max}^{2}/\mu^{2}$, and its inverse is well defined.

Another important property of the twisted mass operator is that its spectrum is independent of the value of $\mu$, being entirely determined by the massless Hermitian Wilson operator $Q_{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}$. As a result, the cost of the eigensolver does not depend on the quark mass, and the same set of deflated modes can be employed for any value of the twisted mass parameter $\mu$. This contrasts with, e.g., the Wilson operator and other discretisations, where the eigenvectors must be recomputed for each quark mass. The inverse of the twisted mass Dirac operator then follows directly from the spectral decomposition of $Q_{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}$ as

$$Q_{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}=D_{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}=\sum_{j=0}^{N_{\rm ev}}\Lambda_{j}\,v_{j}v^{\dagger}_{j}\qquad\Longrightarrow\qquad S(\mu)=D_{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}}^{-1}(\mu)=\sum_{j=0}^{N_{\rm ev}}\frac{1}{\lambda_{j}(\mu)}\,\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}v_{j}v^{\dagger}_{j}\qquad\Longrightarrow\qquad S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}=\sum_{j=0}^{N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\rm ev}}\frac{1}{\lambda_{j}}\,\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}v_{j}v^{\dagger}_{j}\,,$$

(13)

where the modes are ordered by magnitude $|\Lambda_{j}|$ and the IR propagator, $S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}$, is obtained by restricting the sum to the first $N_{\rm ev}^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}$ modes.

We now examine the spectral densities of the massless Wilson Dirac operator, as well as those of the twisted-mass operator at the light quark mass. It has been shown in Refs. [19, 55, 59, 49, 50] that, on sufficiently large-volume ensembles—i.e., where the spectrum is dense and finite-volume effects are negligible—the spectral density of the massless Wilson Dirac operator, $\rho_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}$, integrated over an interval $\left[\Lambda,\Lambda{+}\delta\Lambda\right]$ with $\Lambda$ in the IR region, is proportional to both the interval width $\delta\Lambda$ and the physical volume $V$. Consequently, once these factors are factored out, the normalised density is approximately constant, namely

$$f_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}\left(\Lambda,\delta\Lambda;V\right)\equiv\frac{1}{V\,\delta\Lambda}\int_{\Lambda}^{\Lambda+\delta\Lambda}\!\!\!\!\rho_{\!{\textstyle\mathstrut}{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}}(\tilde{\Lambda};V)\differential{\tilde{\Lambda}}\quad=\quad\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}+df_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}(\Lambda,\delta\Lambda;V)\,,$$

(14)

where $V$ is physical volumes, $\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}$ is a constant related to the chiral condensate $\Sigma$ in the chiral limit [19], see section III.4, and $df_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}$ is a smooth function, providing subleading corrections across physically relevant IR eigenvalue intervals.

We confirm this behaviour using six near-physical-pion-mass $N_{f}{=}2+1+1$ ensembles generated by the Extended Twisted Mass Collaboration (ETMC) at four different lattice spacings and three different physical volumes (see Ref. [64]). Our results are summarised in table 1 and shown in fig. 1. In the left panel, we show the dependence of the normalised density as a function of the eigenvalue norm $\Lambda$. For our ensembles, we determine $\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}$ by performing a linear fit of $f_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}\left(\Lambda,\delta\Lambda;V\right)$ in the $\Lambda/\mu_{\ell}$ interval $\left[1.5,5\right]$ with $\mu_{\ell}$ the simulated light quark mass and $\delta\Lambda/\mu_{\ell}=0.5$. We observe that its values have a mild dependence on the lattice spacing and are approximately constant with respect to the volume.

Given the smoothness of the massless Wilson operator, the density of the spectrum of the twisted operator follows immediately. Wilson eigenvalues with $|\Lambda_{j}|\leq\mu$ are mapped into the region

$$\mu\leq\left(|\lambda_{j}(\mu)|=\sqrt{\Lambda_{j}^{2}+\mu^{2}}\right)\leq\sqrt{2}\,\mu\,,$$

(17)

producing an accumulation of modes in the interval $[\mu,\sqrt{2}\mu]$. For eigenvalues much larger than $\mu$, the spectrum resumes its approximately linear growth in $|\Lambda|$. This pattern is illustrated in the right panel of fig. 1 at the physical light-quark mass for each ensemble given in Table 1. This leads directly to the central observation of this section. As the quark mass $\mu$ increases, all Wilson modes with $\absolutevalue{\Lambda_{j}}\leq\mu$ are compressed into a narrow region of the twisted-mass spectrum and therefore contribute uniformly to its infrared part. In table 1, we report the measured and predicted value for $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}}_{\rm ev}(\sqrt{2}\mu;V)$ according to

$$N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}}_{\rm ev}(\sqrt{2}\mu;V)\equiv\int_{\mu}^{\sqrt{2}\mu}\!\!\!\!\!\!\rho_{\!{\textstyle\mathstrut}{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{TM}$}}}}(\absolutevalue{\lambda};V,\mu)\,\differential{\absolutevalue{\lambda}}\quad=\quad\int_{0}^{\mu}\!\!\!\!\rho_{\!{\textstyle\mathstrut}{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}}(\Lambda;V)\differential{\Lambda}\quad\equiv\quad N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}_{\rm ev}(\mu;V)\quad\simeq\quad\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}\,\mu\,V\,.$$

(18)

Indeed, given the above considerations, the integrated density of eigenvalues in $[\mu,\sqrt{2}\mu]$ of eq. 18 grows linearly with $\mu$.

The previous analysis of the twisted mass operator shows that to keep constant LMA performances, i.e. to deflate exactly a given region of the spectral density of the operator, the cost grows linearly with the quark mass and with the volume. This occurs because all infrared modes contribute comparably due to their similar eigenvalue magnitude induced by the mass shift and increased density with the volume. This scaling can be expressed as

$$N^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}_{\rm ev}\;\propto\;\mu\,V\,.$$

(19)

Although this argument is derived for the twisted mass operator, where the mass shift explicitly accumulates eigenvalues, the conclusion is not specific to this regularisation. For instance, in the Wilson case, the eigenvalues are not similarly clustered but still populate the region near the mass shift threshold [59]. Since different lattice discretisations describe the same infrared physics, the density of modes relevant for low-energy observables must exhibit the same parametric dependence. Therefore, it is natural to expect that the linear scaling with the multiplicatively renormalisable quark-mass parameter (here $\mu$) and the spacetime volume $V$ holds more generally across fermion formulations. This leads to several important consequences as discussed below.

LMA computational cost and memory usage grow as $\mu\,V^{2}/a^{4}$

Given the scaling of $N^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}_{\rm ev}\propto\mu V$ and considering that the cost of lattice operations scales at least with the number of lattice points $V/a^{4}$, the computational cost and memory footprint of exact LMA scale as $\mu\,V^{2}/a^{4}$. This $V^{2}$ scaling renders the method rapidly impractical on large volumes due to both larger resource requirements and poor parallel scaling of algorithms, in particular, multigrid solvers. The LMA further deteriorates as the quark mass increases. In particular, if $N^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}_{\rm ev}$ modes are required to achieve a given performance at the light-quark mass $\mu_{\ell}$, maintaining the same efficiency would require roughly $\mu_{s}/\mu_{\ell}\approx 27$ times more modes at the strange mass and $\mu_{c}/\mu_{\ell}\approx 324$ times more at the charm mass. Given that already at the physical light-quark mass on a $(5\,{\rm fm})^{4}$ volume several hundred modes are needed for a significant gain, extending exLMA to heavier quarks quickly becomes computationally prohibitive. These considerations imply that LMA is most effective at light-quark masses, where the number of required modes remains manageable, and the balance between cost and variance reduction is most favourable.

Contraction costs for all-to-all correlation functions scale with $(N^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}_{\rm ev})^{n}$

Contraction costs can also become rapidly prohibitive with the number of propagators, $n$, involved. While exact eigenvectors enable the construction of all-to-all correlation functions, the associated cost typically scales as $(N^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}_{\rm ev})^{n}$, leading to a steep increase already for a moderate number $n$ of propagators—see e.g. fig. 9 for baryon contractions. As our results will demonstrate, this scaling effectively renders all-to-all baryon correlators ($n=3$) already impractical for light-quark masses and moderate volumes, and, as volume or quark mass increase, fully exact all-to-all correlation functions become generically unfeasible, except for simpler observables, such as single-quark loops with $n=1$. In practice, this limitation can be alleviated by replacing all-to-all contractions with stochastic estimators, since the signal typically saturates at moderate statistics. Other possible strategies not considered here include sparsening [36, 57, 30].

Gain of LMA versus stochastic approaches decreases with $1/\mu^{2}$

Assuming that a stochastic propagator samples the Dirac spectrum approximately uniformly, then the fraction of infrared modes it captures increases with the quark mass, since the IR region relevant for long-distance physics grows proportionally to $\mu$ for a given spectrum. Stochastic estimators, therefore, become progressively more efficient at sampling the physically relevant modes as the mass increases. At the same time, the cost of LMA rises linearly with $\mu$, as the number of eigenmodes that must be treated explicitly increases accordingly. Thus, the two approaches scale in opposite manners, with LMA becoming more expensive, while stochastic estimators improve in quality. As a result, the relative advantage of LMA over stochastic methods rapidly decreases with increasing quark mass, scaling approximately as $1/\mu^{2}$. Consequently, for sufficiently heavy quarks, stochastic approaches become markedly more efficient than LMA.

LMA performance is limited by the heaviest quark mass in the correlation function

The previous remarks apply to correlators with a single quark mass. A more subtle situation arises for mixed correlators involving both light and heavy quarks, with masses $\mu_{\ell}$ and $\mu_{h}$, such as the kaon two-point function. Suppose one deflates $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\ell}$ modes that are determined to be sufficient for achieving a good performance at $\mu_{\ell}$, and does not scale this number for the heavier propagator, i.e., instead of the larger value $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{h}=(\mu_{h}/\mu_{\ell})\,N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\ell}$ one still uses $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\ell}$—this is particularly advantageous for the twisted mass discretization, where the same eigenvectors can be used for all quark masses, see eq. 13. Then the corresponding propagators can be written as

$$S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}},\ell}_{\eta}=S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell}_{\rm ex}+S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{UV}$}}},\ell}_{\eta}\mbox{\quad and\quad}S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}},h}_{\eta}=S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h}_{\rm ex/\eta}+S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{UV}$}}},h}_{\eta}=S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h\subset\ell}_{\rm ex}+S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h\not\subset\ell}_{\eta}+S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{UV}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{UV}$}}},h}_{\eta}\,,$$

(20)

where the infrared and ultraviolet parts are separated and labelled according to whether they are treated exactly (ex) or stochastically ($\eta$). For the heavy propagator, the full infrared region of size $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{h}$ is split into an exact part for the first $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\ell}$ modes, $S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h\subset\ell}_{\rm ex}$, and a stochastic remainder, $S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h\not\subset\ell}_{\eta}$.

Considering then the infrared contribution to a bilinear correlator, one obtains

$$C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell h}_{\rm ex/\eta}\equiv C(S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell}_{\rm ex},S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h}_{\rm ex/\eta})=C(S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell}_{\rm ex},S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h\subset\ell}_{\rm ex})+C(S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell}_{\rm ex},S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},h\not\subset\ell}_{\eta})=C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell(h\subset l)}_{\rm ex}+C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell(h\not\subset l)}_{\eta}\,,$$

(21)

where, in the first equality, the IR contraction is split into two contributions (one purely exact and one mixed); and then, in the second equality, the mixed exact–stochastic term is effectively rendered stochastic since we assume eq. 6 to hold. Consequently, only the light–light part of the correlator benefits from LMA, while the remaining contribution does not. Namely, using eq. 6, the full correlator can be written as

$$C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}},\ell h}_{\rm ex/\eta}=C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}},\ell h}_{\eta}-C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell(h\subset l)}_{\eta}+C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell(h\subset l)}_{\rm ex}\,,$$

(22)

making explicit that only the fully exact part is improved. This means that if all $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{h}$ eigenvalues contribute approximately equally, the improved fraction of the correlator scales as $\mu_{\ell}/\mu_{h}$, since we expect $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\ell}\times N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{h}$ modes to contribute in $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell h}$, but only $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\ell}\times N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\ell}$ have been deflated in $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}},\ell(h\subset l)}$. For example, for a kaon at physical quark masses, this yields a suppression of order $1/27$ of the noise reduction, significantly reducing the overall effects of LMA.

The above considerations strongly indicate that, in modern lattice QCD applications, LMA is beneficial primarily for correlation functions composed exclusively of light-quark propagators, which thus define the focus of this work.

As a by-product of the analysis presented in section III.2, the extracted values of $\rho_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}\equiv Z_{P}\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}$ enable a determination of the renormalised chiral condensate $\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}$ via the Banks–Casher relation [19]. The chiral condensate $\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}$ is related to the disconnected contribution of the scalar insertion, which can be written as

$$\expectationvalue{\left(u\bar{u}+d\bar{d}\right)_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}}(\mu)\equiv Z_{P}\expectationvalue{\left(u\bar{u}+d\bar{d}\right)_{\mathrm{sub}}}(\mu)=-\frac{2\mu}{V}\Tr[D^{-1}(+\mu)D^{-1}(-\mu)]=-\frac{2\mu}{V}\sum_{i}\frac{1}{\Lambda^{2}_{i}+\mu^{2}}\,,$$

(23)

where $\expectationvalue{\left(u\bar{u}+d\bar{d}\right)_{\mathrm{sub}}}$ denotes the properly subtracted (mixing with $\mathds{1}a^{-2}\mu$ and $\mathds{1}\mu^{3}$ removed) scalar light-quark density.

In the infinite-volume limit, the sum can be replaced by an integral. Assuming a constant spectral density up to a cutoff $\Lambda_{\max}$, one obtains

$$\expectationvalue{\left(u\bar{u}+d\bar{d}\right)_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}}(\mu)\simeq-2\,Z_{P}\,\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}\int_{0}^{\Lambda_{\max}}\!\!\!\frac{\mu}{\Lambda^{2}+\mu^{2}}\,\differential\Lambda=-\pi\,Z_{P}\,\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}+\order{\frac{Z_{P}\,\mu}{\Lambda_{\max}}}\,,$$

(24)

where we used

$$\int_{0}^{\Lambda_{\max}}\!\!\!\frac{\mu}{\Lambda^{2}+\mu^{2}}\,\differential\Lambda=\arctan{\frac{\Lambda_{\max}}{\mu}}=\frac{\pi}{2}+\order{\frac{\mu}{\Lambda_{\max}}}\,.$$

(25)

The renormalised chiral condensate is then defined as

$$\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}=-\frac{1}{2}\lim_{\mu\to 0}\expectationvalue{\left(\bar{u}u+\bar{d}d\right)_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}}(\mu)=\frac{\pi}{2}Z_{P}\,\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}^{\,0}\,,$$

(26)

where $Z_{P}\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}^{\,0}$ denotes the spectral density in the chiral limit, i.e. for vanishing sea and valence quark masses. Owing to the relation between the spectral density of the twisted-mass operator and that of the massless Wilson operator, the extracted $\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}$ are already in the valence chiral limit, but correspond to a finite sea quark mass, $m_{\pi}\simeq 140\,$MeV.

Taking the continuum limit of $\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}$ and comparing with independent determinations of $\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}$, we find values in the expected range, see fig. 3, albeit in mild tension with our earlier result from ETMC21 [2]. In that work, $\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}$ was extracted using chiral perturbation theory relations for the pion mass across multiple ensembles, including heavier pion masses and a controlled extrapolation to the chiral limit. The observed discrepancy can therefore be attributed to residual sea-quark mass effects in the present determination, since other lattice errors appear to be well under control, including finite size corrections. Indeed, by comparing three significantly different volumes for the B ensemble yields consistent spectral densities, see table 1.

To quantify these residual sea-quark mass effects, one may employ the chiral expansion of the subtracted quark scalar operator vacuum expectation value [48],

$$\expectationvalue{\left(u\bar{u}+d\bar{d}\right)_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}}(\xi_{\pi})=-2\,\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}\left(1+\xi_{\pi}(4\bar{h}_{1}-\bar{l}_{3})+\mathcal{O}(\xi_{\pi}^{2})\right)\mbox{\quad with\quad}\xi_{\pi}=\frac{m^{2}_{\pi}}{16\pi^{2}f_{\pi}^{2}}\,,$$

(27)

where $\bar{l}_{3}=3.5(3)$ [12] is the next-to-leading-order low-energy constant governing the quark-mass dependence of $m_{\pi}^{2}$, while $\bar{h}_{1}$ denotes a less constrained contact term entering scalar correlation functions. Owing to the relation between the scalar insertion and the spectral density via the trace of modes, this expression directly translates into

$$\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}(\xi_{\pi})=Z_{P}\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}(\xi_{\pi})=\frac{2\,\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}}{\pi}\Big(1+\xi_{\pi}(4\bar{h}_{1}-\bar{l}_{3})+\mathcal{O}(\xi_{\pi}^{2})\Big)\qquad\Longrightarrow\qquad 4\bar{h}_{1}-\bar{l}_{3}=\frac{1}{\xi_{\pi}}\left(\frac{\pi\,\bar{\rho}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}}{2\,\Sigma}-1\right)\,.$$

(28)

This relation suggests two possible strategies. One may either adopt an external input for $\bar{h}_{1}$ to correct for pion-mass effects in the spectral density, or, alternatively, determine $\bar{h}_{1}$ directly from our spectral-density data combined with an independent continuum limit determination of $\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}$ from ETMC21 [2].

Regarding the first approach, a prediction for $H_{1}^{r}$ in three-flavour chiral perturbation theory at $\mu=0.77~\mathrm{GeV}$ is available in Ref. [10], determined in the large $N_{c}$ limit. Converting this to the scale-independent $\bar{h}_{1}$ [48], and assigning a conservative uncertainty of approximately $30\%$ [21], one obtains $\bar{h}_{1}\simeq 4.5(1.3)$. The corresponding $4\bar{h}_{1}-\bar{l}_{3}$ value, represented as a blue dot in fig. 3, can then be used to correct the chiral condensate for chiral effects. The resulting continuum-limit result is

$$\sqrt[3]{\Sigma_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{R}$}}}}}=269.5(4.6)~\mathrm{MeV}\qquad\text{in}~~N_{f}{=}2{+}1{+}1~\text{at 2 GeV scale}\,,$$

(29)

as shown in fig. 3 with very good agreement with our previous ETMC21 determination.

Following the second strategy instead, by comparing to ETMC21 to determine the chiral dependence, we find

$$4\bar{h}_{1}-\bar{l}_{3}=17.4(4.3)\qquad\Longrightarrow\qquad\bar{h}_{1}=5.2(1.1)\,.$$

(30)

We consider a generic meson two-point correlation function, computed using a backwards-running propagator via $\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}$-hermiticity, $S=\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}S^{\dagger}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}$³³3In the case of TM fermions, one has $\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}S^{\dagger}_{r}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}=S_{-r}$, where $r=\pm 1$ is the Wilson parameter. In this section, we adopt the canonical quark basis, where the $f$-quark mass term is $\mu\bar{f}f$ and the $f$-quark reads $S_{r}=\exp{i\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\frac{\pi}{4}}\left[D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}(r)+i\mu\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\right]\exp{-i\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\frac{\pi}{4}}$. , defined as

$$C(t_{1},t_{2},{\bf\it p}_{1},{\bf\it p}_{2},\Gamma_{1},\Gamma_{2})=\Tr\left[T_{{\bf\it p_{2}}}^{t_{2}}\,\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\Gamma_{2}\,S\,T_{{\bf\it p_{1}}}^{t_{1}}\,\Gamma_{1}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\,S^{\dagger}\right]\mbox{\quad with\quad}\,T_{{\bf\it p_{i}}}^{t_{i}}\,({\bf\it x},t)=\delta_{tt_{i}}e^{i{\bf\it x}\cdot{\bf\it p}_{i}},$$

(31)

where $x=({\bf\it x},t)$ denotes the spacetime coordinates. The operator $T_{{\bf\it}}\,\!\!$ acts on the lattice coordinates, projecting the propagator onto the timeslice $t_{i}$ and momentum ${\bf\it p}_{i}$. While the two propagators $S$ may in general have different masses, here we restrict to light quarks, where the benefits from LMA are greatest as discussed in the previous section. For the twisted-mass operator, the propagators can be generated with either equal or opposite signs of the Wilson parameter $r=\pm 1$ multiplying the light-quark twisted mass $\mu$, hence $\pm\mu$ in eq. 7, corresponding either to twisted-mass (TM) or Osterwalder–Seiler (OS) regularisations. Both setups are tested, yielding comparable improvements. Therefore, results are presented for TM or OS regularisation, as specified below.

Regarding momentum boosts, we restrict to the zero-boost case, ${\bf\it p}_{1}={\bf\it p}_{2}={\bf\it 0}$, since non-zero momenta would require dedicated inversions on momentum stochastic sources. Concerning the Dirac structures, several choices vanish at zero momentum, in particular those with unmatched spatial Dirac indices. Also, channels with pseudoscalar or scalar structures ($\Gamma_{i}=\gamma_{5}$ or ${\bf 1}$) rapidly saturate the stochastic noise, and the correlation functions are subject only to gauge noise. In these cases, no significant improvement from LMA is observed, as the correlator on a single configuration is already saturated with a small number of stochastic sources. In the following, we therefore restrict to non-vanishing vector, axial, and tensor structures, where a significant noise-to-signal ratio is observed.

The correlation functions have been computed either stochastically or via all-to-all methods exploiting the exact low-mode eigenvectors, as follows from eq. 6. For the stochastic part, we employ momentum-projected stochastic sources with time and spin dilution,

$$\eta_{\tau,{\bf\it p},\mu}(x,\nu,a)=\delta_{t\tau}\delta_{\mu\nu}e^{\tfrac{i}{2}{\bf\it x}\cdot{\bf\it p}}\,\eta({\bf\it x},a)\mbox{\quad with\quad}\frac{1}{N_{\eta}}\sum_{\eta}\eta^{*}({\bf\it x},a)\eta({\bf\it y},b)\approx\delta_{ab}\delta_{{\bf\it x}{\bf\it y}},$$

(32)

where $\mu,\nu$ are spin indices, and $a,b$ colour indices. Subscripts label independent sources, each requiring a separate inversion. In practice, we invert over all four Dirac components $\mu$, while $\tau$ is varied together with a new stochastic sample $\eta$, and we take ${\bf\it p}={\bf\it 0}$. These sources allow us to approximate the term

$$T_{{\bf\it p_{1}}}^{t_{1}}\,\Gamma_{1}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\approx\frac{1}{N_{\eta}}\sum_{\eta}\sum_{\mu\nu}\eta_{t_{1},{\bf\it p}_{1},\mu}\;(\Gamma_{1}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}})_{\mu\nu}\;\eta^{\dagger}_{t_{1},{\bf\it p}_{1},\nu}\,$$

(33)

in eq. 31 leading to

$$C_{\eta}(t_{1},t_{2},{\bf\it p}_{1},{\bf\it p}_{2},\Gamma_{1},\Gamma_{2})=\sum_{\mu\nu}(\Gamma_{1}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}})_{\mu\nu}\,\phi^{\dagger}_{\eta,t_{1},{\bf\it p}_{1},\nu}\,T_{{\bf\it p_{2}}}^{t_{2}}\,\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\Gamma_{2}\,\phi_{\eta,t_{1},{\bf\it p}_{1},\mu}\mbox{\quad with\quad}\phi_{\eta,t_{1},{\bf\it p}_{1},\mu}=S\,\eta_{t_{1},{\bf\it p}_{1},\mu}\,,$$

(34)

which corresponds to an inner product of the propagated stochastic sources $\phi_{\eta}$. The stochastic IR contribution is obtained analogously as

$$C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\eta}\mbox{\quad same as~\lx@cref{creftype~refnum}{eq:stoch_corr_bil} replacing $\phi_{\eta}$ with\quad}\phi^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\eta,t_{1},{\bf\it p}_{1},\mu}=S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\,\eta_{t_{1},{\bf\it p}_{1},\mu}=\sum_{j=0}^{N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\rm ev}}\frac{1}{\lambda_{j}}\,\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}v_{j}v^{\dagger}_{j}\,\eta_{t_{1},{\bf\it p}_{1},\mu}\,,$$

(35)

where $\eta$ is exactly the same source, but propagated using only the IR part of the spectrum.

The remaining ingredient is the exact all-to-all IR contribution. Explicitly inserting the eigenmode decomposition of $S$ given in eq. 13, in the correlation function yields

$$C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}(t_{1},t_{2},{\bf\it p}_{1},{\bf\it p}_{2},\Gamma_{1},\Gamma_{2})=\sum_{j,k}^{N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{\rm ev}}\frac{1}{\lambda_{j}\,\lambda_{k}^{*}}\Tr\left[T_{{\bf\it p_{2}}}^{t_{2}}\,\Gamma_{2}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\,v_{j}v_{j}^{\dagger}\,T_{{\bf\it p_{1}}}^{t_{1}}\,\Gamma_{1}\gamma_{\scalebox{0.9}{$\scriptstyle 5$}}\,\,v_{k}v_{k}^{\dagger}\right]\,.$$

(36)

From a computational perspective, it is more efficient to evaluate inner products of eigenvectors, obtaining

$$C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}(t_{1},t_{2},{\bf\it p}_{1},{\bf\it p}_{2},\Gamma_{1},\Gamma_{2})=\sum_{j,l=1}^{N_{ev}^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}\frac{1}{\lambda_{j}\,\lambda_{k}^{*}}E_{jk}(t_{1},{\bf\it p}_{1},\Gamma_{1}\gamma_{5})E_{kj}(t_{2},{\bf\it p}_{2},\gamma_{5}\Gamma_{2})\mbox{\quad with\quad}E_{jk}(t,{\bf\it p},\Gamma)=v_{j}^{\dagger}\,T_{{\bf\it p}}^{t}\,\,{\Gamma}\,{v_{k}}$$

(37)

requiring only $N_{ev}^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}(N_{ev}^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}+1)/2$ independent contractions to be computed by exploiting

$$E_{kj}(t,{\bf\it p},\Gamma)=E_{jk}^{*}(t,-{\bf\it p},\Gamma^{\dagger})\,,$$

(38)

and then assembling the all-to-all IR correlation functions during post-processing.

The most critical aspect in optimising LMA performance is determining how many modes to deflate. As discussed in section III.2, for the twisted mass operator, one needs at least as many modes as those within the region $[\mu,\sqrt{2}\mu]$, since all eigenvalues in this interval are of comparable magnitude. For the B64 ensemble, which is used for the testing in this section, this corresponds to roughly 45 modes at the simulated light-quark mass, as reported in table 1.

In algebraic multigrid, a coarse operator $D_{c}$ is defined in terms of prolongation $P$ and restriction $R$ operators as

$$D_{c}\equiv RDP\,,$$

(42)

where $D$ denotes the fine operator. In lattice QCD, $D_{c}$ is also referred to as the little Dirac operator [59]. While this construction generalises straightforwardly to multiple levels, here we restrict the discussion to the finest and coarsest levels only. Accordingly, $D_{c}$ denotes the operator on the coarsest grid, and $R$ and $P$ represent the composite restriction and prolongation operators mapping directly between the finest and coarsest grids.

To define multigrid-based projectors, one typically requires that $R$ and $P$ preserve the identity,

$$\mathbbm{1}_{c}\equiv R\mathbbm{1}P\,,$$

(43)

with $\mathbbm{1}_{c}$ the identity on the coarse space. This leads to the following projectors acting on the fine level:

$$\mathbb{P}\equiv PR\,,\qquad\mathbb{P}_{D,r}\equiv PD_{c}^{-1}RD\,,\qquad\mathbb{P}_{D,l}\equiv DPD_{c}^{-1}R\,,\mbox{\quad and\quad}\mathbb{P}_{D,lr}\equiv\sqrt{D}PD_{c}^{-1}R\sqrt{D}\,,$$

(44)

The operator $\mathbb{P}$ simply projects the identity onto the coarse subspace, while the remaining operators incorporate the coarse inverse and define suitable projections for the quark propagator. For instance, using $\mathbb{P}_{D,r}$ one can decompose the propagator $S=D^{-1}$ as

$$S=\mathbb{P}_{D,r}S+(\mathbbm{1}-\mathbb{P}_{D,r})S\qquad\text{with}\qquad\mathbb{P}_{D,r}S=PD_{c}^{-1}RDD^{-1}=PD_{c}^{-1}R\equiv S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}\,.$$

(45)

An equivalent representation follows from the other projectors, finding

$$S-S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}=(\mathbbm{1}-\mathbb{P}_{D,r})S=S(\mathbbm{1}-\mathbb{P}_{D,l})=\sqrt{S}(\mathbbm{1}-\mathbb{P}_{D,lr})\sqrt{S}\,.$$

(46)

The connection with exact deflation becomes explicit if $R$ and $P$ are chosen as the left and right IR eigenvector matrices of $D$, respectively:

$$R=U^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\,,\qquad P=V^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\,,\mbox{\quad and\quad}\Lambda^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\equiv U^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}DV^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\qquad\Longrightarrow\qquad\mathbb{P}=V^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}U^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\,,\qquad D_{c}=\Lambda^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\,,\mbox{\quad and\quad}S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}=S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}\,.$$

(47)

Thus, the multigrid propagator $S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}$ approaches the IR propagator $S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}$ to the extent that $P$ and $R$ accurately capture the low-mode subspace. When specialising this general AMG framework to lattice QCD, additional symmetry-preserving structures are typically imposed:

• •

  To preserve the sparsity of the Dirac operator, $P$ and $R$ act on space-time aggregates of the lattice, ensuring that $D_{c}$ retains a nearest-neighbor structure. This is also in line with the property of local coherence [59], namely, eigenvectors restricted to aggregates still maintain significant overlap with the low-mode subspace.

• •

  To preserve $\gamma_{5}$-hermiticity and related symmetries, $P$ and $R$ are constructed to act separately on left- and right-handed spin components, enforcing $\gamma_{5}$-compatibility:

  $$\gamma_{5,c}=R\gamma_{5}P\,,\qquad\gamma_{5,c}R=R\gamma_{5}\,,\mbox{\quad and\quad}\gamma_{5}P=P\gamma_{5,c}\,.$$

  (48)

  This allows one to relate left and right eigenspaces and consistently choose

  $$R=P^{\dagger}\qquad\Longrightarrow\qquad D_{c}=P^{\dagger}DP\,.$$

  (49)

  It also follows that $\gamma_{5}$-hermiticity is preserved at the coarse level,

  $$D_{c}^{\dagger}\quad=\quad P^{\dagger}D^{\dagger}P\quad=\quad P^{\dagger}\gamma_{5}D\gamma_{5}P\quad=\quad\gamma_{5,c}P^{\dagger}DP\gamma_{5,c}\quad=\quad\gamma_{5,c}D_{c}\gamma_{5,c}\,,$$

  (50)

  allowing the definition of a Hermitian coarse operator

  $$Q_{c}\equiv D_{c}\gamma_{5,c}=Q_{c}^{\dagger}\,.$$

  (51)

• •

  The $\gamma_{5}$-compatibility is particularly crucial for twisted-mass fermions [4], as it ensures that the twisted-mass term retains its linear form at all multigrid levels:

  $$D_{c}(\mu)\equiv P^{\dagger}(D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}+i\mu\gamma_{5})P=P^{\dagger}D_{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}P+i\mu\gamma_{5,c}=D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}},c}+i\mu\gamma_{5,c}\,,$$

  (52)

  where $D_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{W}$}}}}$ denotes the massless Wilson operator. The coarse operator, therefore, inherits the same structure as the fine twisted-mass operator, i.e., the twisted mass $\mu$ provides a low-mode cutoff, while the eigenvectors remain those of the massless Wilson operator. In practice, this allows a prolongation operator $P$ constructed at a given $\mu$ to be reused for both signs and arbitrary values of $\pm\mu$, which is advantageous when employing multigrid as a preconditioner for twisted-mass simulations and observables [4, 5, 15].

The properties discussed above closely mirror those of exact LMA, but within a more general multigrid framework. It is important to identify, however, where the construction of multigrid preconditioning deviates from the construction of a multigrid-based LMA approach. In particular:

Multigrid preconditioning is dynamic

A central feature of multigrid solvers for lattice QCD is their adaptivity. In a typical multigrid correction, the inverse of $D_{c}$ is computed only to very low accuracy, and this approximate correction is then combined with a smoother. Additionally, the solver dynamically cycles among levels, often via K-cycles, using nested Krylov methods to determine whether to recurse further to coarser grids or return to finer levels. This adaptivity optimises the computational effort by adjusting the coarse-level work according to the current condition number and convergence state. As a result, the effective multigrid correction operator is not fixed, but evolves during the solve. Such strategies necessitate combining the preconditioner with a flexible Krylov solver.

LMA requires static operators

In contrast, LMA demands a strictly static construction where all sources of adaptivity must be removed so that the multigrid-based IR propagator defines a fixed operator. In particular, the inverse of $D_{c}$ must be computed to high (effectively exact) precision, rather than approximated. This is especially challenging for twisted-mass fermions, where the coarse operator is typically highly ill-conditioned; even achieving a relatively loose tolerance (e.g. a $10^{-1}$ relative residual) can already require hundreds of iterations [4]. A practical resolution is to perform exact deflation on the coarsest grid, which is the strategy adopted here. We define

$$S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}},{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}=PS^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{c}P^{\dagger}\qquad\text{with}\qquad S^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}_{c}=\sum_{j=0}^{N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{\rm ev}}\frac{1}{\lambda_{c,j}}\,\gamma_{5,c}v_{c,j}v^{\dagger}_{c,j}\mbox{\quad and\quad}Q_{c}v_{c,j}=\lambda_{c,j}v_{c,j}\,$$

(53)

so that the exact LMA construction is carried out on the coarsest grid. This significantly reduces the cost, both because of the smaller system size and because fewer eigenvectors are required, with the prolongation operator effectively lifting the coarse IR subspace to the fine level.

Different hyper-parameters and tuning

As a consequence, the set of hyperparameters and their tuning strategy differs substantially from standard multigrid preconditioning. Since LMA excludes any dynamic components, only parameters entering the definition of the coarse operator remain relevant, i.e., the aggregation size between levels, the number of null vectors, and the number of IR modes treated exactly on the coarse grid, $N^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{\rm ev}$. As we will see, in practice, the optimal values of these parameters also differ significantly from the corresponding ones that are optimal in multigrid preconditioning. Additionally, the realisation of optimal multigrid LMA will differ significantly with the kind of fermion discretisation, in the same way multigrid preconditioning does for Wilson [26, 13, 46], twisted-mass [4], overlap [27], domain-walls [28], and staggered [29] fermions. Thus, each of these requires dedicated studies.

As given in eq. 53, our strategy combines multigrid operators with exact deflation to construct the propagator $S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}},{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$. A first consequence of this choice is that the identity in eq. 6 no longer holds. This is due to the fact that $S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}},{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$ involves a nested projection, whereas the identity would apply to either the full $S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}}$ or the exact $S^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{IR}$}}}}$. Nevertheless, eq. 6 can still be used as a guiding principle to define an improved estimator by embedding it in a control variates framework:

$$C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}},{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}}_{\eta}(t)\equiv C_{\eta}(t)+\alpha(t)\Big(C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{impr.}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{impr.}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{impr.}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{impr.}$}}}}(t)-C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{\eta}(t)\Big)\,,\qquad\alpha(t)=\frac{{\rm Cov}[C_{\eta}(t),C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{\eta}(t)]}{{\rm Var}[C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{\eta}(t)]}\,.$$

(54)

Here $\alpha(t)$ is chosen to minimise the variance of $C^{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{LMA}$}}},{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}}_{\eta}(t)$ at each time slice. In this formulation, $C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{impr.}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{impr.}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{impr.}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{impr.}$}}}}(t)$ acts as a low-noise control observable, while $C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{\eta}(t)$ provides a stochastic estimator with the same expectation value. Their difference therefore has vanishing mean but non-trivial correlation with $C_{\eta}(t)$, enabling a variance reduction without introducing bias: when $C^{\mathchoice{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}{\scalebox{1.1}{$\scriptscriptstyle\mathsf{MG}$}}}_{\eta}(t)$ is strongly correlated with $C_{\eta}(t)$, a substantial variance reduction can be achieved up to the limit of perfect correlation, where the stochastic estimators cancel. In a setup where eq. 6 holds, one expects $\alpha(t)=1$. However, as shown in the right panel of fig. 7, this is not exactly true, and we find that the optimised $\alpha(t)$ deviates from unity by less than $10\%$ in the large time limit. This would lead to a further reduction in the statistical error, which is, though, only about $2\%$ and thus overall negligible.