Large-time correlation functions in bosonic lattice field theories

Large-time correlation functions have a pivotal role in extracting particle masses from Euclidean lattice field theory calculations, however little is known about the statistical properties of these quantities. In this work, the asymptotic form of the distributions of the correlation functions at vanishing momentum is determined for bosonic interacting lattice field theories with a unique gapped vacuum. It is demonstrated that the deviations from the asymptotic form at large Euclidean times can be utilized to determine the spectrum of the theory.


I. INTRODUCTION
Quantum field theory (QFT) has an essential role in nuclear and particle physics and in condensed matter physics. From the Standard Model of particle physics to topological phase structures in quantum materials, the language of QFT provides a concise and predictive mathematical description. In some cases these descriptions contain a small parameter (coupling) in which an expansion can be performed to derive analytical expressions for relevant physical quantities. However in other cases, the couplings are large and numerical approaches are required to extract physical results; the dominant such approach in many strongly-coupled theories is referred to as lattice QFT (LQFT). The LQFT method involves stochastic Monte-Carlo sampling of the very highdimensional lattice-regulated path integrals that define the correlation functions of the theory in which physical information is encoded. The large Euclidean-time behaviour of correlation functions plays a crucial role in LQFT as it is central to obtaining the energy spectrum of the theory under study. If O(x) is a localised operator built from combinations of the fundamental fields centered around the spacetime point x = (t, ⃗ x) and O(t) ∝ ⃗ x O(x) then, for Euclidean theories satisfying reflection positivity, the bilocal operator C(t) ≡ O(t)O † (0) has (vacuum, |Ω⟩) expectation value where m is the mass of the lowest energy eigenstate that O † (0) |Ω⟩ has an overlap with and Z C is the overlap factor. Therefore, for large enough 1 t, an estimator for m is given bym = ln ⟨C(t)⟩ − ln ⟨C(t + 1)⟩. (2) In practice, to extract m, one calculates the sample mean and sample variance ofm and assumes the validity of * cyunus@mit.edu † wdetmold@mit.edu 1 "Large" implies that t −1 is much smaller than the difference between the masses of the ground state and the first excited state with zero spatial momentum.
In confronting noise issues in the large-time behaviour of correlation functions, it would be useful to have an analytical form for the distribution of C(t) even in simple examples. Such a form would allow for statistical tests of empirical distributions determined in numerical calculations that may diagnose statistical limitations in the empirical sampling. In the present work, it will be demonstrated that in certain cases it is possible to obtain such analytical forms in the large time regime for correlation functions whose source and sink are both constructed to have vanishing spatial momentum. In Sec. II, the analytical form of correlation function distributions for the free real scalar field theory will be derived. In Sec. III, it will be shown that the distribution obtained in Sec. II is valid for correlation functions constructed from local operators for interacting bosonic theories at large times assuming a unique, gapped vacuum. Deviations from the asymptotic form of the distribution are linked to excited states in the spectrum and are controlled by the correlation lengths (or masses) of the theory. Consequently, analysis of these deviations provides a means of determining correlation lengths. Through their more general approach to statistical distributions, these results can potentially provide a path towards more robust determinations of energies and matrix elements in LQFT. 2  In this section, the distributions of correlation functions whose source and sink are both constructed to have vanishing spatial momentum are studied in the setting of a free real scalar field theory. The general form of the distribution is derived before a specific lattice discretisation is chosen and investigated.
Since the free theory momentum modes decouple, all non-zero spatial momentum modes can be factored out trivially. That is, if ϕ(t) is defined as where N t is the lattice size in the temporal direction, V is the spatial volume in lattice units, D(t ′ , t ′′ ) is a discretisation of the continuum Klein-Gordon operator, and [Dϕ] = Nt−1 t=0 dϕ(t) is the integration measure. For the following argument, it is assumed that the discretised D(t ′ , t ′′ ) is • positive-definite, in order that the integral in Eq.
The characteristic function for the composite operator If Q t ′ ,t ′′ = δ t ′ ,0 δ t ′′ ,t + δ t ′ ,t δ t ′′ ,0 is further defined, then it is easy to show that where Note that D and Q are linear operators acting on a vector space of dimension N t and let e t ′ for t ′ = 0, · · · , N t −1 be the basis of unit vectors on each timeslice and for which matrix elements are given as D t ′ ,t ′′ and Q t ′ ,t ′′ respectively. Note that in this basis Qe 0 = e t and Qe t = e 0 and Qe t ′ = 0 for t ′ ̸ ∈ {0, t}. Let (·, ·) be the inner product for which the vectors e t ′ are orthonormal. For the basis , t}. Therefore, to calculate det R(ω), only the subspace spanned by v 0 and v t needs to be considered. Further, the vectors v ± = D 1 2 e ± can be defined where where the ellipsis represents terms that are spanned by v j for j ̸ = 0, t and does not contribute to the determinant, then: Further, for w t ′ = D − 1 2 e t ′ and w ± = D − 1 2 e ± . Then, (w ± , v ± ) = δ ±,± . It follows that where σ, σ ′ = ±. From translational invariance, it also follows that e ± , D −1 e ∓ = e 0 , D −1 e 0 − e t , D −1 e t ± e t , D −1 e 0 − e 0 , D −1 e t = 0 and so R +− = R −+ = 0. Then, by defining 3 it follows that Therefore, the probability distribution of ϕ(t)ϕ(0) is obtained from the characteristic function as: ω-plane The calculation of P ϕ(t)ϕ(0) then reduces to the evaluation of Since ω ± > 0 (Eq. (10)), for x > 0 the integration contour in I(x) can be deformed to the contour C ≡ {i∞ − ϵ → iω + − ϵ} ∪ {iω + + ϵ → i∞ + ϵ} as in Fig. 1. The integral can thus be expressed as: whereĪ Here, K 0 (x) is a modified Bessel function of the second kind. These formulas can be summarised as 4 : The lattice action needs to be specified explicitly to gain more insight into ω ± . A simple discretisation is: . (17) Assuming that N t is odd, the normalised eigenvectors of D are given by: where j = 1, · · · , Nt−1 2 and t indexes the components of the each vector. The corresponding eigenvalues are: Expanding in eigenvectors of D, the vector representing an arbitrary zero momentum field ϕ can be expressed where ⟨. . .⟩ indicates integration over the field ϕ where the integration measure is given by Partially following Ref. [26], in order to take the continuum limit, the infinitesimal time interval ϵ is introduced and the limit t, N t → ∞ is considered such that β = ϵN t and t/N t are fixed. Further, τ ≡ tϵ is defined. The quantity m depends on ϵ (equivalently on N t ) and m(ϵ) should be chosen such that the correlation function decays as exp (−m R τ ) as ϵ → 0 for large τ and β → ∞ (the β → ∞ limit must be taken first), where m R is the renormalized mass. The renormalized field ϕ R = Z ϕ (ϵ)ϕ, where Z ϕ (ϵ) will be chosen such that correlation functions of ϕ R have non-singular behaviour as the continuum limit is taken. Setting k = 2π β j, the renormalized correlation function is: Taking the limit β → ∞ (or equivalently N t → ∞) one therefore obtains: The integration can be performed by defining z = exp (ikϵ) that maps onto the unit circle, giving 5 : with only z − occurring inside the unit circle. Therefore: This also shows that one must chose Z(ϵ) = constant × ϵ + O ϵ 2 for the above equation to be finite for ϵ → 0. Given these constraints, m(ϵ) = m R ϵ and Z ϕ (ϵ) = ϵ are chosen, leading to a renormalised correlation function that is finite in the ϵ → 0 limit.
is coincidental for d = 1 and does not hold in other numbers of spatial dimensions, d. 5 Note that τ ϵ is a non-negative integer, so the integrand is meromorphic.
It is clear from Eqs. (18) and (19) so it is straightforward to show that: The probability density function of the renormalized twopoint function is given as: and since Z(ϵ) = ϵ is chosen, this may be written equivalently as The factor ϵ leads to the modification of Eqs. (7) and (10) as and In the limit ϵ → 0: 6 Note that this is a vector equation: both D −1 et and E, S j and C j have Nt components.
the distribution in Eq. (16) is compared to the numerical distribution of ϕ(t)ϕ(0) for the two-dimensional free real scalar field theory discretised as in Eq. (17). Computations are performed for a lattice of size N s × N t = 20 × 40 for t = 6 and a sample size of N = 10 6 and the resulting distribution is shown in Fig. 2.
The parameters ω ± to use in Eq. (16) for comparison are determined as follows. The lowest moments of the distribution in Eq. (16) are given as: where ⟨x n ⟩ = ∞ −∞ dx P ϕ(t)ϕ(0) (x)x n . After some algebra, these relations can be inverted to give As well as the empirical distribution from the numerical calculations, Fig. 2 shows the distribution of Eq. (16) with the above values of ω ± determined from moments of the empirical distribution. The choice of t = 6 and lattice size are completely arbitrary and equivalently good agreement is seen for all t and for various lattice geometries. The numerical data are clearly well-represented by the expected behaviour and the empirical cumulative distribution function converges to the exact cumulative distribution function given by as the sample size is increased, as seen in Fig. 3 where sample sizes N = 10 2 and 10 6 are used.

