text:
# Multiresolution analysis in statistical mechanics. II. The wavelet transform as a basis for Monte Carlo simulations on lattices

## I Introduction

One of the fundamental challenges in the simulation of very large systems is balancing the competing aims of accuracy, both numerical and physical, and computational efficiency. As the number of degrees of freedom, the number of interactions, and the number of system parameters become large, the computational cost and storage requirements of any algorithm which models these systems rapidly become prohibitive. For a system with $`N`$ degrees of freedom, the complexity of simulation algorithms is typically $`O\left(N^2\right)`$, although this can be reduced to $`O\left(N\right)`$ for methods incorporating both cell and Verlet lists and $`O\left(N^{3/2}\right)`$ for methods including Ewald sums, or increased to $`O\left(N^3\right)`$ or more for quantum methods. Consequently, very large systems are still expensive to simulate, even with efficient algorithms. One way to reduce the complexity of such systems is to employ a coarse-graining method, which systematically reduces the number of degrees of freedom in the system, and thereby the overall complexity of the simulation. Coarse-graining techniques have been developed for both on- and off-lattice calculations. Lattice simulations have generally relied on coarse-grainings based on renormalization group theory, while off-lattice simulations include coarse-graining techniques as varied as united-atom models, mapping to lattice models,  dissipative particle dynamics,  and bead-and-spring models.

At the same time, any benefits obtained from applying a suitable coarse-graining technique must be weighed against the principal difficulty of using such a method—inaccurate numerical results, such as a different critical point. The use of uncontrolled approximations also presents a problem, since even if a simulation converges to give an answer, we often cannot use results from those simulations to provide any insight into the systems about which we are actually interested. In addition, in virtually all cases the transformation is irreversible: we cannot reconstruct our original system after we have “evolved” a coarse-grained system.

The work presented in this paper represents the application of a new technique to enhance the performance of traditional lattice simulations using the wavelet transform method. The theoretical foundations of this technique are outlined in the preceding paper (henceforth denoted as I). First developed as an analysis technique which has found wide acceptance in signal processing, the wavelet transform has not yet been extensively applied to the field of molecular simulations. Exceptions to this include the variational work of Best and Schäfer, applied to statistical field theory, as well as the work of Huang in studying self-similarity in high-energy physics.

We choose to study the class of Ising lattices with pairwise and external field interactions, so that the Hamiltonian can be written in the form

$$=\underset{i}{}h_i\sigma _i\underset{i}{}\underset{k}{}J_{ik}\sigma _i\sigma _k,$$
(1)
where the indices $`i`$ and $`k`$ run over all spins in the system, and the interaction strengths $`h_i`$ and $`J_{ik}`$ can vary with position in the lattice. The nearest-neighbor Ising lattice corresponds to $`J_{ik}=0`$ unless $`\sigma _i`$ and $`\sigma _k`$ are neighboring spins on the lattice, in which case $`J_{ik}=J`$. Our goal is to construct a phase diagram for such systems as efficiently as possible. The original problem would be studied using the Metropolis Monte Carlo (MMC) technique, along with possible improvements such as the Swendsen-Wang or Wolff spin-cluster algorithms. We introduce a new method for studying this problem, using the wavelet transform technique to simulate the system hierarchically, by systematically transforming (1) to provide expressions for potentials at coarser scales.

Since the systems under consideration in this paper are lattices, we present the wavelet transform using the discrete framework, involving filter banks and matrices. In the discussion below, we present equations for one-dimensional systems only; however, the simulations are performed using higher-dimensional wavelets created using the lifting scheme of Sweldens.  The lifting scheme can be applied recursively as needed to achieve representations of the data at multiple length scales. However, we do not move between different physical models for the system, but only change the length scales which the different variables describe. As a result, the wavelet transform can be described as a multiresolution technique; multiresolution techniques comprise a subset of the more general class of multiscale techniques, which can describe any simulation involving multiple modeling steps, possibly involving the use of multiple underlying physical models.  We refer the reader to paper I for further details on the implementation of the wavelet transform with respect to lattice systems.

## II The Wavelet-Accelerated Monte Carlo (WAMC) algorithm

The principal difficulty of performing a wavelet transformation on a lattice system is in working with the discrete set of values that each spin is permitted to take, such as in a spin-$`q`$ Ising model. Because the transformed variables are linear combinations of the original variables, the constraint that the spins of the individual sites on the original lattice must be drawn from the set $`\{q,q+1,\mathrm{},q1,q\}`$ quickly becomes a more complicated constraint on the transformed variables $`\stackrel{~}{s}_i`$ and $`\stackrel{~}{\delta }_i`$. As the system size becomes large, the difficulty of rewriting the spin constraints proves so great that previous investigations of the use of wavelets in statistical field theory ignored Ising models altogether.  Consequently, we would like, if at all possible, to avoid computations involving original states after we have carried out the wavelet transformation. At this point, we take note of the application of wavelets to image compression, where the goal is to reduce the amount of information needed to reconstruct an image. We would like to apply this technique to lattice systems, and reduce the number of degrees of freedom which must be accounted for in our calculations.

We consider our system to be a $`d`$-dimensional regular lattice $``$ with side length $`l`$, so that the size of the lattice is $`N=\left|\right|=l^d`$, and we let a site $`\sigma _i`$on the lattice $``$ be characterized by a “spin” chosen from a finite set $`𝒥`$ of values and by its physical location on the lattice. For the spin-$`\frac{1}{2}`$ Ising model, for example, the set $`𝒥`$ is just $`\{+\frac{1}{2},\frac{1}{2}\}`$ (although for computational convenience this is usually treated as $`\{+1,1\},`$ a convention which we follow below); similarly, for a lattice gas based on a spin-$`1`$ Ising model, $`𝒥=\{0,1,2\}`$ represents the allowed occupation numbers of each lattice site. We then assume that the only physical interactions that occur are either interactions with an external field $`h_i`$ which can vary at each lattice site, or pairwise interactions with the bilinear form $`U(\sigma _i,\sigma _j)=J_{ij}\sigma _i\sigma _j`$, where $`J_{ij}`$ is usually a function only of the spacing between sites $`i`$ and $`j`$. Consequently, the Hamiltonian of the system can be written in the form (1),

$$\beta =\underset{i}{}h_i\sigma _i+\underset{i}{}\underset{j}{}J_{ij}\sigma _i\sigma _j.$$
(2)
For the purposes of our simulations, however, we will find it more convenient to treat the set of spins $`(\sigma _1,\mathrm{},\sigma _N)`$ and the external field $`(h_1,\mathrm{},h_N)`$ as vectors $`𝐮`$ and $`𝐡`$, and the pairwise interaction strengths $`J_{ij}`$ as a matrix $`𝐉`$. Then the Hamiltonian (2) can be written in matrix form as

$$\beta =𝐡^T𝐮+𝐮^T\mathrm{𝐉𝐮}.$$
(3)
This formulation of the problem is similar in spirit to that of graph theory, where the pairwise potential $`J_{ij}`$ is used to generate an adjacency list which specifies which edges interact.  Using (3) as the basis for a Monte Carlo simulation requires the calculation of the change of energy $`\mathrm{\Delta }E_{nm}`$ from microstate $`𝐮_m`$ to microstate $`𝐮_n`$:

$$\mathrm{\Delta }E_{nm}=𝐡^T\left(𝐮_n𝐮_m\right)+\left(𝐮_n𝐮_m\right)^T𝐉\left(𝐮_n𝐮_m\right).$$
(4)
If moves are restricted to changes of single spin flips, then only a single entry of $`𝐮_n𝐮_m`$ is nonzero, and therefore the calculation (4) reduces to a dot product, instead of a matrix multiplication.

As described in I, the action of the wavelet transform is to insert between each product in (3) or (4) the identity matrix in the form $`𝐈=𝐖^T𝐖,`$ where $`𝐖`$ is the wavelet transform which maps data from one scale to the next coarser scale, containing half as many data points. The resulting expressions rewrite the Hamiltonian in terms of wavelet-transformed averages and differences, with downsampling needed to reduce the number of variables from $`2N`$ to $`N`$. As before, the wavelet transform can be iterated by applying it to successive sets of averages, leading after $`K`$ iterations to Hamiltonians of the form:

