Hybrid Monte Carlo Simulation with Fourier Acceleration of the N = 2 Principal Chiral Model in two Dimensions

Motivated by the similarity to QCD, specifically the property of asymptotic freedom, we simulate the dynamics of the SU(2) × SU(2) model in two dimensions using the Hybrid Monte Carlo algorithm. By introducing Fourier Acceleration, we show that critical slowing down is largely avoided and increases the simulation e ffi ciency by up to a factor of 300. This yields numerical predictions at a precision exceeding that of existing studies and allows us to verify the onset of asymptotic scaling.


Introduction
Hybrid Monte Carlo (HMC) was first introduced in [1] and developed to allow simulations of quantum chromodynamics (QCD) including fermions.However, like other Monte Carlo methods, it suffers from critical slowing down as a critical point such as a second or higher order phase transition is approached.In general, the number of steps required to obtain a statistically independent configuration, the autocorrelation time τ, scales as a power law of the correlation length ξ or the lattice spacing a. Namely, where z quantifies the severity of slowing down and typically z ∼ 2. Employing Fourier Acceleration (FA) causes the system's long and short wavelength modes to evolve at a more equal rate [2], therefore reducing the autocorrelation time significantly.Gauge invariant theories such as QCD mix the high and low frequency modes, causing the effect of FA to be masked.Nevertheless, FA has been applied successfully in SU(2) [3] and SU(3) [4] gauge theory simulations.Present attempts to bypass this issue in the HMC context include a FA inspired version of Riemannian Manifold HMC [5] and employing a softly gauge-fixed FA HMC algorithm [6].The 2-dimensional chiral model presents a semi-realistic context to illustrate the advantage of using FA in HMC as it is an asymptotically free (but non-gauge invariant) theory which develops massive states dynamically and is analytically understood [7,8,9,10].For the N = 2 case, we show that FA achieves its goal: a significant reduction of z compared to unaccelerated HMC.The paper is organised as follows.Firstly, we introduce the HMC algorithm and its extension by FA.The chiral model is discussed in Section 3. We contrast numerical results of the internal energy density against analytical expansions before computing the β function and verifying the onset of asymptotic scaling.In Section 5, the reduction in critical slowing is quantified.

The Hybrid Monte Carlo Algorithm and Fourier Acceleration
For a Euclidean action S , the HMC algorithm allows us to efficiently generate field configurations ϕ according to the target distribution exp(−S (ϕ)).This is achieved by constructing a Hamiltonian from the action and a kinetic term, featuring auxiliary momenta p conjugate to the degrees of freedom of the theory.By solving Hamilton's equations, one obtains a candidate for the Markov chain which, based on the Metropolis update [11], has a high acceptance rate while exploring distant regions of the (ϕ, p) phase space.We choose standard Gaussian momenta, such that the Hamiltonian reads where sums over lattice sites are implied by the dot notation.We solve the equations of motion (EOM) numerically using the leapfrog scheme for a fixed unit trajectory length.Subject to this constraint, and as HMC performs optimally with an acceptance rate R of 0.65 [12], the number and size of integration steps are calibrated such that R ∈ [0.6, 0.75].

Modified Dynamics
As the lattice spacing a is reduced, increasingly highfrequency modes are contained in the simulation, imposing a progressively strict upper bound on the integration step size.For a fixed number of integration steps, many sweeps are required to yield a notable change in the physical low-energy modes resulting in large autocorrelation times.One possibility to reduce the magnitude of this is to modify the dynamics.Denoting the action kernel by K, the new Hamiltonian becomes such that the modified dynamics are governed by Comparing the former (relating the time evolution of the degree of freedom to the conjugate momenta) with the relation between position and momentum in classical mechanics suggests interpreting K −1 as a mass.For a kernel involving derivatives, this mass is momentum-dependent such that the evolution of high momentum modes is slowed down while slowly evolving modes are accelerated.The overall more uniform evolution speed reduces the degree of critical slowing down.While the inversion of K in configuration space is possible in principle, the computational cost proliferates with the lattice size, motivating us to go to Fourier space where the kernel is diagonal and thus trivially inverted.Denoting the Fourier transform by a tilde, the EOM for the degrees of freedom used in practice is

