Origin of the sharp dependence of Parisi's fluctuation-dissipation noise-reduction method on the spin-component number in the two-dimensional O($n$) sigma model

In lattice simulations, the connected two-point correlation function of an operator $g$ with an operator $f$ can be estimated from the linear response of $\langle g\rangle$ to a small source $h$ coupled to $f$ in the action, via the fluctuation-dissipation relation $\langle g f\rangle_c = \lim_{h\to 0}(\langle g\rangle_h - \langle g\rangle_0)/h$ (up to the overall sign fixed by the convention for coupling $h$). Parisi (Nucl. Phys. B 205, 337 (1982)) showed that, in the two-dimensional O($n$) nonlinear sigma model simulated by a discrete Langevin equation, this response ('subtraction') estimator reduces the statistical error on the long-distance spin-spin correlation by a factor of order $10^6$ at equal computer time for $n=4$, while for $n=3$ it gives no improvement over direct averaging; the method is exact and optimal in the free-field limit $n\to\infty$. He remarked that an analytic study of this surprisingly sharp dependence on the number of spin components 'would be welcome'. Forty years later, the surrounding agenda of the source paper - efficient computation of plaquette-plaquette (glueball) correlators in pure gauge theory, stochastic (Langevin) quantization of SU(3), momentum-dependent one-particle-irreducible functions, and per-configuration propagator-matrix estimators - has been largely realized or superseded by multilevel sampling, adjoint-source methods, and modern distillation technology, but no analytic or controlled numerical account of the $n$-dependence of the response estimator's efficiency exists. The open problem is the analytic study Parisi requested: to explain quantitatively the sharp $n$-dependence of the response estimator's efficiency in the 2D O($n$) sigma model - the reported dichotomy of no gain at $n=3$ versus a gain of order $10^6$ at $n=4$ - and to recover the free-field optimality at $n\to\infty$, while engaging the recorded reliability caveats on the original measurement.

Background

Monte Carlo simulation is the standard nonperturbative tool for lattice-regulated quantum and statistical field theories: one generates field configurations $\Phi$ distributed according to the Boltzmann weight $e^{-A(\Phi)}$ for a Euclidean action $A$ and estimates observables by ensemble averages. The quantities of most physical interest - masses, mass gaps, correlation lengths - are extracted from two-point correlation functions, which decay exponentially with separation $|x|$ as $e^{-m|x|}$, while the variance of a direct product estimator $g(0)g(x)$ is governed by the lightest states that can propagate between the two insertions. For correlators of composite operators the variance typically decays much more slowly than the signal (the 'signal-to-noise problem', articulated for lattice QCD by Parisi, Phys. Rep. 103 (1984) 203), so the relative statistical error grows exponentially with separation and the computer time needed for fixed precision grows as the square of the required precision.