$$\beta \stackrel{~}{H}=\left(\stackrel{~}{𝐡}^{\left(K\right)}\right)^T\stackrel{~}{𝐮}^{\left(K\right)}+\left(\stackrel{~}{𝐮}^{\left(K\right)}\right)^T\stackrel{~}{𝐉}^{\left(K\right)}\stackrel{~}{𝐮}^{\left(K\right)},$$
(5)
where in (5) the $`\stackrel{~}{𝐮}^{\left(K\right)}`$ represent “block spins” whose values are determined by wavelet averaging over some well-defined region of the original system. The Hamiltonians (3) and (5) have the same formal structure, so that Monte Carlo simulations of the two systems are essentially identical. The only modifications needed to simulate a coarse-grained Hamiltonian are the ability to select new microstates $`\stackrel{~}{𝐮}_i^{\left(K\right)}`$ which are generated through wavelet transformations of the original microstates $`𝐮_i,`$ and the elimination of unwanted degrees of freedom from (5). It should be noted that in (5), the elements of $`\stackrel{~}{𝐮}^{\left(K\right)}`$ are not restricted to the same values as in the original system, but are free to take on any value which is consistent with the wavelet transform applied to the system.

From above, we saw that for even Hamiltonians $`\stackrel{~}{H}\left(\stackrel{~}{𝐱}\right)`$, we should have that $`\delta _{\stackrel{~}{H}}=0`$ for any wavelet difference $`\delta ,`$ where $`_{\stackrel{~}{H}}`$ denotes the ensemble average weighted by the Hamiltonian $`\stackrel{~}{H}`$. As a “worst-case scenario” for our method, we shall assume not only that $`\delta _{\stackrel{~}{H}}=0,`$ but also that any terms in the Hamiltonian (5) containing fluctuation terms can be neglected as well. This assumption allows us to reduce the size of $`\stackrel{~}{𝐉}^{\left(K\right)}`$ from $`N\times N`$ to $`2^{Kd}N\times 2^{Kd}N`$, where $`d`$ is the lattice dimensionality. Consequently, instead of performing calculations involving all of the original variables $`𝐱`$ which describe the state of our system, we consider functions only of local averages of our original variables. However, we anticipate that this simplification of the interactions present in the system will have a significant impact on the thermodynamic behavior of the resulting system; we illustrate these effects below.

To generate the new microstates $`\stackrel{~}{𝐮}_i^{\left(K\right)}`$, we need an estimate for the probability distribution $`p\left(\stackrel{~}{u}_i^{\left(K\right)}\right)`$ which describes the individual sites in the coarse-grained lattice. Determining the correct distribution for a given $`\stackrel{~}{u}_i^{\left(K\right)}`$ would require a detailed simulation of the original system. An alternative, ignoring the effect of neighboring block spins, would be to perform an exact enumeration of the spins within a block, which is possible only for the smallest of block spins. Since we would like to apply this method to systems of arbitrary size, we want to avoid both of these options. Therefore, we simulate a sublattice with the same dimensions as $`\stackrel{~}{u}_i^{\left(K\right)}`$, ignoring physical interactions with the rest of the system by using either free or periodic boundary conditions. Using the standard Metropolis acceptance criterion, we compute distributions for the properties of the small lattice, such as the magnetization. Then, according to the matrix formulation described in Sections II and III of Paper I, since the wavelet transform defines a single block spin $`\stackrel{~}{u}_i^{\left(K\right)}`$ as a linear function of the individual spins at level $`K1`$ which it replaces, we can use the linearity properties of probability distributions to convert the distribution of the properties directly into a distribution for the block spin $`\stackrel{~}{u}_i^{\left(K\right)}`$. Finally, using the distribution for the block spin $`\stackrel{~}{u}_i^{\left(K\right)}`$ as a starting point, we perform a Monte Carlo simulation on the system of block spins defined by the Hamiltonian (5).

Although (3) and (5) are structurally the same, we cannot impose a one-to-one correspondence between the states in the configuration space of (4) and the states in the configuration space of (5). Consequently, the thermodynamic information obtained from the two will not necessarily be identical; as we have shown in paper I, there is under fairly broad conditions a loss of entropy associated with the application of coarse-graining to a system. We can ensure that the detailed balance condition for the simulation based on (5) is satisfied for the new simulation by requiring

$`{\displaystyle \frac{\alpha \left(\stackrel{~}{𝐮}_m^{\left(K\right)}\stackrel{~}{𝐮}_n^{\left(K\right)}\right)}{\alpha \left(\stackrel{~}{𝐮}_m^{\left(K\right)}\stackrel{~}{𝐮}_n^{\left(K\right)}\right)}}={\displaystyle \frac{p\left(\stackrel{~}{𝐮}_n^{\left(K\right)}\right)}{p\left(\stackrel{~}{𝐮}_m^{\left(K\right)}\right)}}e^{\beta \left(\left(\stackrel{~}{𝐮}_n^{\left(K\right)}\right)\left(\stackrel{~}{𝐮}_m^{\left(K\right)}\right)\right)},`$

where $`\alpha \left(mn\right)`$ is the probability of accepting a move from microstate $`m`$ to microstate $`n,`$ and $`p\left(m\right)`$ is the probability of selecting microstate $`m`$ as determined from simulations on finer-grained lattices at lower scales.

## III Theoretical performance of WAMC versus traditional MC

The wavelet transform is a hierarchical method which can be applied iteratively to a system to obtain successively coarser descriptions of a system. To describe the operation of the wavelet transform on a lattice model, we need to introduce some notation based on the various length scales in the problem. In the original problem, the applicable length scales are the lattice spacing $`l`$, the correlation length $`\xi ,`$ and the total lattice size $`L`$. Applying the wavelet transform method once increases the lattice spacing by some factor $`a`$, so that the ratios of correlation length to lattice spacing and of system length to lattice spacing each decrease by $`a`$. If we apply the wavelet transform $`m`$ times in succession, the corresponding factor becomes $`a^m`$.

We perform the simulation in a series of $`K`$ stages, where the length scales at each stage are functions of the length scales at the previous stages. The initial simulation is performed on a sublattice of the original problem, with lattice size $`L^{\left(1\right)}<L`$, where the superscript denotes the first stage of the simulation. The lattice spacing of the first stage is the same as in the original problem, so we define $`l^{\left(1\right)}=l`$. At each subsequent stage of the simulation, the lattice spacing of the $`k`$th stage is defined by the recursive relation $`l^{\left(k\right)}=L^{\left(k1\right)}l^{\left(k1\right)}`$. Since $`l^{\left(1\right)}`$ is fixed to be the lattice spacing of the original lattice, the adjustable parameter in this relation is the lattice size $`L^{\left(k\right)}`$ of each stage. If we assume that the lattice is the same length in all dimensions at every stage, a single variable in stage $`k`$ is a block variable representing the $`\left(L^{\left(k1\right)}\right)^d`$ variables simulated in stage $`k1`$.

Assuming that the lattice is the same length in all directions both in the original problem and at every stage in the wavelet-transformed problems, there are $`N_t=L^d`$ lattice variables in the original problem, and $`N^{\left(k\right)}=\left(L^{\left(k\right)}\right)^d`$ lattice variables in the $`k^{th}`$ stage of the wavelet-transformed problem. However, each variable in stage $`k`$ is a block variable representing the average behavior of the $`\left(L^{\left(k1\right)}\right)^d`$ variables in a block at stage $`k1`$, so the number of total degrees of freedom represented at stage $`k`$ is $`N_t^{\left(k\right)}=_{i=1}^kN^{\left(i\right)}`$, where $`N_t^{\left(K\right)}=N_t.`$ The number of simulated degrees of freedom is $`N_s=N_t`$ for traditional Metropolis Monte Carlo (MMC), but $`N_s=_{i=1}^KN^{\left(i\right)}`$ for WAMC. Because the running time of Monte Carlo simulations is usually linear in the number of degrees of freedom being simulated, the advantage of coarse-graining the system using a wavelet transform becomes evident. For example, consider an “original problem” of simulating a cubic lattice with $`256`$ Ising variables on a side. If we divide the original problem into two stages consisting of cubes of $`16`$ Ising variables on a side, we reduce the original problem of analyzing $`256^3=16,777,216`$ variables to the simpler problem of analyzing $`2\left(16^3\right)=8192`$ variables. Although it is more difficult to produce a trial configuration in a simulation of the wavelet-transformed problem than in a simulation of the original Ising lattice problem, this is more than offset by the reduction in the number of degrees of freedom being simulated.

## IV Results

For the purposes of comparison, our “experimental” systems are two-dimensional Ising models of size $`32\times 32`$, where we have run both MMC simulations on the full lattice, and WAMC simulations at a variety of resolutions; we shall denote these resolutions using the notation $`(x,y)`$, where $`x`$ indicates the length of the block size simulated in the first stage to estimate the probability distribution $`p\left(\stackrel{~}{𝐮}^{\left(K\right)}\right)`$ to be used in the second stage, and $`y`$ denotes the number of blocks on a side of the lattice in the second stage of the simulation.