Distribution of Momenta
Inserting K −1 in the kinetic term changes the distribution of the momenta and requires the sampling to be done in Fourier space by the aforementioned inversion-related problem.We will only consider real degrees of freedom such that for a D dimensional cubic lattice of side length L, the momentum Fourier modes p k are subject to Hermitian symmetry Consequently, only L D of the 2L D real and imaginary components are independent.Organising these into a real-valued object Πk such that Hence, the components Πk are sampled from independent Gaussians with zero mean and standard deviation L D / K−1 k .Reconstructing p from Π and taking the inverse Fourier transform yields the required real space momentum samples, following the modified kinetic term of Equation 3. In Appendix A we describe the construction of Πk explicitly for the D = 2 case.

The SU(2) × SU(2) Chiral Model
Let σ i be the Pauli matrices such that the SU(2) matrixvalued field ϕ(x) in the fundamental representation is parameterised by α α α(x) ∈ R 3 and can be expressed as with α = |α α α|, α α α = α α α/α and where we defined the parameter vector c µ = (c 0 , c i ) for c 0 = cos α, c i = α α α i sin α.The unitarity constraint and the det ϕ = 1 requirement impose the condition allowing us to interpret the 4 real scalar fields c µ as the degrees of freedom constraint to the S 3 group manifold.For a lattice with sites x and with basis vectors µ (both being D = 2 dimensional vectors), we choose the action of the principal chiral model as where β is defined in terms of the coupling constant g 0 via β = 4/(Ng 2 0 ).The second line motivates us to choose the action kernel as

Equations of Motion
To increase readability, we largely suppress the dependence on the Hamiltonian time t, meaning ϕ x ≡ ϕ(x, t), but for clarity write the discrete EOM using the more explicit notation.From the parameterisation in Equation 8, the momenta p x ≡ p(x, t) conjugate to the degrees of freedom α i (x, t) may be defined through the left time derivative of ϕ x : For the leapfrog integrator, the momenta are defined on a time lattice shifted by half a step size ϵ relative to the field.The discrete analogue of Equation 12is which, to linear order, can be written as the exponential update The N = 2 model allows us to evaluate the exponential exactly such that the dynamics remain on the group manifold.The finite machine precision introduces a small error which requires us to re-enforce Equation 9 occasionally.For our simulations, we found that doing so every 10 4 trajectories is sufficient.
For the Hybrid algorithm, we introduce a Gaussian kinetic term and deduce the EOM for the momenta by imposing energy conservation.Namely and, upon using Equations 10 and 12, Relabelling x → x + µ in the second term yields Based on Equation 17, we thus identify such that the momenta are updated according to

Fourier Acceleration
Using Equation 11and the Fourier representation of the Dirac delta, one quickly finds To assure that the inverse is defined for all k, the acceleration mass parameter M is introduced, yielding For free theories, taking M as the physical mass fully avoids critical slowing down.Considering an asymptotically free theory, where the modes decouple as the coupling tends to zero, motivates choosing this kernel.While the optimal value of M for interacting theories is less clear [2], the physics is independent of M, allowing us to fix M = 1/10 until Section 5.There we show that, provided 1/M is of the same order as ξ measured in units of the lattice spacing, the simulation efficiency only weakly depends on the precise value of M. For a two-dimensional lattice, the modified Hamiltonian used in practice is therefore given by Based on the general discussion in Section 2.1, modifying the dynamics changes the EOM for the degrees of freedom to The field is then updated according to