III. LARGE-TIME CORRELATORS
The above results for a free real scalar field theory can approximately describe the statistical behaviour of correlation functions for many operators at large time in a broad class of bosonic lattice field theories. Consider a theory in D = d + 1 dimensions and a local operator 8 O(t, ⃗ x) such that: • The Euclidean action of the theory is real, • There is a unique, translationally invariant, gapped vacuum |Ω⟩, • O(t, ⃗ x) is covariant under temporal and spatial translations, • The vacuum expectation value of O(t, ⃗ x), ⟨Ω|O(t, ⃗ x)|Ω⟩, vanishes.

IfŌ(t) is defined as:
where ⃗ x runs over all lattice sites on the given time slice, then as t, N t → ∞: with ω > 0. Here, |m⟩ is the zero-momentum eigenstate with the smallest energy such that and m is its energy.
To see this, the joint probability distribution P O(t),O(0) (u, v) that O(t) takes value u and O(0) takes value v will be considered. This is defined as such thatŌ By integrating over all remaining variables, Eq. (43) can be expressed as: where here and henceforth the abbreviationŌ i ≡Ō i (0) will be used and PŌ 1 ,··· ,Ō N (u 1 , · · · , u N ) is the joint probability density of events in which eachŌ i takes the value u i and is normalised such that If the box sizes are larger than a few correlation lengths, theŌ i are independent of each other up to exponentially small effects proportional to e −ml , where m is the mass gap specified after Eq. (41) and l is the distance between centers of neighbouring boxes. Consequently, Such an approximation becomes exact in the infinite volume limit, and with N and l both taken to infinity 9 : wherePŌ i (k i ) is the Fourier transform of PŌ i (u i ). In the fourth line, that fact thatPŌ i (k) =PŌ 1 (k) by translational invariance has been used. Note thatPŌ 1 (0) = 1 as PŌ 1 (·) is a normalized probability density function. Similarly,P ′Ō 1 (0) = 0 as Ō 1 = 0 by assumption, and settingP ′′ O1 (0) = −σ 2 O1 ,PŌ 1 (λ) can be expressed as: This implies that as N → ∞: and by integrating over λ it is clear that: Therefore in the infinite-volume limit, the joint probability density PŌ (t),Ō(0) (u, v) is given as: (50) As a consequence, the distribution of the product O(t)Ō(0) in the same limit is: This expression reduces to Eq. (12) with ω ± = 1 , so finally from Eq. (16) one obtains 10 : .
(52) To test the validity of Eq. (52), interacting ϕ 4 theory in two dimensions is investigated numerically following Ref. [28]. The action for this theory is: where i labels the sites, b, u, K are couplings andμ labels the directions. Through the rescalings the action be rewritten as where ⟨ij⟩ indicates summation over all pairs of neighbouring points. A schematic illustration of the phase diagram corresponding to the above action is given in the Fig. 4 in terms of exponentials of the couplings θ and χ.
Defining the one-parameter path through the coupling space calculations have been performed for a lattice size of N s × N t = 40 × 40 and sample size N = 10 6 for s = 0.05 × k where k ∈ {2, · · · , 10} using a publicly available code [29]. For these values of s, the system is in the disordered phase, and larger s values correspond to smaller renormalised mass values, m, closer to the critical line. The parameter σ 2 O in Eq. (52) is determined through Eq. (37). To quantify how well Eq. (52) describes the numerical data, the total variance between the empirical distribution E t (q) at time t and the expected asymptotic distribution of Eq. (52) is calculated: It is expected that the total variance vanishes in the large t, β limit. The logarithm of the total variance versus t is shown in Fig. 5 where it is seen that the total variance decreases as t is increased until it reaches a plateau value that appears to be independent of s. Note that, as is seen in Eq. (52), 1 is the dominant term in the expansion of PŌ (t)Ō(0) (q) at large times and arises from the contributions of the vacuum state. The subleading terms in this expansion are due to excited states with vanishing spatial momentum. 11 It follows that: wherem is the mass of the second excited state with zero spatial momentum, ρ(µ) is the density of states, and A m and A µ are t-independent quantities. In this expression, it is assumed that 0 ≪ t ≪ β/2 so that effects of the finite temporal extent can be ignored. Considering the total variance as a function of time, it is seen that where T th (t) is the infinite sample size limit of T (t). At values of t for which A m (q) dominates the remaining contributions, one expects to find linear behaviour with slope −m. Such behaviour is seen in Fig. 5 for some values of t, but at larger t, a constant behaviour is seen. The above expansion assumes the distribution lim V →∞ PŌ (t)Ō(0) (q), while numerical calculations determine E t (q), that deviates from the exact distribution at finite statistical sampling. Since e −mt A m (q) vanishes as t increases, the deviation of E t (q) from the exact distribution will be larger than e −mt A m (q) at large times, invalidating the above expansion. The main contribution to this deviation is expected to be due to the finite sample size, since convergence to normality is generically very robust if mN s ≫ 1. Additionally, since t ∼ N t /2 in the figure, effects of the finite temporal extent may need 10 It must be stressed that the limit Nt → ∞ must be taken before the limit t → ∞ while the limit V → ∞ can be interchanged with the limits t → ∞ and Nt → ∞. 11 SinceŌ(0) is invariant under spatial translations, δ(Ō(0) − u) is also invariant. Therefore ⟨Ψ|δ(Ō(0) − u)|Ω⟩ is non-vanishing only if |Ψ⟩ has vanishing spatial momentum. This proves that PŌ (t),Ō(0) (u, v) can be expanded in eigenstates with vanishing spatial momentum and same is the true for PŌ (t)Ō(0) (q).
to be accounted for. This expectation is numerically confirmed for the system under consideration in Fig. 6 where results of calculations are shown for s = 0.5 for sample sizes N ∈ {10 4 , 10 5 , 10 6 , 10 7 }. Since Eq. (59) depends on the mass of the lowest energy state with the correct quantum numbers (Eq. (41)), the time-dependence of the total variance between the empirical distribution and the asymptotic expectation can be used to extract the corresponding mass, m. Estimates of m from the exponential falloffs seen in Fig. 5 are shown in Fig. 7 as a function of s. Estimates from fits to the time dependence of the correlation function ⟨C(t)⟩ are also shown. The fits to extract these masses are discussed in Appendix A. As can be seen in Fig. 7, the masses extracted from the total variance and the time-dependence of the correlation function itself are consistent although the total variance provides a less precise estimate for the parameters used in this study.