### A Order parameter

Usually, the property of greatest interest in a simulation of a lattice system is the order parameter $`\eta `$. For spin systems, $`\eta `$ is generally taken to be the magnitude of the average magnetization, so that $`\eta =m`$. \[For XY and Heisenberg models, and other models where spins are oriented, we generally consider only the magnitude of the average vector $`\eta =m=\left|𝐦\right|`$.\] Generally, this is a very simple property to compute, since the value of the order parameter is constantly updated during the course of the simulation, and is thus always available.

For the $`32\times 32`$ Ising model, the results of a MMC simulation, as well as $`(4,8)`$\- and $`(8,4)`$-WAMC simulations are shown as Figure 1. The primary difference in the curves for the three cases is that as the coarse-graining process decreases the number of degrees of freedom in the final stage of the simulation, the location of the Curie temperature, indicating onset of spontaneous magnetization, increases and the steepness of the curve below the Curie temperature decreases. This result is consistent with our findings for average absolute magnetization $`\left|m\right|`$ from analytical models, discussed in paper I. In the present case, we note further that we achieve agreement between the different models not only in the low-temperature region, but also in the high-temperature regime $`TT_c`$. The differences in the intermediate regime can be attributed largely to the difference in behavior that results from the use of the wavelet transform to move from the original Hamiltonian (3) to a coarse-grained Hamiltonian (5). Additionally, the increased noise in the WAMC results at intermediate and high temperatures arises because of the approximations used for the probability distributions $`p\left(𝐮^{\left(K\right)}\right)`$ at the second stage of the simulation. The relative lack of noise in the MMC results stem in part from the fact that the Metropolis technique leads to non-ergodic sampling of phase space as temperature increases, as the simulation tends to cycle through a limited number of states.

### B Internal Energy

Plotting the internal energy $`U`$ as a function of the temperature, we obtain curves that follow the same general pattern outlined in paper I. As illustrated in Figure 2, at low temperatures, the internal energy, as computed for the $`32\times 32`$ model using standard MC as well as $`(4,8)`$\- and $`(8,4)`$-WAMC simulations, is in exact agreement for all methods. This occurs because only a few microstates of the system, corresponding to states that have all spins aligned, are actually observed by the system, and the wavelet transform preserves the energy of these states exactly. All three eventually reach an average internal energy of zero, but exact agreement is only expected in the infinite-temperature limit, when the difference in energy levels between microstates becomes unimportant. For intermediate temperatures, as before, the disagreement is a result of the change in form of the Hamiltonian that results from neglecting local correlations. Also, we note that for WAMC the “noise” in the internal energy increases both with increasing proximity to the “observed” critical point of the system as well as with increasing coarse graining. The additional coarse graining yields a Hamiltonian with reduced numbers of energy levels, since the energy of a block spin is defined here to be a function only of its overall magnetization, and not of its internal magnetization fluctuations; the reduced number of discrete energy levels yields noisier data.

### C Fluctuation properties

Fluctuation properties are useful for locating critical points since in the vicinity of a critical point, the magnitude of fluctuation properties is known to diverge as $`t^\alpha `$, where $`t\left|TT_c\right|/T_c`$.  Thus, a rapid increase in the value of a fluctuation property such as the heat capacity at constant external field, $`C_H=\left(E^2E^2\right)/k_BT^2`$, with respect to temperature can be used to estimate the critical temperature of a system. However, use of the wavelet transform leads to a decrease in the magnitude of the heat capacity, since the coarse-graining leads to smaller variances in the distribution of the energy $`U`$. Consequently, the maximum value of the heat capacity $`C_{\mathrm{max}}`$ decreases as a function of the number of degrees maintained in the problem.

In Figure 3, the heat capacity is shown as a function of dimensionless temperature $`k_BT/J`$ for the same systems as for the order parameter and internal energy measured above. We see that the location of the maximum of the heat capacity does increase, as expected. Although the relative maximum of the heat capacity obtained from the $`(4,8)`$\- and $`(8,4)`$-WAMC simulations appears to be identical, they differ by about 3 percent. Moreover, the actual value of the maximum is not as important as its existence and its location as a function of $`k_BT/J`$ and of the resolution of the model.

### D Scaling results

One application of the wavelet-accelerated MC method is to provide an upper bound for locating phase transitions. Running multiple simulations, at different levels of resolution, one can determine for each level of resolution the approximate phase transition temperature $`T_p\left(N_s\right)`$, where $`N_s`$ is the number of degrees of freedom (here, block spins) in the given system. From these data, one can extrapolate a scaling relationship of the form

$$\left(T_pT_p^{}\right)N_s^\gamma ,$$
(6)
where $`\gamma `$ is the corresponding “scaling” exponent, and $`T_p^{}`$ is the phase transition temperature for the untransformed model. The estimate obtained for the scaling exponent $`\gamma `$ depends upon the technique used to calculate the transition temperature for a given system—for example, estimating divergence of the heat capacity versus the onset of spontaneous magnetization. Our simulations suggest that $`\gamma `$ is typically between $`0.20`$ and $`0.25`$, with lower values obtained from divergence of heat capacity than from the onset of spontaneous magnetization.

As explained in I, using a relationship like (6) to estimate the phase transition temperature will usually lead to an overestimate of the phase transition temperature. This is a consequence of the underestimation of entropy that occurs through the reduction of the size of configuration space as a result of the wavelet transform. Estimates of $`k_BT_p`$/$`J`$, as determined by (6) for the two-dimensional Ising model considered here typically varied between about $`2.7`$ and $`2.9`$, which is an error of approximately 25 percent from the theoretical value of $`2.27`$ provided by the Onsager solution, but only about 20 percent from the results determined by the traditional MC simulations, which gave $`k_BT_p/J2.35`$. Although these errors are somewhat sizable, it is useful to note that the total computation time required to obtain the estimate using scaling laws is at least an order of magnitude smaller than the computation time required to perform a direct simulation on the original system. Thus, if computational time is at a premium, an effective approach may be to use the wavelet transform method to provide an upper bound for $`T_p`$, and then perform a direct simulation for the parameter space with temperatures below $`T_p`$.

### E Decorrelation time

Another important measure to study is the time required for decorrelated samples. It is well known that in the vicinity of the critical point, traditional Monte Carlo algorithms experience so-called “critical slowing-down.”  The Monte Carlo aspect of the WAMC algorithm does not vary from traditional Metropolis Monte Carlo, so we expect that the performance of the two algorithms should be similar, when measured near their respective critical temperatures.

To compare the two methods, we generated $`2^{24}=\mathrm{16\hspace{0.17em}777\hspace{0.17em}216}`$ new configurations for the $`32\times 32`$ Ising model at the critical point using traditional Metropolis Monte Carlo, as well as for the $`(4,8)`$-, $`(8,4)`$-, and $`(16,2)`$-WAMC models. To determine the correlation time, we used the blocking technique of Flyvbjerg and Petersen.  The results are shown as Figures 4 and 5 for the MMC and $`(16,2)`$-WAMC models, respectively. The salient feature in the graph is the onset of a plateau in the value of the variance of the energy $`\sigma _E^2`$; according to the method of Flyvbjerg and Petersen, this indicates that configurations separated by a distance of $`2^x`$ steps are statistically independent, where $`x`$ is the number of blocking transformations that have been performed. For the MMC model, we find that $`x13`$ or $`x14`$ provides a decent estimate; for the $`(16,2)`$-WAMC model, $`x=16`$ is a good estimate for the index. \[For the $`(4,8)`$\- and $`(8,4)`$-models (not shown), $`x=15`$ is a reliable estimate.\] In each case, this indicates that between $`2^{14}=\mathrm{16\hspace{0.17em}384}`$ and $`2^{16}=\mathrm{65\hspace{0.17em}536}`$ steps are required between independent configurations. Thus, we conclude that there is no degradation of performance near the critical point of a WAMC simulation, relative to traditional MMC simulations.

## V Analysis

### A Measured performance comparison