Strong and Weak Coupling Expansions
Considering a D-dimensional lattice of volume V and free energy density f for a general SU(N) × SU(N) model, we define the internal energy e as Through the standard observation of statistical physics Expansions for e in powers of the coupling for the SU(N) × SU(N) principal chiral model are discussed in detail in the literature with the results summarised in [7,13,14].After accounting for a relative factor of 4 in the definition of β, we extract the strong coupling expansion and the weak coupling expansion (30) The latter is obtained by computing loop corrections to the mean field result, resulting in the numerical constants Q 1 = 0.0958876 and Q 2 = −0.067.
The numerical result of Equation 28 and the analytic predictions of Equations 29, 30 are plotted in Figure 1 against different values of β for L = 16.Simulating 5000 trajectories for each data point was sufficient to obtain numerical results whose precision exceeds that of existing studies which use a local updating algorithm [9,15].To account for autocorrelations, we quote the standard error on the mean, corrected by the square root of the integrated autocorrelation time τ, as the statistical uncertainty δ of our measurements.Specifically, for M observations, and ρ(t) denotes the autocorrelation function while the cut K is chosen following [16].The error bars in Figure 1 have been enhanced by a factor of 10 for better readability.In their region of applicability, the analytic predictions are well matched by the numerical results.However, in the range β ∈ [0.5, 1.0] the expansions lose applicability and the lattice simulation offers corrections of up to 6%.

Correlation Function
Due to the taken periodic boundary conditions, the two-point correlation function of ϕ with correlation length ξ is symmetric and given by The computation of the correlation function becomes expensive for large lattice sizes, motivating us to employ an efficient approach using Fourier transforms described in Appendix B and based on the close relationship between correlations and convolutions.Additionally, we normalise the correlation to its value at d = 0 such that any constants in Equation 33 may be ignored.Again, we compute the error of the correlation function as the autocorrelation time corrected standard error on the mean.Figure 3 shows the numerical results for β = 0.8667 on a L = 64 lattice, the fitted Equation 32 with the associated reduced χ 2 , and the inferred correlation length.

Asymptotic Scaling and the β Function
Asymptotic scaling requires the dimensionless ratio of any dimensionful quantity to the appropriate power of the regular- isation scale Λ to approach a constant as the coupling strength tends to zero.The β function is known at three-loop accuracy using lattice regularisation with where G 1 = 0.04616363 [7,8].Integrating Equation 34 and defining the inverse correlation length ξ, measured in units of the lattice spacing, as the mass m yields The remaining part of this section is dedicated to testing the stability of Equation 36 with β for the N = 2 case against the mass gap scaling prediction computed in [10].In order to do so, we evaluate the integral in Equation 36 numerically using the three-loop β function and compute the correlation length ξ for various values of β based on 10 5 trajectories (rejecting the first 2000 as burn in) using FA HMC with a fixed acceleration mass of 1/10.The results are presented in Table 1.For each value of β, the lattice size L was chosen such that L/ξ ≳ 10 to reduce systematic errors due to finite size effects as suggested by [7].The fitting range for  Equation 32 was determined manually to yield χ 2 r ≈ 1 while excluding the noise-dominated regions of the correlation function data.Towards the lower end of the surveyed range of β, the data becomes quickly dominated by noise, such that only a small number of data points are suitable for the fitting, resulting in χ 2 r smaller than 1.The collected data of Equation 36 is plotted against β in Figure 4 with the continuum limit prediction of Equation 37 indicated by the dashed line.The overall convergence to this constant is fairly poor, overshooting it by up to 12% at intermediate values of β.This is consistent with the results of the N = 3, 6, 9, 15 cases and believed to be caused by a dip in the lattice beta function [7,8].The behaviour near the upper bound of the considered β range suggests the onset of scaling.

Critical Slowing Down Elimination
As suggested in [17], using the susceptibility χ as the observable to measure the effect of Fourier Acceleration is beneficial as it is dominated by low momentum modes and thus particularly prone to critical slowing down.Using any other variable should, however, result in the same conclusions but possibly with a less dramatic reduction in z.For a lattice of volume L 2 , the susceptibility of a single configuration is defined as and can be viewed as the average point-to-point correlation of that configuration.In practice, performing the double sum over lattice sites quickly becomes the bottleneck for the analysis as L gets large.It is therefore vital to compute χ via FFTs and the computationally much more efficient form given in Equation B.4.

