Universal Dynamics of Independent Critical Relaxation Modes

Scaling behavior is studied of several dominant eigenvalues of spectra of Markov matrices and the associated correlation times governing critical slowing down in models in the universality class of the two-dimensional Ising model. A scheme is developed to optimize variational approximants of progressively rapid, independent relaxation modes. These approximants are used to reduce the variance of results obtained by means of an adaptation of a quantum Monte Carlo method to compute eigenvalues subject to errors predominantly of statistical nature. The resulting spectra and correlation times are found to be universal up to a single, non-universal time scale for each model.

One of the most remarkable characteristics of critical behavior is universality. For instance, it is generally accepted that upon approach of a critical point, the correlation length diverges with a power law and that the exponents are universal in the sense that they depend only on the qualitative features characterizing the direction of approach of the critical point. The definition of the correlation length can be based on any of various correlation functions, the most obvious ones of which are the order parameter and energy autocorrelation functions. The nature of the definition shows up in the amplitude of the power law: Within a given universality class, this amplitude is of the form mA, where A is universal, but depends on the observables used in the definition of the correlation length; m is non-universal, yet independent of these detail, and a function only of the representative of the universality class.
In terms of the spectrum of the transfer matrix one can define an infinity hierarchy of correlation lengths. Indeed, the spectrum of the transfer matrices, which has been studied extensively for two-dimensional systems, displays such universal amplitudes. 1 For statics the transfer matrix generates translations in space. The Markov matrix is its analog for translations in time in for stochastic dynamics. From this perspective, the analog of the correlation length is the correlation time and in this Letter we address universality of the critical point amplitudes that describe the vanishing of gaps of the spectrum of the Markov matrix and the slowing-down of the independent modes.
The dynamic critical exponent z links the divergences of the correlation length and time, ξ and time τ , via the relation τ ∼ ξ z . For a finite system of linear dimension L, with its thermodynamic fields fixed at the critical values of the infinite system, this becomes τ ∼ L z , since the correlation length is limited by the size of the system. We denote by 1 = λ L0 > λ L1 ≥ . . . the eigenvalues of the Markov matrix, and by ∆ Li = 1 − λ Li the spectral gaps. For the associated correlation times one has τ −1 Li = L d ln λ Li and the latter are expected to display the following scaling for L → ∞. Here A i and z are universal, and m t is a non-universal, metric factor. [2][3][4] Furthermore, d is the dimensionality of the system, which enters because we consider a Markov matrix that evolves the system in a local sense only. More specifically, we study single-spin-flip dynamics in Ising models defined on a square lattice of linear dimension L with nearest and next-nearest neighbor couplings K and K ′ and periodic boundaries. We focus on models described by three ratios β = K ′ /K, namely β = − 1 4 , 0, 1, the opposite-, nearest-, and equivalent-neighbor models.
The dynamics is generated by the Markov matrix P , element P (S ′ , S) of which defines, given a configuration S, the conditional probability of a transition to a configuration S ′ . If S and S ′ differ by more than one spin, P (S ′ , S) = 0. If both configurations differ by precisely one spin, where H denotes the spin Hamiltonian. The diagonal elements P (S, S) follow from the conservation of probability S ′ P (S ′ , S) = 1.
We compute the spectrum of P by means of a method used previously for a single eigenstate 5 generalized to several dominant eigenvalues of the Markov matrix. This method was introduced by Ceperley and Bernu in the context of quantum Monte Carlo methods 6 . A crucial element in our approach is the construction of optimized trial states, and for this purpose we generalize ideas of Umrigar et al. 7 .
The condition of detailed balance is used to define a stochastic process that has the Boltzmann distribution exp[−H(S)/kT ] ≡ ψ B (S) 2 as its stationary state. Consequently, P (S ′ , S) ≡ 1 ψ B (S ′ ) P (S ′ , S)ψ B (S) is symmetric in S and S ′ . We write eigenvectors of the transformP in the functional form ψ (±) (S)ψ B (S) defined in Eqs. (12) and (13) of Ref. 5. That is, in the first place we restrict ourselves to translationally and rotationally modes, which are even or odd under spin inversion, and, secondly, the trial functions are written as linear combinations of monomials of the magnetization and other long-wavelength Fourier transforms of the spin configuration.
We generalize to simultaneous optimization of multiple trial states, a powerful method 7-9 of optimizing a single many-parameter trial function. The latter is done by minimization of the variance of the configurational eigenvalue: Suppose that ψ T (S, p) is the value of the trial function ψ T for configuration S and some choice of the parameters p to be optimized. The configurational eigenvalue λ(S, p) of a spin configuration S is defined by where the prime indicates matrix multiplication byP , i.e., f ′ (S) ≡ S ′P (S, S ′ )f (S ′ ) for arbitrary f . The optimal values of the variational parameters are obtained by minimization of the variance of λ(S, p), estimated by means of a Monte Carlo sample. We refer to Ref. 5 for details and mention only one of the key features of this method: in the ideal case, i.e., for an exact eigenstate ψ T , the variance vanishes if it were to be computed exactly but also if one employs an approximate Monte Carlo expression. A similar zero-variance principle holds for the method of simultaneous optimization of several trial states to be discussed next.
For simplicity of presentation we first generalize the above method to the more general ideal case in which one can exactly compute m eigenvalues of the Markov matrixP . Suppose we have m basis states ψ Ti , i = 1, . . . , m and again M spin configurations S α , α = 1, . . . , M sampled from ψ B 2 . The case we consider is ideal in the sense that we assume that these states ψ Ti span an m-dimensional invariant subspace ofP . In that case, by definition there exists a matrixΛ of order m such that Again, the prime on the left-hand side of this equation indicates matrix multiplication bỹ P . If M is large enough,Λ is for all practical purposes determined uniquely by the set of equations (4) and one findsΛ and where Z is an arbitrary normalization constant; again, the prime indicates matrix multiplication byP . In the non-ideal case, the space spanned by the m basis states ψ Ti is not an invariant subspace of the matrixP . In that case, even though Eq. (4) generically has no true solution, Eqs. (5) and (6) still constitute a solution in the least-squares sense, as may be verified by solving the normal equations. If a set of states span an invariant subspace, so does any non-singular set of linear combinations. In principle, the optimization criterion should have the same invariance. The spectrum of the matrixΛ has this property, which suggests that one subdivide the sample in subsamples and compute the variance of the local spectrum over these subsamples. In practice, however, precisely this invariance gives rise to a near-singular non-linear optimization problem. Therefore, to avoid slow or no convergence, we add a contribution to the above least-squares merit function to ensure that the basis states themselves are good approximate eigenstates, rather than just their linear combinations, and we use an iterative optimization procedure: First a combination of the single and multi-eigenstate merit functions is used, and finally the resulting approximate eigenstates are optimized one at a time using the single-state procedure only. Unfortunately, this method is capricious and often we proceed by trial and error.
The variational states can be used directly only to obtain results with systematic errors, but these can be suppressed by the quantum Monte Carlo projection method introduced by Ceperley and Bernu 6 . Define generalized matrix elements For t = 0, Eqs. (6) are Monte Carlo estimators for these matrix elements, apart from the inconsequential normalization constant Z. One can view the matrix elements for t > 0 as having been obtained by the substitution |ψ Ti →P t/2 |ψ Ti , which implies that spectral weights of "undesirable" states lower down in the spectrum are reduced. The matrix elements in Eqs. (7) are the following time auto-and cross-correlation functions of the Markov process generated by the matrix P : ψ Ti (S 0 )ψ Tj (S t ) ψ B 2 and ψ Ti (S 0 )ψ ′ Tj (S t ) ψ B 2 , where S 0 and S t are spin configurations separated in time by t single-spin flips.
It should be noted that in the limit of vanishing statistical error, each eigenvalue estimate obtained by the above method is bounded from above by the corresponding exact eigenvalue. The reader is referred to the work of Ceperley and Bernu in Ref. 6 for further details and references. The systematic error decreases for increasing projection time t while the statistical error increases. An optimal intermediate t has to be chosen, which yields biased estimators and some uncertainty in the reliability of statistical error estimates.
Of the three Ising-like models investigated here, the critical point is exactly known only for the nearest-neighbor model, where it occurs at K = K c (0) = 1 2 ln(1 + √ 2). The critical points of the two crossing-bond models -K c (1) = 0.1901926807(2) and K c (− 1 4 ) = 0.6972207(2)-were determined by means of a transfer-matrix technique combined with finite-size scaling 10 . This analysis confirmed with a high precision that the two crossing-bond models with belong to the static Ising universality class.
Monte Carlo averages were taken over 1.2 × 10 8 spin configurations, for system sizes in the range 5 ≤ L ≤ 20. For the nearest-neighbor model these samples were separated by a number of Monte Carlo steps per spin equal to one for L = 5 and increasing quadratically to ten for L = 20. For the other systems these numbers where multiplied by the appropriate scale factors. These surprisingly short intervals are possible because the convergence of the eigenvalue estimates as a function of projection time t in Eqs. (7) is governed by lower-lying Markov matrix eigenvalues. These are much smaller than the largest odd eigenvalue, which usually determines the relaxation rate. For the system size L = 5, the Monte Carlo results for the largest odd eigenvalues of the three models were compared with numerically exact results 5 . The consistency of both types of results confirms the validity of our numerical procedures.
As noted before for the largest odd eigenvalue of the nearest-neighbor model 5 , the high statistical accuracy of the Monte Carlo estimates of the eigenvalue is due to the accuracy of the approximation of the eigenvector of the Markov matrix by the optimized trial vectors. The present Monte Carlo results for the largest odd eigenvalues of the nearest-neighbor models agree with those of Ref. 5. The new data are based on statistical sample smaller by a factor of about 7, but the current trial vectors had more variational freedom.
For finite system sizes L we expect corrections to the leading scaling behavior τ L ∼ L z . Following Ref. 5, we assume corrections proportional to even powers of 1/L: where the series is truncated at order n c . Although we cannot exclude other powers in 1/L, we have used Eq. (8) to analyze the Monte Carlo autocorrelation times.
Results of such fits with n c = 3 are presented in Table I. The smallest systems do not fit Eq. (8) well for this value of n c . However, the residuals decrease rapidly when L 0 the smallest system size included in the fit is increased. The smallest acceptable value of L 0 , as judged from the χ 2 criterion, is also included in Table I.   TABLE I. Universality of the dynamic exponent z. Results of least-squares fits for the dynamic exponent for three Ising-like models and for five distinct relaxation modes, identified in the first column: ok refers to odd mode number k and ek refers to the corresponding even mode. Subsequent pairs of columns list L 0 , the smallest system size included in the fit, and the resulting estimates of z(β) for three ratios β = K ′ /K. Estimated errors are shown in parentheses. To account for possible lack of convergence as a function of projection time t and flaws in Eq. (8), two standard errors are quoted. The estimates of z obtained from the largest odd eigenvalues for the three models shown in Table I Table I as a confirmation of universality of the dynamic exponent for different models and modes of relaxation.
Correlation-time amplitudes were obtained from least-squares fits using Eq. (8) with z fixed at 13 6 , which happens to be close to the most accurate results in Table I. Subsequently, the non-universal metric factors m t were computed by fitting to Eq. (1). Defining m t (1) ≡ 1, we found m t (− 1 4 ) = 2.391 ± 0.002 and m t (0) = 1.5572 ± 0.0005. Table II shows results of the fits. Figure 1 is a semi-logarithmic plot of the effective, size-dependent amplitudes A Li (β) ≡ τ Li L −z /m t derived from the spectral gaps of the Markov matrices of the opposite-, nearest-, and equivalent-neighbor Ising models, β = − 1 4 , 0 and 1. Finally, we note that if one suppresses all but the magnetization dependence of the optimized trial functions, the number of nodes of the resulting functions equals the number of the corresponding eigenvalue counted from the top of the spectrum, which is in agreement with the odd-even alternation shown in Tab. II. This research was supported by the (US) National Science Foundation through Grants DMR-9214669 and CHE-9625498 and by the Office of Naval Research. This research was conducted in part using the resources of the Cornell Theory Center, which receives or received major funding from the National Science Foundation (NSF) and New York State, with additional support from the Advanced Research Projects Agency (ARPA), the National Center for Research Resources at the National Institutes of Health (NIH), IBM Corporation, TABLE II. Universality of correlation-time amplitudes. Results of least-squares fits for the finite-size amplitudes for three Ising-like models and for five distinct relaxation processes. The first column and the ones labeled L 0 are as in Table I. The columns labeled A i (β) contain the amplitudes defined in Eq. (1) for three interaction ratios β = K ′ /K with metric factors m t as given in the text. Estimated errors, as defined in Table I, are shown in parentheses. The difference A i (1) − A i (β) divided by its error is denoted by r. ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ ✸ + + + + + + + + + + + + + + + + ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷ ✷