In comparing the performance of the standard Monte Carlo algorithm to the wavelet-accelerated Monte Carlo algorithm, we performed $`5\times 10^5`$ lattice passes on a $`32\times 32`$-lattice on a 733 MHz Pentium II: for the standard Monte Carlo algorithm, this meant that, on average, $`5\times 10^5`$ attempts were made to flip each spin. For the WAMC algorithm, $`5\times 10^5`$ attempts were made to flip a spin on a given level. The results are summarized in Table 1. We see that the $`(4,8)`$\- and $`(8,4)`$-simulations, which have the smallest total number of lattice sites ($`80`$), perform the fastest; however, even the $`(16,2)`$-simulation, which has one-fourth as many variables ($`260`$) as the $`32\times 32`$ standard Monte Carlo simulation ($`1024`$), finishes in less than $`8`$ per cent of the time required for the latter simulation. As the system size increases, the computational efficiency achieved by breaking down the system into multiple stages, all of relatively equal size, becomes even greater: for a $`128\times 128`$-lattice, the performance gain increases from a factor of approximately 25 to a factor of approximately 50, when we compare the $`(32,4)`$\- and $`(16,8)`$-WAMC simulations to the standard MC model. However, for the $`(8,16)`$-WAMC model, the complexity of assigning one of $`65`$ possible values to each of $`256`$ variables according to the correct probability distribution becomes comparable to that of the original problem, so that in fact standard MC runs in roughly a factor of 3 faster than the $`(8,16)`$-model.

Similar results are also observed for calculating the phase diagram of a $`64\times 64`$-Ising lattice, which is shown in Figure 6 as a plot of average magnetization as a function of temperature and external field strength using an $`(8,8)`$-model, for temperatures between $`T=0.5`$ and $`T=5.0`$, and for field strengths between $`h=1`$ and $`h=1`$. The phase diagram reproduces the essential features of the original two-dimensional ferromagnetic Ising lattice, such as the phase separation at $`h=0,`$ although the exact shape differs from the results obtained via a standard Metropolis Monte Carlo simulations. However, the plot based on WAMC calculations is created approximately 40 times faster than would a comparable plot using standard MMC.

### B Comparison with renormalization group methods

Our observations also indicate that the accuracy of the wavelet-accelerated Monte Carlo simulations depends on the relative proximity to an “attractive fixed point” of the physical model in parameter space. Borrowed from renormalization group theory, these attractive fixed points represent the limiting behavior of the system under various conditions, such as the zero- and infinite-temperature limits and the limits of zero and infinite external field. As we approach these limiting cases, the approximations made in obtaining our wavelet-transformed Hamiltonian become increasingly less significant.

Combining these observations allows us to design an on-line fine-tuning algorithm for the coarse graining of our system: the further away from the critical point of the parameter space, the smaller the number of degrees of freedom $`N^{\left(K\right)}`$ necessary to simulate the system must be. Thus, if we keep track of changes in fluctuation properties such as the heat capacity or the magnetic susceptibility as we change the system parameters $`(h,T)`$, we can get an estimate of our relative distance to the critical point. If we are sufficiently far away from the critical point, we can choose either to increase the number of stages $`K`$ that we simulate, or we can look at more degrees of freedom at lower stages by increasing $`N^{\left(1\right)},N^{\left(2\right)},\mathrm{},N^{\left(K1\right)}`$. As we approach the critical point, we can either decrease the number of stages $`K`$ or include more degrees of freedom at higher stages by reducing $`N^{\left(1\right)},\mathrm{},N^{\left(K1\right)}`$.

### C Sources of error

In general, the source of our errors can be traced to the assumption that local fluctuation terms could be reasonably ignored in our coarse-grained Hamiltonian (5). This naive but otherwise useful assumption yields correct thermodynamic behavior when the overall physics of the system is particularly simple: in the low- and high-temperature regimes, for instance, when the number of observed microstates is small or when the differences between observed microstates is inconsequential. For more complicated behaviors, as found at intermediate temperatures and above all in the vicinity of a critical point, the use of this assumption has a drastic effect on both the phase space of the system, which in turn affects all thermodynamic properties of the system, including the internal energy and entropy of the system, as well as fluctuation properties of the system.

More complicated methods for dealing with fluctuation terms have significant drawbacks associated with them. Treating the fluctuation terms just like the block averages maintained in $`\stackrel{~}{𝐮}^{\left(K\right)}`$ means that the wavelet-transformed Hamiltonian is no simpler than the original Hamiltonian, which affords few advantages in computational time. Likewise, other approaches, such as parametrizing the probability distribution for the elements of $`\stackrel{~}{𝐮}^{\left(K\right)}`$ using a property like the energy $`E,`$ introduce new functional dependencies which cannot be taken into account using the wavelet transform. Thus, we sacrifice one of the major advantages of the method—moving from one level to another is achieved exclusively through use of the wavelet transform. Thus, the most promising avenue for dealing with fluctuation terms is to develop a probability distribution for the fluctuation terms via the same approach used to determine probability distributions for the local averages. Then, using the probability distribution for the fluctuation terms, we can treat the discarded terms of the Hamiltonian as a noise term which can be used to restore some of the entropy that was lost as a result of the coarse-graining \[see paper I for more details\]. However, we have presented our results here to show what can be achieved under “worst-case” conditions, without the use of inverse coarse-graining methods.

Many coarse-graining techniques control errors by fitting the parameters of a new Hamiltonian to ensure agreement with some known structural information about the system, such as the radial distribution function. For lattice systems, this iterative approach is reflected in renormalization group theory, and notably the Wilson recursion method, which finds the fixed points of the system. As formulated, WAMC creates a coarse-grained Hamiltonian by truncating the Hamiltonian obtained after application of the wavelet transform. As an improvement to this, it should be possible to use the wavelet transform to determine which terms will appear in the Hamiltonian, and then determine the appropriate parameters to ensure the best fit for some desired property of the system using an iterative approach.

### D Constructing an adaptive algorithm for MC using the wavelet transform

As pointed out above, the wavelet transform method tends to produce overestimates for the critical point of the system; therefore, if we start with the high-temperature limit of our algorithm and slowly reduce the temperature in our simulation, we can observe the movement toward the critical point by watching various fluctuation parameters, such as the heat capacity $`C_H=\left(E^2E^2\right)/k_BT^2`$. Near the critical point, we expect to see a rapid increase in the value of $`C_H`$, consistent with the expected logarithmic divergence observed in the limit of finite-size systems.  If we use the onset of this logarithmic divergence as an indicator, we can then “step down” and use a finer lattice including more degrees of freedom. This system will naturally better reflect the physics of our system, particularly in the vicinity of the critical point. We expect that very near the critical point, we will have to simulate the system at the original scale, since this will be effectively the only level which accurately represents the underlying behavior of the system. However, the region of parameter space where this is necessary is relatively small compared to the complete parameter space. This is especially true when we consider that as we proceed below the critical temperature of the system, the logarithmic divergence of $`C_H`$ will also vanish. As a result, as we move increasingly far away from the critical point, we begin to approach the other fixed-point behaviors associated with the low-temperature limits of the system. Since these are reasonably well-preserved using the wavelet transformation, we can safely return to increasingly coarse-grained descriptions of our system as the simulation proceeds past the critical point.

As an example, we compute the spontaneous magnetization curve for a $`64\times 64`$ Ising lattice in the temperature range $`0.5T10.0`$, with $`\mathrm{\Delta }T=0.05`$, and choosing as our refinement criterion $`\mathrm{\Delta }C_H/\mathrm{\Delta }T0.5`$, until we reach the finest scale, corresponding to the original problem. We begin by coarse-graining the system to an $`(8,8)`$-model, where we find that the criterion is triggered only at $`T=5.1`$; we then continue with a $`(4,16)`$-model, down to $`T=4.0`$, at which point the refinement criterion is exceeded. Refining once more, we proceed with a $`(2,32)`$-model until $`T=3.4`$, at which point the threshold is again crossed. Since the next refinement is the original problem, we proceed at this level of resolution until we have passed the critical point, so that $`\mathrm{\Delta }C_H/\mathrm{\Delta }T`$ is positive. As a coarsening criterion, we select for simplicity the opposite of the refinement criterion, $`\mathrm{\Delta }C_H/\mathrm{\Delta }T0.5.`$ Using this criterion, we find that we coarsen the model to the $`(2,32)`$-, $`(4,16)`$\- and $`(8,8)`$-models at temperatures of $`T=1.75`$, $`T=1.65`$, and $`T=1.55`$, respectively. The rapid coarsening of the model results from the higher estimates of the critical point in the coarsened models. Since we are well past the critical point, we expect changes in the heat capacity as a function of temperature to be relatively small, and thus it is possible to obtain accurate results from a relatively coarse model. Computationally, the time required to create this diagram was only 28 per cent that required to perform a standard Metropolis Monte Carlo simulation with the same number of steps. Moreover, in the regions that were not simulated using MMC, the computation time required was just 8 per cent of the time required for MMC. The resulting plot of magnetization versus temperature, shown as Figure 7, compares favorably to the analytical solution of Onsager, which is also shown.