Calibrating the Acceleration Mass
To optimise the performance of FA HMC through the choice of M, we conducted a grid search at fixed β = 1.1333 based on the cost function cost(M) = simulation time acceptance rate as the measurement error of an observable scales as √ τ.
Figure 5 details our conclusion that the conventional choice M = 1/ξ yields close to optimal acceleration, being roughly 10% slower than the extremum of the scanned values of M.
We adopt this choice henceforth and show in the following that this still offers a substantial speed-up compared to unaccelerated HMC.

Results
For each value pair β, L from Table 1 we simulated 10 5 trajectories (5000 burn in) using both standard HMC and FA HMC to deduce the susceptibility integrated autocorrelation time τ χ in either case.Due to terminating the sum in Equation 31, the error of τ χ scales as τ χ /M, such that for a fixed number of observations M, the uncertainty in τ χ increases towards the continuum.The data, as well as the fitted power law of Equation 1, are plotted in Figure 6.
For FA HMC, we inferred the dynamical critical exponent z to be 0.21 ± 0.01, offering a significant reduction compared to standard HMC with z = 1.777 ± 0.037.Even though critical The integrated autocorrelation time of the susceptibility against the system correlation length ξ on log-log axes for standard HMC (blue) and FA HMC (red) simulations.Fitted power laws are superimposed with the value of the critical exponent presented in the legend.Each data point corresponds to 10 5 trajectories (5000 burn in) with values of L and β as specified in Table 1.slowing down was not fully eliminated, employing Fourier Acceleration allows us to gain a factor of roughly 320 in the susceptibility autocorrelation time at β = 1.4,ξ ≈ 91.25.Consequently, one can expect that HMC requires 300 times more trajectories than FA HMC to achieve comparable accuracy, translating to significantly larger CPU time requirements.In terms of the cost function, this corresponds to a factor of around 25 as shown by plotting their ratio in Figure 7.We find the cost advantage scales roughly as ξ 0.78 .Indeed, as the functional dependency of the cost function is dominated by √ τ and as the simulation time and acceptance rate are approximately the same for both algorithm choices, the same scaling is also found from half the difference of the gradients in Figure 6.

Conclusions
The main conclusion from this study is the significant reduction of the severity of critical slowing down when introducing Fourier Acceleration to the Hybrid Monte Carlo simulations of the principal chiral model.The value of the acceleration mass parameter was shown to be not critical and the natural choice of setting it equal to the inverse of the correlation length yields close to optimal acceleration.The simplicity of this and the parallels between the principal chiral model and physical gauge theories such as QCD further encourages investigating Fourier Acceleration in the context of the latter.In addition, we computed the mass over Λ ratio to confirm the asymptotic scaling prediction and find our results in agreement with other SU(N) × SU(N) studies.To achieve these objectives, we developed a simulation and analysis package which can serve as a foundation for further studies of SU(N) × SU(N) principal chiral models and is available in [18].It will be interesting to see how the performance of FA HMC compared to standard HMC develops against N.Moreover, we implemented an efficient approach to compute (auto-)correlation functions, resulting in a negligible time requirement for data analysis.While in this article for simplicity we have tried to keep the acceptance approximately constant to be able to compare different runs, a further interesting avenue would be to investigate the influence of the mass parameter on the acceptance.achieves this.The factors of √ 2 are required as due to Hermitian symmetry Hence Equation 7 follows.