The two-dimensional O($n$) nonlinear sigma model is a standard laboratory for lattice methodology: each site carries a unit spin $\sigma_i \in S^{n-1}$ ($i=1,\dots,n$ spin components) with nearest-neighbor coupling $-\beta^{-1}\sum_{\langle xy\rangle}\sigma(x)\cdot\sigma(y)$; for $n\ge 3$ the model is asymptotically free in two dimensions with a dynamically generated mass gap $m$, and in the formal limit $n\to\infty$ it becomes a free theory of $n-1$ massive components. Parisi simulated it by stochastic quantization: instead of a Markov chain over configurations, one evolves a fictitious Langevin time $t$ according to a stochastic differential equation whose stationary ensemble is $e^{-A}$, so that long Langevin-time averages replace ensemble averages; the unit-length constraint on the spins is enforced by a Brownian motion on the sphere (projecting the drift and adding the appropriate curvature term to the noise, whose per-step variance is $(2+c\epsilon)/\epsilon$ with step size $\epsilon$ and constant $c=2$ in the source's convention).

The response (or 'subtraction') estimator exploits the fluctuation-dissipation relation: perturb the action to $A(\Phi) + h\, f(\Phi)$; then, with the Boltzmann weight $e^{-A}$, the connected correlator is $\langle g f\rangle_c = \lim_{h\to 0}\big(\langle g\rangle_0 - \langle g\rangle_h\big)/h$ (the source writes this relation with the opposite overall sign convention, which is irrelevant for the error analysis), where each expectation value is a Langevin-time average. This is a control-variate-type estimator: the source-induced shift of $\langle g\rangle$ may have far smaller fluctuations than the product $g f$ itself. A precursor for systems with fermions is the method of Scalapino and Sugar (Phys. Rev. Lett. 46 (1981) 519), in which correlation functions are obtained from the response to an external source coupled to a bilinear. Parisi's companion paper (Nucl. Phys. B 180 (1981) 378) proved optimality of the estimator in the free-field case, and the 1982 paper implemented it for the 2D O($n$) model, measuring the correlation of the column-summed spin $S(l)=\sqrt{n/N}\sum_r \sigma^1(l,r)$ (summing the first spin component over a time-slice column of width $N$, which increases the overlap with the single-particle state).

The empirical result and the puzzle: on a $30\times 30$ lattice at $\beta=0.4$ with $n=4$, source strength $h=10^{-4}$ and relaxation time $\tau=100$ Langevin steps per measurement, the estimator reached statistical errors of order $10^{-9}$ on the last visible point of the correlation function, where an equal-CPU direct Monte Carlo average gives errors of order $10^{-3}$ - a gain of about $10^6$ in precision, or $10^{12}$ in computer time, consistent with the CDC 7600 rounding limit of $10^{-10}$. Similar gains were found for other lattice sizes, other $\beta$, and larger $n$ (values $n=4,5,8$ on lattices from $6\times 6$ to $30\times 30$), and the method is exact in the free limit $n\to\infty$. But for $n=3$ there was no improvement over the direct average. Parisi wrote: 'It is not easy to understand why the method does not work for n=3... It is however surprising that going from n=3 to n=4 one gains a factor 10^6; an analytic study of this effect would be welcome.'

The source paper's broader agenda concerned gauge theories. For composite operators it found the response method only 'slightly better' than direct averaging (for the energy-energy correlator the estimator is contaminated by terms $e^{-2m|x|} + \langle\sigma\rangle^2 e^{-m|x|}$ from the finite-run magnetization), and for pure gauge theories the plaquette-plaquette correlator (the elementary square Wilson loop, whose correlation measures glueball masses - bound states of the gauge field in a theory without matter) would be swamped by slowly decaying non-gauge-invariant link-link contributions. Parisi also noted that computing one-particle-irreducible (1PI) correlators - those immune to such disconnected contaminations, obtained formally by a Legendre transform of the generating functional - is 'not easy to do unless we stay at zero external momentum', and proposed, for $\Phi^4$-type actions reducible to quadratic form via a Symanzik/Hubbard-Stratonovich transformation (an exact auxiliary-field rewriting that trades the quartic self-interaction for a quadratic form coupled to an auxiliary field) with operator $D=-\Delta+m^2+g\sum_i\Phi_i^2+(-g)^{1/2}\sigma$, a per-configuration evaluation of matrix elements of $D^{-1}$ as a compositeness-free estimator. He closed: 'No very efficient method has been found in this paper for computing the plaquette-plaquette correlation function. I hope that such a method exists and that it will be found by someone else.'

That wish was granted by a different route. Variational smearing and multihit/link-averaging (Parisi-Petronzio-Rapuano, Phys. Lett. B 128 (1983) 418) cut the variance of Wilson-loop measurements; stochastic quantization of non-abelian gauge theories was systematically developed in the 1980s (review: Damgaard and Hüffel, Phys. Rep. 152 (1987) 227), making Langevin dynamics on SU(3) link variables routine. Above all, multilevel (two-level) Monte Carlo exploiting the locality of the action (Lüscher and Weisz, JHEP 09 (2001) 010) achieves an exponential reduction of statistical errors for loop correlation functions in pure SU(3) gauge theory: the four-dimensional lattice is sliced, an ensemble of sub-lattice configurations is generated independently in each slice, and the correlator factorizes into products of sub-lattice averages that can be recombined at no extra error. Meyer (JHEP 01 (2003) 048) gave a general theory of multilevel performance tied to the low-energy spectrum and measured Wilson-loop pair correlators with constant relative errors; Cè, Giusti and Schaefer (Phys. Rev. D 93 (2016) 094507) extended the scheme to theories with fermions via domain decomposition; Barca et al. (Phys. Rev. D 110 (2024) 054515) obtained exponential error reduction for glueball correlators in pure-gauge SU(3) (updated for quenched QCD in PoS LATTICE2024 (2025) 062), while Kitching-Morley and Jüttner (PoS LATTICE2021 (2022) 133) analyzed the saturation of the multilevel gain near criticality. Momentum-dependent matrix elements are now obtained from adjoint sources (Michael, Nucl. Phys. B 259 (1985) 58) and direct three-point-function computations at general momenta (e.g. Skullerud and Kizilersü, JHEP 09 (2002) 013), and per-configuration propagator-matrix estimators are standard through distillation (a framework restricting the quark fields to a fixed low-rank subspace of each configuration, from which all needed propagator data are obtained as small per-configuration matrices; Phys. Rev. D 80 (2009) 054506). What no later line picked up is the one question Parisi explicitly flagged as lacking an analytic study: the sharp $n$-dependence of the response estimator itself.

Problem Statement

For the two-dimensional O($n$) nonlinear sigma model (unit spins $\sigma_i \in S^{n-1}$, nearest-neighbor coupling, inverse coupling $\beta$) on finite square lattices, consider the connected correlation function of the column-summed spin $S(l)=\sqrt{n/N}\sum_r \sigma^1(l,r)$, estimated by the response (subtraction) estimator $\lim_{h\to 0}\big(\langle S(l)\rangle_h-\langle S(l)\rangle_0\big)/h$ (equal to the connected correlator up to the overall sign fixed by the source-coupling convention), where a source $h\,S(1)$ is added to the action and expectation values are measured as long-time averages of a correct stochastic dynamics (e.g. the projected discrete Langevin equation of the source, or any Markov chain with the same stationary ensemble). The reported facts are: at equal computer budget the estimator's statistical error is smaller than that of the direct product estimator by a factor of order $10^6$ for $n=4$ on a $30\times 30$ lattice at $\beta=0.4$ (with similar gains for $n=5,8$ and other $\beta$ and lattice sizes), while for $n=3$ there is no gain; in the formal limit $n\to\infty$ the theory is free and the estimator is known to be optimal. The open problem, as posed by the source, is the analytic study of this surprisingly sharp dependence: provide an analytic, quantitative account of the estimator's variance (or of the error-gain factor) as a function of $n$, comparing different $n$ at fixed inverse coupling $\beta$ and fixed lattice size (a comparison regime of the solver's own choosing; the reported gains span several $\beta$ and lattice sizes and the conditions of the $n=3$ null are not recorded in the source). The account is to establish: (i) the $n$-dependent variance or error-gain factor; (ii) the origin of the reported $n=3$ versus $n=4$ dichotomy - either a physical mechanism whose effect changes sharply between $n=3$ and $n=4$, or a smooth $n$-dependence, including the possibility that the $n=3$ null is an artifact of the recorded measurement caveats; and (iii) recovery of the free-field optimality at $n\to\infty$ as a consistency requirement. A controlled numerical study of the estimator's efficiency may support the account. The account must engage the recorded reliability caveats on the original measurement (single finite source strength without $h\to 0$ extrapolation, undocumented error estimation, and estimator-variant dependence). The scope is the response/subtraction estimator in the form reported by the source - the uncorrelated difference of $\langle S(l)\rangle_h$ and $\langle S(l)\rangle_0$ under a single source $h\,S(1)$ - for two-point spin correlation functions of the 2D O($n$) sigma model; no restriction on $n$, $\beta$, or lattice size beyond what is stated is implied.

References

1. G. Parisi, 'Correlation functions and computer simulations (II)', Nuclear Physics B 205 (1982) 337-344, DOI 10.1016/0550-3213(82)90360-1, https://doi.org/10.1016/0550-3213(82)90360-1
2. G. Parisi, 'Correlation functions and computer simulations', Nuclear Physics B 180 (1981) 378-384, DOI 10.1016/0550-3213(81)90056-0, https://doi.org/10.1016/0550-3213(81)90056-0
3. D. J. Scalapino and R. L. Sugar, 'Method for performing Monte Carlo calculations for systems with fermions', Physical Review Letters 46 (1981) 519-521, DOI 10.1103/PhysRevLett.46.519, https://doi.org/10.1103/PhysRevLett.46.519
4. G. Parisi, R. Petronzio and F. Rapuano, 'A measurement of the string tension near the continuum limit', Physics Letters B 128 (1983) 418-420, DOI 10.1016/0370-2693(83)90930-9, https://doi.org/10.1016/0370-2693(83)90930-9
5. G. Parisi, 'The strategy for computing the hadronic mass spectrum', Physics Reports 103 (1984) 203-211, DOI 10.1016/0370-1573(84)90081-4, https://doi.org/10.1016/0370-1573(84)90081-4
6. O. Martin, 'Correlated subtractions for computer simulations', Nuclear Physics B 251 (1985) 425-438, DOI 10.1016/0550-3213(85)90270-6, https://doi.org/10.1016/0550-3213(85)90270-6
7. P. de Forcrand and C. Roiesnel, 'Refined methods for measuring large-distance correlations', Physics Letters B 151 (1985) 77-80, DOI 10.1016/0370-2693(85)90826-3, https://doi.org/10.1016/0370-2693(85)90826-3
8. C. Michael, 'Adjoint sources in lattice gauge theory', Nuclear Physics B 259 (1985) 58-76, DOI 10.1016/0550-3213(85)90297-4, https://doi.org/10.1016/0550-3213(85)90297-4
9. M. Namiki, I. Ohba, K. Okano, M. Rikihisa and S. Tanaka, 'New procedure of numerical simulation of long-range correlations and energy gaps by stochastic quantization: fixed potential model and O(3) nonlinear sigma-model', Progress of Theoretical Physics 76 (1986) 501-511, DOI 10.1143/PTP.76.501, https://doi.org/10.1143/PTP.76.501
10. P. H. Damgaard and H. Hüffel, 'Stochastic quantization', Physics Reports 152 (1987) 227-398, DOI 10.1016/0370-1573(87)90144-X, https://doi.org/10.1016/0370-1573(87)90144-X
11. E. Abdalla, M. C. B. Abdalla, N. A. Alves and C. E. I. Carneiro, 'Numerical simulations of the O(3) and CP$^1$ models using the Langevin equations and the Metropolis algorithm', Physical Review D 41 (1990) 571-580, DOI 10.1103/PhysRevD.41.571, https://doi.org/10.1103/PhysRevD.41.571
12. J. I. Skullerud and A. Kizilersü, 'Quark-gluon vertex from lattice QCD', Journal of High Energy Physics 0209 (2002) 013, DOI 10.1088/1126-6708/2002/09/013, https://doi.org/10.1088/1126-6708/2002/09/013
13. M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Joó, H.-W. Lin, D. G. Richards and K. J. Juge, 'Novel quark-field creation operator construction for hadronic physics in lattice QCD' (the distillation method), Physical Review D 80 (2009) 054506, DOI 10.1103/PhysRevD.80.054506, https://doi.org/10.1103/PhysRevD.80.054506
14. M. Lüscher and P. Weisz, 'Locality and exponential error reduction in numerical lattice gauge theory', Journal of High Energy Physics 0109 (2001) 010, arXiv:hep-lat/0108014, DOI 10.1088/1126-6708/2001/09/010, https://doi.org/10.1088/1126-6708/2001/09/010
15. H. B. Meyer, 'Locality and statistical error reduction on correlation functions', Journal of High Energy Physics 0301 (2003) 048, arXiv:hep-lat/0209145, DOI 10.1088/1126-6708/2003/01/048, https://doi.org/10.1088/1126-6708/2003/01/048
16. M. Cè, L. Giusti and S. Schaefer, 'Domain decomposition, multilevel integration, and exponential noise reduction in lattice QCD', Physical Review D 93 (2016) 094507, arXiv:1601.04587, DOI 10.1103/PhysRevD.93.094507, https://doi.org/10.1103/PhysRevD.93.094507
17. L. Altenkort, A. M. Eller, O. Kaczmarek, L. Mazur, G. D. Moore and H. T. Shu, 'Lattice QCD noise reduction for bosonic correlators through blocking', Physical Review D 105 (2022) 094505, arXiv:2112.02282, DOI 10.1103/PhysRevD.105.094505, https://doi.org/10.1103/PhysRevD.105.094505
18. B. Kitching-Morley and A. Jüttner, 'A numerical and theoretical study of multilevel performance for two-point correlator calculations', Proceedings of Science LATTICE2021 (2022) 133, arXiv:2201.01610, DOI 10.22323/1.396.0133, https://doi.org/10.22323/1.396.0133
19. L. Barca, S. Schaefer, F. Knechtli, J. A. Urrea-Niño, S. Martins and M. Peardon, 'Exponential error reduction for glueball calculations using a two-level algorithm in pure gauge theory', Physical Review D 110 (2024) 054515, arXiv:2406.12656, DOI 10.1103/PhysRevD.110.054515, https://doi.org/10.1103/PhysRevD.110.054515
20. L. Barca, J. Finkenrath, F. Knechtli, M. J. Peardon, S. Schaefer and J. A. Urrea-Nino, 'Update on two-level sampling for glueball observables in quenched QCD', Proceedings of Science LATTICE2024 (2025) 062, arXiv:2501.17988, DOI 10.22323/1.466.0062, https://doi.org/10.22323/1.466.0062
21. T. Harris, 'Variance-reduction strategies for lattice QCD', Proceedings of Science LATTICE2025 (2026) 009, arXiv:2605.00643, DOI 10.22323/1.518.0009, https://doi.org/10.22323/1.518.0009