## VI Conclusions

The WAMC algorithm can dramatically reduce the time required to calculate the thermodynamic behavior of a lattice system; the trade-off for these savings in time is in the accuracy of the results obtained, a general feature of coarse-graining techniques. The error that results is a function of position in parameter space: the results obtained are generally accurate in the vicinity of fixed attractors of the system, and decrease as one approaches critical points of the parameter space. Near critical points, deviations from results performed on the original lattice system are the result of coarse-graining the Hamiltonian by eliminating local fluctuation terms. Consequently, this suggests that a hierarchical simulation which uses fluctuation properties such as the heat capacity $`C_H`$ to gauge proximity to critical points in “real time” would yield dramatic savings in the computation time of the behavior of a lattice system over a wide region of phase space, as regions of space close to fixed attractors would be simulated at a very coarse scale, with full-detailed simulations reserved only for regions of parameter space close to critical points.

## VII Acknowledgments

Funding for this research was provided in part by a Computational Sciences Graduate Fellowship (AEI) sponsored by the Krell Institute and the Department of Energy.

Table 1. Performance comparison for Metropolis Monte Carlo (MMC) versus Wavelet-Accelerated Monte Carlo (WAMC) algorithms

| Simulation | Time for $`5\times 10^5`$ passes (s) |
| --- | --- |
| $`32\times 32`$ MMC | 1824.22 |
| $`(4,8)`$-WAMC | 64.5781 |
| $`(8,4)`$-WAMC | 55.1094 |
| $`(16,2)`$-WAMC | 144.516 |


format: text

text:
# How glasses explore configuration space

## I Introduction

Glassy systems spend long times out of equilibrium, i.e. in regions of configuration space which have vanishing Boltzmann probability. The evolution towards more and more optimized regions of phase space is so slow, that the fast degrees of freedom react on short time scales as if they were in equilibrium in the frozen background of the slow variables. On longer time scales however the off-equilibrium nature of the glassy phase shows up in the phenomenon of aging: the response to small perturbations, as well as the correlation functions, depend on the “age” of the glass, i.e. the time spent in the low temperature phase . The dynamics becomes slower and slower as the age becomes larger and nevertheless, even on the largest scales observed there is no tendency to stop: the system eventually wanders away from any finite region of phase space.

In the last years, following some developments of spin glass mean field theory , many theoretical , numerical  and more recently experimental papers  have focused on the study of off-equilibrium fluctuation-dissipation relations. In , it has been shown that the fluctuation-dissipation relation found in mean-field analysis of glassy relaxation can be used to define effective temperatures, higher then the one of the heat bath governing the heat exchanges between slow degrees of freedom.

The analysis of the fluctuation-dissipation relation has revealed deep relations between aging dynamics and the nature of the free-energy landscape . In this paper we would like to review a picture of glassy dynamics based on the analysis of these relations .

Two kinds of glassy systems emerge from experimental and theoretical analysis. One can identify two kind of glassy behaviors

* A first class, of which spin glasses are an example, where the asymptotic value of extensive quantities, like internal energy, and magnetization seem not to depend on the cooling rate or other differences in the cooling procedure.
* A second class, the one of structural glasses, in which the apparent asymptotic value of these quantities depends strongly on the cooling rate, and remain strongly different from the equilibrium values for the the largest time scales which can be probed.

Aging phenomena are common to both families. During aging the linear response to external perturbations displays age dependent behavior and anomalies with respect to equilibrium, which is conform to the fluctuation-dissipation theorem. Many efforts have been devoted to relate these anomalies to the properties of the underlying energy -or free-energy- landscapes .

A first step towards the comprehension of this relation can be achieved idealizing the first situation assuming that for long times all extensive “one time observable” (1TO), or more precisely the ones that can be written as sum of local quantities, are close to equilibrium. We call this situation ultra-long time limit and observe that in principle a) this limit does not imply the absence of aging in two time observables such as responses and correlation functions, b) usual nucleation arguments imply that in systems with regular short range interactions this limit is achieved in a finite, i.e. volume independent, time. Of course this last observation remains a matter of principle in structural glasses, where this equilibration time is exceedingly long, and one needs to understand dynamics on much shorter time scale.

For notational simplicity, in the theory we are going to expose, we will use the language of magnetic systems. Our considerations however will very general and can immediately applied e.g. to glass forming systems of classical particles in interaction. In order to illustrate the ideas we discuss a spin glass system where the spins $`S_x`$ interact through a Hamiltonian

$$H(𝐒)=\underset{𝐱<𝐲}{}𝐉_{𝐱,𝐲}𝐒_𝐱𝐒_𝐲.$$
(1)
on the D-dimensional square lattice and short range interactions $`J_{x,y}`$. The spins will be for simplicity assumed to be Ising variables $`S_x=\pm 1`$. We will consider the dynamical setting of a fast quench from high to low temperature at an instant marking the origin of the time axes.

## II Aging in the ultra-long time limit

The classical aging experiments concentrate on the study of the response of the system to an external perturbation in the linear regime. Consider the effect of an external magnetic field $`h`$, corresponding to a perturbation term in the Hamiltonian $`H^{(1)}(𝐒)=𝐡_𝐱𝐒_𝐱`$, acting from the quenching time $`0`$ to a waiting time $`t_w`$. At the linear order in $`h`$, the magnetization $`M(t,t_w)=\frac{1}{N}_xS_x(t)_h`$ at later times $`t`$ can be written as $`M(t,t_w)=h\chi (t,t_w)`$, where, introducing the instantaneous response $`R(t,s)=\frac{\delta M(t)}{\delta h(s)}`$, the susceptibility can be written as $`\chi (t,t_w)=_0^{t_w}𝑑sR(t,s)`$. This is a quantity usually measured in aging experiments and exhibit scaling, non-time translation invariant behavior. In usual equilibrium regimes, this function is related to the correlation function $`C(t,t_w)=\frac{1}{N}_{x,y}S_x(t)S_y(t_w)`$ by the fluctuation-dissipation relation $`\chi (tt_w)=\beta C(tt_w)`$. Off-equilibrium the same relation has in general no reason to be valid. One usually defines the fluctuation-dissipation ratio (FDR)

$$\beta \stackrel{~}{x}(t,t_w)=\frac{\chi (t,t_w)/t_w}{C(t,t_w)/t_w}$$
(2)
or $`x(t,C)`$ defined by the relation $`x(t,C(t,t_w))=\stackrel{~}{x}(t,t_w)`$. In equilibrium conditions the FDR is equal to unity and $`\chi (t,t_w)=\beta C(tt_w)`$. Spin glass mean field theory , has suggested that in aging systems $`x(t,C)`$ tends to a non trivial limit $`x(C)`$ taking the limit $`t,s\mathrm{}`$ for fixed values of $`C(t,s)`$. In these systems the FDR appears to have some universal character. It is independent of the particular response/correlation function probed. The function $`x(C)`$ appears to have mathematically the properties of a probability function. The ratio $`T/x(C)`$ can be considered an effective temperature governing heat exchanges on degrees of freedom reacting on the time scales specified by the value of the correlation . Such non generic behavior can be understood in the ultra-long time limit, resulting in a rather detailed statistical description on how glassy systems visit configurational space in aging dynamics. In the limit in which the FDR becomes time independent, one time observables (e.g. energy and magnetization) become asymptotically close to their final value. We will call this limit “ultra-long time limit” in reference to the fact that structural glasses are always observed on much shorter time scales. Let us stress that in short range systems with regular interactions the asymptotic values of the 1TO must be the equilibrium one. The ultra-long time limit while it could be appropriate in spin glasses where the values of 1TO’s seem not to depend on cooling procedure, but it is certainly not an appropriate description for structural glasses, where the 1TO are strongly out of equilibrium. In this context we can still think of this limit as a useful conceptualization to study aging in a simplified situation and getting hints on the statistical principles that govern glassy dynamics even far from this ideal situation.

## III Equilibrium ergodicity breaking