Appendix B. Correlations via Fourier Transforms
Consider the two-point correlation function of a generic observable O with observations O m and average ⟨O⟩.Defining f m ≡ O m − ⟨O⟩ and denoting the discrete convolution by * shows that correlations can be written as a convolution in the observation chain index m: The cross-convolution theorem then yields which can be computed efficiently using FFTs.In this study, the typical number of elements to correlate is O 10 5 for which we measured the FFT approach to be six orders of magnitude faster than the naive summation.However, care must be taken when using FFTs as their product yields the circular convolution of the original signals rather than the linear one.For physical correlations, the former automatically accounts for the periodic lattice boundary conditions and is therefore desirable.On the other hand, for autocorrelations in the chain of measurements, a linear convolution is required.This can be obtained by simply padding the original data with sufficiently many zeros such that the circular mixing with these leaves the original data unchanged.The padding is usually done in powers of two as FFT routines work most efficiently with such a number of elements.Once all Fourier transforms have been taken, one must truncate the introduced pads to get the final autocovariance function result.
To employ this approach for the computation of the wall-to-wall correlation in Equation 33, two observations are required.with Φ m (i) = p ϕ m (v, i) and where the summation over α = 0, 1, 2, 3 in Euclidean space is implied.Note that for each value of α, the remaining sum is a convolution in the variable i of the scalar quantity Φ m and is thus efficiently computed using FFTs.Similarly, we write Equation 38 as Due to the lattice periodicity and since i and j both index columns of the lattice, one can always express i as a shifted version of j such that the inner bracket is simply the convolution already encountered in Equation B.3.

Figure 1 :
Figure 1: The internal energy density against β based on simulating 5000 trajectories (10% burn in and L = 16) per data point using Fourier accelerated HMC with a fixed acceleration mass of M = 1/10.Error bars are enlarged by a factor of 10.The strong (blue) and weak (red) coupling expansions are superimposed.

Figure 2 :
Figure 2: Composition of the wall-to-wall correlation function at separation d = 3 in terms of point-to-point correlations.

)
To obtain ξ, we exploit this symmetry by fitting C(d) to the average correlation function data at d and L − d on d ∈ [0, L/2], effectively doubling the statistics.We compute the correlation function based on wall-to-wall correlations which are defined as the average point-to-point correlations contained by two walls in the lattice with separation d.This is schematically shown in Figure2and has the benefit that the number of wall pairs P in each field configuration is large, increasing the statistics for C(d).We omit the Hamiltonian time dependence and expand the lattice site x into its coordinates (v, i).The mth field configuration evaluated at the lattice site (v, i) is thus denoted by ϕ m (v, i).At a fixed value of i, averaging over the row indices v and w of the point pairs in Figure 2 yields one measurement of the wall-to-wall correlation at separation d for the currently considered configuration.Denoting the ensemble average as ⟨. . .⟩ m , the wall-to-wall correlation becomes

Figure 3 :
Figure3: The wall-to-wall correlation function inferred from 10 5 trajectories (5% burn in) at β = 0.8667 on a L = 64 lattice and plotted on a logarithmic y scale.A cosh function is fitted with the inferred correlation length and the associated χ 2 per degree of freedom presented in the legend.

Figure 4 :
Figure 4: The mass over Λ ratio of Equation 36 against β.Each data point is the result of fitting the correlation length to the averaged correlation function data based on 10 5 trajectories (2000 burn in) using the lattice size L specified in Table. 1.The continuum prediction of Equation 37 is shown as a dashed line.

Figure 5 :
Figure 5: The cost function of Equation 39 for different values of the acceleration mass, normalised to the M = 1/ξ choice.Each data point corresponds to 1.5 × 10 4 trajectories (with 10% burn in) for L = 160, β = 1.1333.

χFigure 6 :
Figure6: The integrated autocorrelation time of the susceptibility against the system correlation length ξ on log-log axes for standard HMC (blue) and FA HMC (red) simulations.Fitted power laws are superimposed with the value of the critical exponent presented in the legend.Each data point corresponds to 10 5 trajectories (5000 burn in) with values of L and β as specified in Table1.

Figure 7 :
Figure7: Ratio of the standard HMC and FA HMC cost functions against the system correlation length on a log-log scale with a fitted power law superimposed.The simulation parameters are identical to those used in Figure6.

Table 1 :
The fitted correlation length, its error, and the associated χ 2 r for different pairings of β and L. The lattice size was chosen conservatively to minimise finite size effects.All simulations employ the FA HMC algorithm with the acceleration mass fixed to 1/10 and run for 10 5 trajectories, rejecting the first 2000 as burn in.