IV. SUMMARY AND OUTLOOK
In this work, the statistical behaviour of correlation functions of bosonic lattice field theories at large Euclidean time-separations has been investigated. The ex- act distribution of bilocal correlation functions in free scalar field theory was determined and numerical Monte-Carlo calculations were seen to converge to this distribution. It was also shown that the distributions of many correlation functions in a general class of bosonic lattice field theories approach the same universal distribution in the large-time limit. Numerical tests also confirmed this behaviour.
Extension of these results to other phenomenologically relevant theories such as lattice Quantum Chromodynamics is possible and may help in diagnosing the signalto-noise problem that plagues calculations of many quantities. A more thorough understanding of noise in Monte-Carlo sampling of lattice field theories along directions analogous to those pursued here, may lead to improved strategies for extracting physical information. In particular, tests of empirical distributions against the expected asymptotic distribution of correlation functions at large time may build confidence that a given level of sampling is sufficient for robust physical results to be determined. A deeper exploration of this direction is left to future work.
where Q q (·) is the quantile function, whose argument is a bootstrap set. A 68% confidence interval is given as Q 1 6 ({E}) − δEf , Q 5 6 ({E}) + δEf where {E} = {E f |f = 1, · · · , F } andf = argmax f W f . In Fig. 8, the calculated values of log⟨C(t)⟩(t) and T (t) are shown as a function of t for each s ∈ {0.10, 0.15, · · · , 0.50}. The resulting masses extracted from the two data sets are shown in Fig. 7 as a function of s.