We would like to discuss how a non trivial FDR in the ultra-long time limit relates to the phenomenon of ergodicity breaking in the equilibrium distribution. Breaking of ergodicity, or absence of ergodicity at equilibrium means that the weight of the equilibrium distribution is concentrated on more then one disjoint phase space regions, called ergodic component, where the ergodic property holds separately. If the system is prepared in one of these regions it will never get out from it. A familiar example is provided by systems in presence of a first order phase transition, where different phases are equally dynamically stable. In that case suitable “external field” can be used to project on the different phases. The situation we would like to describe, that of an ideal glass, can be expected to be different. In disordered systems, the possible different ergodic components dominating at low temperature can be expected to be as disordered as the configurations of the high temperature phase and no external field is available to select any of them. A suitable description of that situation has been provided by spin glass mean field theory  where ergodicity breaking is characterized statistically through the comparison of the different ergodic components. One first define a measure of similarity among configuration, and then studies its statistical distribution induced by the canonical measure and the eventual quenched disorder. For spin glasses it is natural to define the overlap between two spin configurations $`𝐒=(S_1,\mathrm{},S_N)`$ and $`𝐒^{}=(S_1^{},\mathrm{},S_N^{})`$ as the normalized scalar product:

$$q(𝐒,𝐒^{})=\frac{1}{N}\underset{i}{}S_iS_i^{}$$
(3)
The overlap probability function (OPF), is defined as

$$P(q)=\overline{\frac{1}{Z^2}\underset{𝐒,𝐒^{}}{}\mathrm{exp}\left(\beta (H(𝐒)+H(𝐒^{}))\right)\delta (qq(𝐒,𝐒^{}))}.$$
(4)
where the over-line denotes the average over the quenched variables of the system (if any). If ergodicity holds the OPF equal to a single delta function. Any difference from this simple form is sign of ergodicity breaking. Conversely, a single delta function is not necessarily associated to ergodicity. The different components $`\alpha =1,2,3,\mathrm{}`$ appear in the equilibrium distribution each with a probability weight $`w_\alpha `$ reflecting its free-energy difference with the lowest state. It turns out that the OPF can be non-trivial only if the participation ratio $`_\alpha w_\alpha ^2`$ is non-vanishing in the thermodynamic limit. This excludes e.g. the case of an extensive number of ergodic component where $`w_\alpha \mathrm{exp}(N\mathrm{\Sigma })`$.

Despite its origin in the context of spin glasses, the concept of overlap, and the corresponding probability function, can be suitably to more general glassy systems. A possible definition in the case of structural glasses and some applications are discussed in .

## IV How Aging relates to equilibrium

Linear response theory allows to relate the OPF of the canonical probability and the FDR according to the relation :

$`P(q)={\displaystyle \frac{\mathrm{d}x(q)}{\mathrm{d}q}}.`$ (5)
The argument that leads to (5) is rather formal, and consists in finding a family of observables, expressed as sum of local quantities, whose average values relate respectively in equilibrium and in dynamics to the moments of the OPF and the derivative of FDR with respect to the correlation. As the dynamical averages should tend to the equilibrium ones, one deduces the identity between the OPF and the derivative of the FDR. The arguments runs as follow. Let us consider a set of operators which translate the lattice along a direction $`\widehat{ı}`$ parallel to one of the coordinate axes:

$$𝒯_k^{(p)}(x)=x+\frac{kL}{p}\widehat{ı};k=1,\mathrm{},p1.$$
(6)
We denote by $`𝒮_1`$ the set in which the coordinate in the direction $`\widehat{ı}`$ takes the values $`1,2,\mathrm{},L/p`$. We define our family of perturbations to have the form

$$H_p(𝐒)=\underset{x𝒮_1}{}J_x^{(p)}S_xS_{𝒯_1^{(p)}(x)}\mathrm{}S_{𝒯_{p1}^{(p)}(x)},$$
(7)
where the couplings $`J_x^{(p)}`$ are independent, identically distributed Gaussian variables, with zero mean and variance $`E_JJ_x^2=p`$. Despite the apparent long range character of the interactions, the perturbations $`H_p`$ can be seen as short range observable in a folded space, and as such they belong to the class of observables whole long time limit values tend according to our hypothesis to the equilibrium ones. Let us consider now a system evolving with Hamiltonian $`H+hH_p`$ These perturbations can be seen as short range perturbation and one can expect their effect to be described by linear response theory for small $`h`$. If we consider the expected values, in off-equilibrium dynamics and at equilibrium of the perturbation $`H_p`$, following simple mathematical manipulations implied by linear response theory and by self-averaging properties with respect to the variables $`J_x^{(p)}`$, one finds respectively for the dynamic and equilibrium averages

$`H_p(t)_{dyn}`$ $`=`$ $`\beta h{\displaystyle _0^1}𝑑qpq^{p1}x(t,q)`$ (8)
$`H_p_{eq}`$ $`=`$ $`\beta h(1{\displaystyle _0^1}𝑑qq^pP(q))`$ (9)
In the ultra-long time limit, $`H_p(t)_{dyn}H_p_{eq}`$ for all $`p`$. Correspondingly, $`x(t,q)x(q)`$ and integrating (8) by parts we find the relation (5). This relation relates the possibility of a persistent non-trivial FDR in the aging dynamics to to ergodicity breaking in the equilibrium measure, and as we will discuss, can be taken as the starting point for an analysis of how aging systems visit configuration space. The relation can be generalized to other important features which have been identified theoretically as possible in aging dynamics. One of this feature is ultrametricity which at equilibrium implies a hierarchical organization of the ergodic components. In dynamics this property means that the times $`t(C,t_w)`$ necessary for the correlation $`C(t,t_w`$ to take the value $`C`$ and defined by the relation $`C(t(C,t_w),t_w)=C`$ verify for $`C_1<C_2`$ the relation: $`t(C_2,t_w)/t(C_1,t_w)\genfrac{}{}{0pt}{}{}{t_w\mathrm{}}0`$, in other words that decreasing values of the correlation correspond to increasingly long relaxation times. Reasoning similar to the ones that lead to (5) allow to conclude that surprisingly in a given systems static ultrametricity and dynamic ultrametricity either are both present or both absent. Recently the relation (5) has also been generalized to local quantities in . A more complete discussion of the relation (5) should qualify the validity of the use of linear response in presence of ergodicity breaking, which involves the commutation of the thermodynamic limit and $`h0`$ in statics, as well as the long time limit and $`h0`$ in dynamics. This property has been called “stochastic stability”, in reference to the fact that the commutation of limits is possible whenever the introduction of weak random perturbations of the kind (7) has not major effects on the statistical distribution of the states relevant for the equilibrium probability and the long time off-equilibrium dynamics. On physical ground one expects stochastic stability to be valid and linear response theory to have a range of validity in glassy systems. Violation of this properties should result in dynamical crossovers for $`t\mathrm{}`$ and $`h0`$ which are in principle experimentally observable. So if in some systems one could measure the asymptotic value of the FDR this would lead informations on the equilibrium free-energy landscape.

## V Physical picture

The relation (5), which links off-equilibrium dynamics to the properties of the free-energy landscape relevant at equilibrium, has been discussed in the previous section on a purely formal ground: we have supposed equilibration of 1TO together validity of LRT and obtained the relation from straight mathematical analysis. The scope of this section is to discuss its physical origin. In order to simplify the discussion, we will consider the simplest model of aging dynamics with a non-trivial FDR. We will suppose that the dynamics will consists in a short time equilibrium like part and a long time aging part characterized by a unique age dependent relaxation time $`\tau (t_w)`$ which grows and tends to infinity for large $`t_w`$. Notice that this is the assumption commonly used to fit aging response data both in structural glasses  and in spin glasses . In this picture the correlation functions can be decomposed in a stationary short time part and an aging long time part as

$$C(t,t_w)=C_{st}(tt_w)+C_{ag}\left(\frac{tt_w}{\tau (t_w)}\right)$$
(10)
as usual we define $`q_{EA}`$ as the value that separate the stationary from the aging part of the correlation so that $`C_{st}(tt_w)=C_{st}^{}(tt_w)q_{EA}`$ is a monotonically decreasing function equal to $`1q_{EA}`$ for $`tt_w=0`$ and tending to zero for $`tt_w\mathrm{}`$, while $`C_{ag}(\frac{tt_w}{\tau (t_w)})`$ is equal to $`q_{EA}`$ for $`\frac{tt_w}{\tau (t_w)}=0`$ and tends to $`0`$ for $`\frac{tt_w}{\tau (t_w)}\mathrm{}`$. Correspondingly we will assume that the FDR will take the value 1 in the stationary domain ($`q_{EA}C(t,t_w)1`$) and the constant value $`x<1`$ in the aging regime $`C(t,t_w)q_{EA}`$, so that the linear susceptibility takes the form

$$\chi (t,t_w)=\beta C_{st}(tt_w)+\beta xC_{ag}\left(\frac{tt_w}{\tau (t_w)}\right)$$
(11)
. The relation between correlation and response is often visualized plotting parametrically $`\chi (t,t_w)`$ vs. $`C(t,t_w)`$ for fixed $`t_w`$. The present case corresponds to having for large $`t_w`$ two straight lines respectively of slope $`\beta `$ for $`q_{EA}<C(t,t_w)1`$ and $`\beta x`$ for $`0C(t,t_w)<q_{EA}`$. This form of the correlation and response is the one found in models of the $`p`$-spin model family . Numerical simulations of glass-forming systems indicate that the two regime behavior is a good approximation for the finite $`t_w`$ behavior of the response not too close to the estimated value of $`q_{EA}`$, although in that case the parameter $`x`$ depends on $`t_w`$ . In these cases however 1TO are far from equilibrium, differently from what we suppose here.

The form (11) corresponds to an equilibrium distribution where the function $`P(q)`$ has the simple form

$`P(q)=(1x)\delta (qq_{EA})+x\delta (q).`$ (12)
A straightforward computation shows that the observables $`H_p`$ of the previous section are expressed according to (9) as

$`H_p=\beta h(1q_{EA}^p+xq_{EA}^p).`$ (13)
The equilibrium regime at short times implies that the system spends its time in regions of the configuration space where it has the time of approximately equilibrate before relaxing further on a much longer scale. This observation can be put at the basis of decomposition where the dynamical variables, say the instantaneous values of the spins, are written as the sum of a fast and a slow contribution.

$$S_x(t)=[S_x(t)m_x(t)]+m_x(t),$$
(14)
where $`[S_x(t)m_x(t)]`$ is the fast component and the slow component $`m_x(t)`$ can be simply defined as

$$m_x(t)=\frac{1}{\tau (t_w)}_t^{t+\tau (t_w)}duS_x(u).$$
(15)
The set of the fast variables that correspond to the same slow variables define metastable regions of configuration space, called “quasi-states” that have a life time of the order of $`\tau (t_w)`$. This notion can be related to the potential energy landscape notions of basin or more precisely metabasins  in the inherent structure picture , with whom the quasi-states could be identified.

We can at this point rather naturally define thermodynamic quantities: to each quasi-state labeled by the values of the slow variables $`𝐦=(m_1,\mathrm{},m_N)`$ we can associate a free-energy $`f(𝐦)`$ and we can define a “configurational entropy” $`\mathrm{\Sigma }(f)`$ as the logarithmic multiplicity of quasi-states as a function of their free-energy $`f`$.

For long but finite times the quasi-states dominating the dynamics will present a small but extensive free-energy difference with the lowest state. The parameter $`q_{EA}`$ has the meaning of self-overlap within a quasi-state: $`q_{EA}=\frac{1}{N}_xm_x^2(t)`$. The equilibrium distribution on the other hand concentrates on the ergodic components, or true states, which have finite free-energy differences with the ground state and infinite life time. As we have stressed, if the function $`P(q)`$ is non trivial the multiplicity of such states grows less then exponentially with $`N`$. On the other hand, for the number of quasi-states with free-energy difference $`N\mathrm{\Delta }f`$ with the ground state, one can expect: $`𝒩(f)\mathrm{exp}(N\rho \mathrm{\Delta }f)`$.

Let us discuss the effect on the system of a small perturbation and concentrate for simplicity on the case of $`p=1`$, corresponding to an external magnetic field. As discussed first in , the effect of the field on the equilibrium distribution will be twofold. On one hand within any given ergodic component, with given average magnetization $`M_{av}`$, higher weight will be given to configurations with higher magnetization difference with the average $`MM_{av}`$. On the other, in the selection of the components higher weight will be given to components with higher average magnetization. The actual value the magnetization will take can be then understood as a two step free-energy minimization process. The first step is a single component free-energy minimization. This is done according to the usual equilibrium relation within a component $`\beta h(MM_{av})=\frac{1}{N}\beta [H_1^2_{s.c.}H_1_{s.c.}^2]=\beta h(1q_{EA})`$, where with $`_{s.c.}`$ we have denoted the average on a single component. The second step is the choice of the components with the appropriate average magnetization. The probability that a given unperturbed component has average magnetization $`M_{av}`$ is for small values of $`M_{av}`$ Gaussian with zero average and variance $`Nq_{EA}`$. Taking then into account the independence of the perturbation from the original Hamiltonian, we have that the number of quasi-states with free-energy difference $`\mathrm{\Delta }f`$ with respect to the unperturbed ground state and average magnetization $`M_{av}`$ is given by

$$𝒩(f,M_{av})\mathrm{exp}(N[\rho \mathrm{\Delta }fM_{av}^2/2q_{EA}])$$
(16)
In presence of the perturbation the actual value of $`M_{av}`$ will be that one minimizing $`\mathrm{\Delta }fhM_{av}`$, respecting the constraint that the number of states is higher then zero: $`\rho \mathrm{\Delta }fm^2/2(1q_{EA})0`$. One finds that the minimization is achieved in fact for $`\rho \mathrm{\Delta }fM_{av}^2/2(1q^p)=0`$ thus finding, comparing with (13) that $`\rho =\beta x`$.

The equivalence (9) tells that the same selection criteria should be valid asymptotically in off-equilibrium dynamics. The way the energy distributes among the different degrees of freedom in dynamics is asymptotically the same as at equilibrium. Free-energy not only governs equilibrium, but also dynamics: quasi-states of equal free-energy are selected asymptotically with equal probability. Notice that this conclusion do not depend on the observables we have considered in the analysis, according to our line of reasoning the same factor $`x`$ appears in the anomalous response to any perturbation. Under these circumstances one may wonder how the system can always be microscopically far from the true equilibrium components as the off-equilibrium behavior of correlation and response implies. This relates to the abundance of quasi-states extensive the free-energy difference $`N\mathrm{\Delta }f`$, which in order to have a non trivial FDR must be exponentially large. In that case a small $`\mathrm{\Delta }f`$ while implies small difference for the 1TO from the equilibrium values, it implies big microscopic differences between the quasi-states visited and the ergodic components.

This considerations, based on time scale separation and vicinity of 1TO to their equilibrium value does not of course directly apply to structural glasses. However, based on the analysis of the case of the $`p`$-spin model where 1TO block to values different from the equilibrium ones  one can hypothesize that the previous quasi-state selection principle could also be valid far from the ultra-long time limit. This is comforted by numerical simulations which have reported the existence of well defined FDR very far from equilibrium which are fairly constant on the aging regime  and do not appear to depend on different quantities which have been studied .

In this regime we probe in a region of parameters corresponding to finite configurational entropy, and as 1TO are strongly out of equilibrium, no connection can be made between the measurable values of the FDR and the equilibrium OPF. However, thanks to equiprobability also in this case the FDR can be related to some static characterization of the configuration space. At a given time, corresponding to a quasi-state free-energy $`f`$ the system will have reached with equal probability one of the $`𝒩(f)\mathrm{exp}(N\mathrm{\Sigma }(f))`$ possible quasi-states, where the configurational entropy $`\mathrm{\Sigma }(f)`$ is an increasing function of $`f`$. In order to rationalize the relation (13) one has just to suppose that the dynamics is dominated by the relaxation of the configurational entropy towards smaller values, and that a small perturbation in the Hamiltonian do not change its rate of reduction. In other words, at a given time $`t_w`$, the configurational entropy will take the same value $`\mathrm{\Sigma }(t_w)`$ both in absence and in presence of the perturbation. Writing then the total free-energy as an unperturbed term $`f`$ plus a perturbation term $`\beta hM_{av}`$ one can write that $`\mathrm{\Sigma }(t_w)=\mathrm{\Sigma }(f)M_{av}^2/2(1q_{EA})`$. Minimizing the total free-energy $`f=f\beta hM_{av}`$ with this constraint one recovers the formula (13) with the relation

$$\beta x=\frac{\mathrm{\Sigma }}{f}$$
(17)
during aging the effective temperature $`1/(\beta x)`$ verifies with the configurational entropy the same relation that the true temperature verifies at equilibrium with the total thermodynamic entropy. This kind of relation has been recently used to argue in favor of the the “Edwards measure” for lattice gas systems  and granular material under shear . Let us notice that even if this is not indicated explicitly in the notation, now the various parameters ($`f`$, $`q_{EA}`$ etc.) are slow functions of time. The equiprobability hypothesis, if in the ultra-long time limit it comes rather naturally from 1TO equilibration, is much harder to justify in when these quantities are out of equilibrium. At present they are not clear the physical principles that would explain it, although it has been tentatively related to a chaotic property with respect to the thermal noise  according to which two “clones”, generated doubling a given system at the time $`t_w`$ and evolving later with different realization of the thermal noise, would give rise to divergent trajectories.

We would like to stress that the analysis we have presented implies equivalence between the relaxation of field induced perturbations after removal of the field and the regression of spontaneous fluctuation  in a way that generalizes the classical Onsager argument for equilibrium systems . Both kind of conditions in fact imply free-energy minimization with a configurational entropy constraint, and within the two step relaxation model we have considered lead to the relation (11).

Up to now we have considered as a reference setting a simple quench from high to low temperature. We can generalize the present dynamical picture to more complex thermal histories retaining equiprobability of quasi-states with equal free-energy. We would get in this case a dynamical picture in agreement with multi-temperature thermodynamics described in , where the macroscopic state of a glass is specified by a suitable number of history dependent effective temperatures. The model discussed here corresponds to just an effective temperature in addition to the external one. According to direct numerical measurement of the FDR in glass-forming models this corresponds to a good approximation to the situation found in model liquid systems . Good indications in favor of the present picture come from numerical simulations of glass forming liquids. In ref.  the relation (17) has been tested with positive answer starting from a direct measure of the FDR and an estimate of the configurational entropy as the logarithmic number of inherent structure with given energy. Other interesting numerical evidence has been found in the Kob-Andersen kinetic model in  and in the related problem of granular materials under shear in . Experimentally there are clear evidence that violations of FDT are present on long time scales , but if these correspond to effective temperatures is not clear at the present stage.

## VI Conclusions

In this paper we have discussed a possible picture on how glassy system visit configuration space based on the analysis of the fluctuation dissipation relations during aging. The ultra-long time limit, possibly relevant for the case of spin glasses, gives us a case where we can understand in a good detail the meaning of anomalous response and effective temperatures defined from off-equilibrium fluctuation dissipation relations. The asymptotic value of anomalous response is related to the existence of a non trivial equilibrium OPF $`P(q)`$, which in turn implies ergodicity breaking in the equilibrium distribution. The meaning of the equivalence can be found in a dynamical principle of selection of the quasi-states in the asymptotic limit: quasi-states with equal free-energy are selected with equal probability in the dynamical process. Under these circumstances where the system is macroscopically close to equilibrium, the possibility of always being microscopically far from equilibrium and continue to age relates to the abundance of states, which in order to have a non trivial FDR must be exponentially large as soon in the free-energy difference $`\mathrm{\Delta }f`$ is finite.

This leads to a dynamical picture that can be generalized to the case in which the system is macroscopically out of equilibrium, assuming equiprobability of quasi-states of equal free-energy, and a rate of decrease of the configurational entropy independent of possible small perturbations. This hypothesis predicts the existence of FDR (effective temperatures) related to the growth of the configurational entropy with quasi-state free-energy and therefore independent of the particular couple of correlation and response measured. The experimental verification of this property, and the clarification of the physical principles leading to it are open problems for future research.

## VII Acknowledgments

The present contribution is largely based on papers  written in collaboration respectively with M. Mezard, G. Parisi and L. Peliti and with M.A. Virasoro, whom I warmly thank.


text:
# Comment on ”Vector and axial-vector mesons at finite temperature”.

## Abstract

It is shown, that the correlators of vector and axial-vector currents at finite temperature $`T`$ in order $`T^2`$ reduce to mixing of vacuum $`VV`$ and $`AA`$ correlators with universal mixing coefficient given by the parameter $`\epsilon =T^2/6F_\pi ^2`$, unlike the claim in the recent paper by Mallik and Sarkar .

In a recent paper  Mallik and Sarkar considered vector and axial-vector correlators at finite temperature in order $`T^2`$. Earlier, the same problem was treated in our paper . Mallik and Sarkar criticize the calculation method used in  and claim, that the results obtained in that paper are wrong. In this comment we show that this is not true.

Let us remind shortly the basic results of .

The thermal correlators ($`J=V,A`$)

$$C_{\mu \nu }^J(q,T)=id^4xe^{iqx}\underset{n}{}n|T\{J_\mu ^a(0_),J_\nu ^a(x)\}exp[(\mathrm{\Omega }H)/T]|n$$
(1)
were considered at low $`T`$ in order $`T^2`$. Here $`e^\mathrm{\Omega }=\underset{n}{}n|exp(H/T)|n`$ and the sum is over the full set of eigenstates of Hamiltonian $`H`$. The currents $`V_\mu ^a`$ and $`A_\mu ^a`$ are isovector vector and axial currents. The validity of the chiral effective theory (CET) – a low energy consequence of QCD – was taken as granted. The explicit form of the Hamiltonian was not specified. It was assumed, that pions are massless and current algebra relations were exploited. In the sum over $`n`$ in order $`T^2`$ it is sufficient to account for the contribution of two lowest states – vacuum and one-pion. Then the factor $`T^2`$ arises simply from the one-pion phase space. The contributions of higher excited states are suppressed as $`T^k`$$`k>2`$ or exponentially, $`e^{m_h/T}`$. It was found  that in order $`T^2`$ (at $`q^2=Q^2,Q^2T^2`$)

$$C_{\mu \nu }^V(q,T)=(1\epsilon )C_{\mu \nu }^V(q,0)+\epsilon C_{\mu \nu }^A(q,0)$$
(2)
$$C_{\mu \nu }^A(q,T)=(1\epsilon )C_{\mu \nu }^A(q,0)+\epsilon C_{\mu \nu }^V(q,0),$$
(3)
where $`\epsilon =T^2/6F_\pi ^2`$, $`F_\pi =91MeV`$ is the pion decay constant. Therefore, in order $`T^2`$ the thermal $`V`$ and $`A`$ correlators reduce to their vacuum values, but in the vector correlator there appears the admixture of the axial one and vice versa. The positions of the poles (and other singularities) of thermal correlators are the same as at $`T=0`$ (up to appearance of $`A`$ poles in $`V`$ channel and vice versa). The residues of the poles are shifted in a universal way given by the parameter $`\epsilon `$. In ref. more detailed arguments were presented and the results were generalized to correlators of other currents (baryon, etc.).

In their consideration of the same problem, Mallik and Sarkar obtained the equation, equivalent to (2) with a difference, that mixing coefficients are not universal and depend on the form of low-energy Lagrangian. They assumed that the terms containing $`\rho `$, $`a_1`$, $`\omega `$ and $`f_1`$ fields are added to the standard CET Lagrangian. Then, in their calculation of $`C_{\mu \nu }^V(q,T)`$ (2), besides the factor $`1\epsilon `$ in front of $`C_{\mu \nu }^V(q,0)`$ there appears a term proportional to $`g_1^2`$, where $`g_1`$ is the $`\rho w\pi `$ coupling constant. However, embedding of the $`\rho w\pi `$ interaction term into effective low-energy Lagrangian is not legitimate. As can be easily proved (using eq.(2.22) of ) this term violates $`SU(2)\times SU(2)`$ chiral symmetry, leads to appearance of additionaal term in the expression for axial current and, consequently, violates current algebra relations. Moreover, famous low-energy theorems, which follow from current algebra (Goldberger-Treiman, Adler, Gross-Lelwellyn-Smith, Weinberg sum rules etc.) would also be violated if such a term persists. Therefore, the disagreement of the Mallik and Sarkar results with eq.2 comes entirely from illegitimate form of the low-energy Lagrangian used in . If the $`\rho w\pi `$ interaction term in the Lagrangian is omitted, then the results of  exactly coincide with eq.2, found in .

Now, a few words about the criticizm of our paper , given in . It was claimed in , that the self-energy diagrams on Fig.1 (Fig.1b in , Fig.1c,d in ) were ignored in . In fact, there are two types of such diagrams: 1) Fig.1a, where the pion is virtual. Such diagram corresponds to vacuum correlator, is $`T`$-independent and is accounted in $`C_{\mu \nu }(q,0)`$ in (2),(3); 2) Fig.1b, where pion comes from thermal bath and is on the mass shell. For a pion on mass shell its interaction vertex is proportional to momentum (Adler theorem). Therefore, the contribution of the Fig.1b diagram is $`T^4`$ and may be omitted. (We repeat here the arguments of ref. for clarity).

To conclude: the results, no correlator mixing, obtained in  are valid only for the model Lagrangian used in , but not for QCD and its consequense - low energy chiral effective theory (CET). On the other hand our paper  based on QCD and CET shows correlator mixing at low temperature, The thermal $`V`$ and $`A`$ correlators mixing in order $`T^2`$ \- eq. (2),(3) was confirmed in recent paper . Hence the claim of the authors in  that our paper  is ”incorrect” is not correct.