format: text

text:
# The Approach to Ergodicity in Monte Carlo Simulations

## I Introduction

A goal of Monte Carlo (MC) simulations in statistical mechanics  is the calculation of ensemble mean values of thermodynamic quantities. Ensemble mean values are multidimensional integrals over configuration space

$$U=𝑑𝐱P(𝐱)U(𝐱),$$
(1)
where $`P(𝐱)`$ is the probability of finding a system in the state defined by $`𝐱`$, and the functional form of $`P(𝐱)`$ depends on the ensemble investigated. MC simulations usually generate a sampling of configuration space $`\left\{𝐱_k\right\}_{k=1}^K`$ by the use of a stochastic process with stationary probability $`P(𝐱_k)`$. The quantity $`U`$ evaluated at $`𝐱_k`$ is the output of the simulation $`U(𝐱_k)=U_k`$, and its arithmetic mean value $`\overline{U}`$, in principle, must approach the ensemble mean value. In this paper we refer to the set of configurations generated in a Monte Carlo simulation as a time sequence, and we study the behavior of these temporal sequences $`\left\{U_k\right\}`$ and their arithmetic mean, to understand better how MC simulations approach ergodic behavior. It is important to emphasize that there are two time variables to consider. The time variable $`k`$ labels the separate configurations generated in a Monte Carlo walk. Variations of properties with $`k`$ provide information about the short-time behavior of a MC simulation. The time variable $`K`$ labels the total length of the MC walk, and variations of computed properties with $`K`$ provide information about the convergence of the simulation on a long time scale.

Given an infinite time, the stochastic walker in a MC simulation visits every allowed point in configuration space. Ergodic behavior is reached when the length of the walk is sufficiently long to sample configuration space appropriately. In practice, this does not mean that the space has been densely covered but that every region with non-negligible probability has been reached. In such a case we can say that the simulation is effectively ergodic or that it has reached the ergodic limit.

For a finite walk, in the event of broken ergodicity , phase space is effectively disconnected. The different disconnected regions (called components) are separated by barriers of zero effective probability. If a stochastic walker starts its walk in one of these regions, it may not cross the barriers within the time of the simulation. If the simulation length is increased, some barriers may become accessible for the walker and phase space is better sampled. We can conclude that a time $`\tau `$ exists such that, for simulation lengths shorter than $`\tau `$, the walker becomes trapped in one of the phase space components. For simulation lengths much larger than $`\tau `$, phase space is effectively covered by the walker.

In this study we imagine a system having more than one time scale $`\tau _1\tau _2\mathrm{}\tau _\mathrm{\Lambda }`$. In a Monte Carlo simulation each scale comes from stochastic processes with different correlation times. A precise definition of the correlation times for Monte Carlo processes is given in Section III, but for the moment we can think of these correlation times as identical to physical time scales of the system under study. To understand these time scales more fully, it is useful to focus on an example. Prototypical of systems having such disparate time scales are atomic and molecular clusters. Typical cluster potential surfaces have many local minima separated by significant energy barriers.  The local minima can be grouped into basins of similar energies, with each basin separated from other basins again by energy barriers. At short Monte Carlo times a cluster system executes small amplitude oscillations about one of its potential minima. We can think of these vibrational time scales as the shortest time scales that define a cluster system. As the simulation time is extended the system eventually hops between different local minima within the same basin. The time scale for the first hops between local minima can be considered the next shortest time scale for the simulation. At still longer Monte Carlo times, the system hops between different energy basins defining yet another time scale for the simulation. This grouping of time scales continues until the longest time scale for a given system is reached. At Monte Carlo times that are long compared to this longest time scale, the simulation is ergodic.

Consider a system with several time scales as mentioned above. If the length of the simulation is smaller than the smallest correlation time, the walker may become trapped in an effectively disconnected region and the sampling of phase space is incomplete. By increasing the time, the memory of the initial condition in the sampling decreases as the walker crosses to other previously unreachable regions. These oscillations and hoppings can be modeled by a superposition of stochastic processes with different correlation times. These processes with non-zero correlation times are known as colored noise processes (as opposed to zero correlation time white noise processes).  From the study of the autocorrelation functions of a stochastic model defined using these colored noise processes, we can verify that, at a fixed run length $`K`$, there exist two different groups of processes; those that contribute to the autocorrelation function with terms that decay like $`1/k`$ (called diffusive processes), and those that contribute to the autocorrelation function with terms that decay slower than $`1/k`$ (called non-diffusive processes). When the time of the simulation is increased, some non-diffusive processes at shorter run lengths, start to contribute to the autocorrelation function like diffusive processes. After the walk length reaches the largest correlation time $`\tau _\mathrm{\Lambda }`$, all processes contribute to the autocorrelation function with terms that decay like $`1/k`$. At this point, the simulation is at the diffusive regime and effective ergodicity has been reached. A principal goal of this work is to investigate the way in which the MC output $`\left\{U_k\right\}`$ reaches the diffusive limit (i.e. the ergodic limit) by studying the properties of autocorrelation functions under changes of scale in time, $`KbK`$ with $`b>1`$. By time scaling it is possible to infer the decay law of the non-diffusive contributions with respect to the total simulation time $`K`$. The functional dependence of the non-diffusive contributions on the parameter $`b`$ that is used to scale $`K`$ is determined empirically. We have found the decay law so determined to be a particularly valuable method of concluding when a simulation can be considered ergodic. Unlike previous studies  that only have investigated the behavior of certain autocorrelation functions in the ergodic regime, by focusing on the approach to ergodic behavior we have a more careful monitor of the onset of ergodicity. Once the non-diffusive contributions have decayed to a point where they are too small to be distinguished from zero to within the fluctuations of the calculation, we can say that the ergodic limit has been reached.

The autocorrelation functions we use to measure the approach to the ergodic limit are based on one of the probes of ergodicity developed by Thirumalai and co-workers , and is often called the energy metric. The energy metric has been proposed as an alternative to other techniques  (like the study of the Lyapunov exponents ) for the study of ergodic properties in molecular dynamics (MD) simulations. The metric has been used to study the relative efficiency of MC simulation methods as well.  The MC metric as used in the current work can easily be extended from the energy to other scalar observables of the system.

We present two key issues in this paper. First, from the knowledge of the decay law of the non-diffusive contributions to the MC metric, we infer how long a simulation must be to be considered effectively ergodic. Second, once the ergodic limit is reached, we can compare the results from different numerical algorithms to measure relative efficiencies. Because the outcomes of MC simulations are noisy, we have found it useful to separate diffusive and non-diffusive terms in the MC metric with a Fourier analysis so that we can neglect the high frequency components of the noise. This technique has given reproducible results.

To test the match between the stochastic model and actual Monte Carlo simulations, we examine the approach to ergodic behavior in simulations of Lennard-Jones clusters. Recently  we have studied the thermodynamic properties of Lennard-Jones clusters as a function of temperature using both J-walking  and parallel tempering methods.  Both simulation techniques require an initial high temperature that must be ergodic when Metropolis Monte Carlo methods  are used. If the Metropolis method does not give ergodic results at the initial high temperature, systematic errors propagate to the lower temperatures in J-walking and parallel tempering simulations, and the results can be flawed or meaningless. In most Monte Carlo simulations of clusters at finite temperatures,  the clusters are defined by enclosing the atoms within a constraining potential about the center of mass of the system. The constraining potential is necessary because clusters at finite temperatures have finite vapor pressures, and the association of any one atom with the cluster can be ill-defined. From experience  we have found that if the radius of the constraining potential and the initial high temperature are not both carefully chosen, it can be difficult to attain ergodicity with Metropolis methods. A key concern then is the choice of constraining radius and the choice of initial temperature. We verify the stochastic model by investigating Monte Carlo simulation results as a function of the temperature and the size of the constraining potential.

The contents of the remainder of this paper are as follows. In Section II we motivate the studies that follow by examining numerally the behavior of a set of Monte Carlo simulations of a 13-particle Lennard-Jones cluster. This cluster system is used to illustrate the results throughout this paper. In Section III we introduce the stochastic model based on a continuous time sequence. In Section IV we extend the model to discrete time sequences characteristic of actual Monte Carlo simulations. In Section V we test the discrete stochastic model with applications to Lennard-Jones clusters and in Section VI we summarize our conclusions. Many of the key derivations needed for the developments are found in two appendices.

## II An Example Calculation

Before discussing the major developments of this work, it is useful to understand the nature of the problem we are attempting to solve by examining some numerical results on a prototypical system. We take the 13-particle Lennard-Jones cluster defined by the potential function

$$V(𝐱)=4\epsilon \underset{i=2}{\overset{N}{}}\underset{j=1}{\overset{i1}{}}\left[\left(\frac{\sigma }{r_{ij}}\right)^{12}\left(\frac{\sigma }{r_{ij}}\right)^6\right]+\underset{i=1}{\overset{N}{}}V_C(\stackrel{}{x}_i,R_c)$$
(2)
where $`\epsilon `$ and $`\sigma `$ are the standard Lennard-Jones energy and length parameters, $`N`$ is the number of particles in the cluster (13 in the present case), $`r_{ij}`$ is the distance between particles $`i`$ and $`j`$

$$r_{ij}=|\stackrel{}{x}_i\stackrel{}{x}_j|,$$
(3)
and $`V_C`$ is the constraining potential discussed in Sec. I

$$V_C(\stackrel{}{x}_i,R_c)=\{\begin{array}{cc}0\hfill & \left|\stackrel{}{x}_i\stackrel{}{X}_c\right|<R_c\hfill \\ \mathrm{}\hfill & R_c<\left|\stackrel{}{x}_i\stackrel{}{X}_c\right|\hfill \end{array},$$
(4)
where $`\stackrel{}{X}_c`$ is the coordinate of the center of mass of the cluster and $`R_c`$ is the radius of the constraining sphere. The 13-particle Lennard-Jones cluster has a complex potential surface with many minima separated by significant energy barriers,  and ergodicity problems associated with the simulation of properties of this system are well-known.  We now consider a Metropolis MC simulation of the average potential energy of the system in the canonical ensemble at temperature $`k_BT/\epsilon =0.393`$($`k_B`$ is the Boltzmann constant). This average potential energy $`\overline{V}_k`$ is defined by

$$\overline{V}_k=\frac{1}{k}\underset{k^{}=1}{\overset{k}{}}V_k^{}$$
(5)
and is displayed in the upper panel of Fig. 1 as a function of the walk length $`k`$ for 20 independent simulations each initialized from a random configuration. Over the maximum time scale $`K`$ of the walks, it apparent that the potential energy averaged over each independent walk has not converged to the same result. Such unreproducible behavior is symptomatic of a simulation not yet at the ergodic limit.

Fig. 1: The upper panel shows the “time evolution” of $`\overline{V}_k`$ (in units of $`\epsilon `$) for $`M=20`$ independent experiments. The lower panel shows $`d_k`$ (in units of $`\epsilon ^2`$) vs. $`k`$ for the experiments of the upper panel. $`R_c`$ has been set to $`4\sigma `$ and $`k_BT/\epsilon =0.393`$. At least two basins with different energies are present. Clearly, $`d_k`$ goes to a constant when $`k`$ is increased within the total time scale of the simulation.

At the ergodic limit (i.e. for the maximum walk length $`K`$ greater than that included in Fig. 1) the averages displayed in the upper panel of Fig. 1 must approach the same value for each walker. Using related ideas developed elsewhere,  the extent to which the walks approach the same limit can be measured in terms of a metric $`d_k`$ defined by

$$d_k=\frac{2}{M(M1)}\underset{i=2}{\overset{M}{}}\underset{j=1}{\overset{i1}{}}\left[\overline{V}_k^{(i)}\overline{V}_k^{(j)}\right]^2,$$
(6)
In Eq. (6) $`M`$ represents the number of independent walks, and $`\overline{V}_k^{(i)}`$ is the average potential energy computed in walk $`i`$ at MC time $`k`$. The metric measures the energy fluctuations in the walk as a function of the walk length. For an ergodic simulation, the metric must decay to zero. For the 20 simulations of the 13-particle Lennard-Jones cluster, the metric as a function of $`k`$ is plotted in the lower panel of Fig. 1. Rather than asymptotically approaching zero, over the short length of the walk displayed here, $`d_k`$ has decayed to a constant, and as discussed later in this paper, over the time scale of this simulation, $`d_k`$ can be qualitatively represented by the function

$$d_k=\frac{A_K}{k}+B_K$$
(7)
where $`A_K`$ and $`B_K`$ are coefficients that are dependent on the total walk length $`K`$. As $`K`$ is increased to a time where the walk is ergodic, $`B_K`$ must decay to zero. Major goals of this work are to understand how $`B_K`$ decays and to use the discovered decay law to determine the onset of ergodic behavior. Our approach is to introduce first a continuous stochastic model of a simulation followed by a discrete model more clearly linked to actual MC studies.

## III Stochastic Model

We have discussed in the introduction how the output of MC simulations can be considered to be a combination of stochastic processes with different time scales, and how the contributions to autocorrelation function from these processes can vary when the length of the simulation is enlarged. Here we present a continuous time model for the stochastic processes that occur in a simulation. Even though a MC simulation occurs in a discrete time (each MC point represents a time unit), we find that the continuous model helps to understand better the ideas used in the modeling of the MC output.

In this section the ensemble mean value is used to find the expression for the autocorrelation functions of the model. Although in actual numerical calculations the ensemble mean is replaced by a mean over a finite number of independent experiments, the results obtained here give information about the limit of an infinite sample.

The stationary process used to sample space is a stochastic process. We assume the output of the MC simulation can be modeled by a linear superposition of stochastic processes with different correlation times $`\tau _{\mathrm{}}0`$,

$$U(t)=U_c+\sqrt{\mathrm{\Gamma }_0}\xi (t)+\underset{\mathrm{}=1}{\overset{\mathrm{\Lambda }}{}}\sqrt{\mathrm{\Gamma }_{\mathrm{}}}g_{\mathrm{}}(t/\tau _{\mathrm{}}),$$
(8)
where $`U_c`$ is a constant, the random variable $`\xi (t)`$ represents white noise processes with zero correlation time ($`\tau _0=0`$), and the $`\{g_{\mathrm{}}(t/\tau _{\mathrm{}})\}`$ are stochastic processes with correlation times $`\tau _{\mathrm{}}>0`$. $`\xi (t)`$ and $`g_{\mathrm{}}(t/\tau _{\mathrm{}})`$ have units of the square root of time, and $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_{\mathrm{}}`$ are constants with units of $`U^2/t`$. If $`U`$ is chosen to be the the $`x`$-coordinate of a particle, $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_{\mathrm{}}`$ have units of a diffusion constant. Consequently we refer to these constants as generalized diffusion coefficients. The white noise process has the following properties

$`\xi (t)`$ $`=`$ $`0`$ (9)
$`\xi (t)\xi (t^{})`$ $`=`$ $`\delta (tt^{}),`$ (10)
and the remaining colored noise processes are assumed to satisfy

$`g_{\mathrm{}}(t/\tau _{\mathrm{}})`$ $`=`$ $`0`$ (11)
$`g_{\mathrm{}}(t/\tau _{\mathrm{}})g_{\mathrm{}}(t^{}/\tau _{\mathrm{}})`$ $`=`$ $`{\displaystyle \frac{1}{\tau _{\mathrm{}}}}f_{\mathrm{}}\left({\displaystyle \frac{|tt^{}|}{\tau _{\mathrm{}}}}\right)`$ (12)
so that they represent processes with a memory $`f_{\mathrm{}}`$. Even though correlations between processes with different correlation times may be non-zero, we assume the processes to be independent, i.e.

$`g_{\mathrm{}}(t/\tau _{\mathrm{}})g_{\mathrm{}^{}}(t^{}/\tau _{\mathrm{}^{}})`$ $`=`$ $`g_{\mathrm{}}(t/\tau _{\mathrm{}})g_{\mathrm{}^{}}(t^{}/\tau _{\mathrm{}^{}})`$ (13)
$`=`$ $`0\mathrm{}\mathrm{}^{}`$ (14)
$`\xi (t)g_{\mathrm{}}(t^{}/\tau _{\mathrm{}})`$ $`=`$ $`\xi (t)g_{\mathrm{}}(t^{}/\tau _{\mathrm{}})`$ (15)
$`=`$ $`0t\mathrm{and}t^{}.`$ (16)
The memory function is assumed to be a continuous function that depends only on the distance between $`t`$ and $`t^{}`$ disregarding the time origin (stationary condition). The memory function represents the correlation between two times of the process $`g_{\mathrm{}}`$. In our model we impose the condition

$$\frac{t}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{t}{\tau _{\mathrm{}}}\right)<_0^t𝑑t^{}\frac{1}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{t^{}}{\tau _{\mathrm{}}}\right)<\mathrm{}.$$
(17)
The scope and implications of the leftmost inequality are explored in Appendix A. In Appendix A we also examine the conditions $`f_{\mathrm{}}`$ must satisfy in order to yield contributions to autocorrelation function that decay more weakly than $`1/t`$. We now assume that this inequality can be taken as a bound to possible maxima of $`f_{\mathrm{}}`$ appearing at $`t>0`$. The rightmost inequality enables us to assume $`f_{\mathrm{}}`$ is normalized

$$_{\mathrm{}}^{\mathrm{}}𝑑t\frac{1}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{|t|}{\tau _{\mathrm{}}}\right)=1.$$
(18)
We have identified here the time scale $`\tau _{\mathrm{}}`$ with the correlation time of the stochastic process $`g_{\mathrm{}}`$. This identification is valid if

$$_{\mathrm{}}^{\mathrm{}}𝑑t\frac{|t|}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{|t|}{\tau _{\mathrm{}}}\right)=\tau _{\mathrm{}},$$
(19)
which implies that the behavior of $`f_{\mathrm{}}`$ at large $`t`$ must be $`𝒪\left(t^{(2+ϵ)}\right)`$, or smaller.

In addition, by the properties of the ensemble mean value, we have that for all real $`\lambda `$

$`0`$ $``$ $`\left[g_{\mathrm{}}(t/\tau _{\mathrm{}})+\lambda g_{\mathrm{}}(t^{}/\tau _{\mathrm{}})\right]^2`$ (20)
$``$ $`g_{\mathrm{}}(t/\tau _{\mathrm{}})^2+2\lambda g_{\mathrm{}}(t/\tau _{\mathrm{}})g_{\mathrm{}}(t^{}/\tau _{\mathrm{}})+\lambda ^2g_{\mathrm{}}(t^{}/\tau _{\mathrm{}})^2`$ (21)
$``$ $`{\displaystyle \frac{1}{\tau _{\mathrm{}}}}\left\{f_{\mathrm{}}(0)+2\lambda f_{\mathrm{}}\left({\displaystyle \frac{|tt^{}|}{\tau _{\mathrm{}}}}\right)+\lambda ^2f_{\mathrm{}}(0)\right\}.`$ (22)
Equation (22) must be true for all $`\lambda `$. Therefore, the discriminant of the polynomial in $`\lambda `$ must be non-positive

$$4\left[f_{\mathrm{}}\left(\frac{|tt^{}|}{\tau _{\mathrm{}}}\right)f_{\mathrm{}}(0)\right]\left[f_{\mathrm{}}\left(\frac{|tt^{}|}{\tau _{\mathrm{}}}\right)+f_{\mathrm{}}(0)\right]0.$$
(23)
Consequently, $`f_{\mathrm{}}(0)=\mathrm{max}\{f_{\mathrm{}}(x)x0\}`$. Other properties of $`f_{\mathrm{}}`$ are studied in Appendix A.

The ensemble mean value $`U`$ is time independent. The ensemble mean value of the noise processes is zero. Therefore, $`U_c`$ must be equal to $`U`$. Processes defined by Eq. (8) have two different components, uncorrelated white noise and correlated processes with correlation time $`\tau _{\mathrm{}}`$. Because the goal of the simulation is the calculation of the ensemble mean $`U`$ by the analysis of the time series, we study the behavior of the temporal mean $`\overline{U}(t)`$

$`\overline{U}(t)`$ $`=`$ $`{\displaystyle \frac{1}{t}}{\displaystyle _0^t}𝑑t^{}U(t^{})`$ (24)
$`=`$ $`U+{\displaystyle \frac{1}{t}}\sqrt{\mathrm{\Gamma }_0}W(t)+{\displaystyle \frac{1}{t}}{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{\Lambda }}{}}}\sqrt{\mathrm{\Gamma }_{\mathrm{}}}G_{\mathrm{}}(t/\tau _{\mathrm{}}),`$ (25)
where $`W(t)`$ is a Wiener process,

$`W(t)`$ $`=`$ $`{\displaystyle _0^t}𝑑t^{}\xi (t^{})`$ (26)
$`W(t)`$ $`=`$ $`0`$ (27)
$`W(t)W(t^{})`$ $`=`$ $`t_<,`$ (28)
with $`t_<=\mathrm{min}(t,t^{})`$, and $`G_{\mathrm{}}(t/\tau _{\mathrm{}})=_0^t𝑑t^{}g_{\mathrm{}}(t^{}/\tau _{\mathrm{}})`$.

Equation (25) implies that the evolution of the temporal mean $`\overline{U}(t)`$ has the same structure as $`U`$, with an uncorrelated term and terms with tailed correlation functions.

The autocorrelation function of the process $`\overline{U}`$ at times $`t`$ and $`t^{}`$ is defined by

$`\kappa (t,t^{})`$ $`=`$ $`\left(\overline{U}(t)U\right)\left(\overline{U}(t^{})U\right)`$ (29)
$`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_0}{tt^{}}}W(t)W(t^{})+{\displaystyle \frac{1}{tt^{}}}{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{\Lambda }}{}}}\mathrm{\Gamma }_{\mathrm{}}G_{\mathrm{}}(t/\tau _{\mathrm{}})G_{\mathrm{}}(t^{}/\tau _{\mathrm{}}),`$ (30)
where we have used Eqs. (14) and (16) to neglect terms involving processes with different correlation times.

Because we have assumed $`f_{\mathrm{}}`$ is a continuous function, $`f_{\mathrm{}}`$ reaches its maximum and minimum value within any closed interval considered. The $`\mathrm{}`$th non-diffusive contribution to $`\kappa (t,t^{})`$

$$\frac{1}{tt^{}}G_{\mathrm{}}(t/\tau _{\mathrm{}})G_{\mathrm{}}(t^{}/\tau _{\mathrm{}})=\frac{1}{tt^{}}_0^t𝑑t_1_0^t^{}𝑑t_2\frac{1}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{|t_1t_2|}{\tau _{\mathrm{}}}\right),$$
(31)
is bounded

$`{\displaystyle \frac{1}{t_<t_>}}{\displaystyle _0^{t_<}}𝑑t_1{\displaystyle _0^{t_>}}𝑑t_2{\displaystyle \frac{1}{\tau _{\mathrm{}}}}f_{\mathrm{}}\left({\displaystyle \frac{t_{min}}{\tau _{\mathrm{}}}}\right){\displaystyle \frac{1}{tt^{}}}`$ $`G_{\mathrm{}}(t/\tau _{\mathrm{}})G_{\mathrm{}}(t^{}/\tau _{\mathrm{}})`$ $`{\displaystyle \frac{1}{t_<t_>}}{\displaystyle _0^{t_<}}𝑑t_1{\displaystyle _0^{t_>}}𝑑t_2{\displaystyle \frac{1}{\tau _{\mathrm{}}}}f_{\mathrm{}}(0)`$ (32)
$`{\displaystyle \frac{1}{\tau _{\mathrm{}}}}f_{\mathrm{}}\left({\displaystyle \frac{t_{min}}{\tau _{\mathrm{}}}}\right){\displaystyle \frac{1}{tt^{}}}`$ $`G_{\mathrm{}}(t/\tau _{\mathrm{}})G_{\mathrm{}}(t^{}/\tau _{\mathrm{}})`$ $`{\displaystyle \frac{1}{\tau _{\mathrm{}}}}f_{\mathrm{}}(0),`$ (33)
where $`t_>=\mathrm{max}(t,t^{})`$, and $`t_{min}`$ is the time at which $`f_{\mathrm{}}`$ reaches its minimum value in the closed interval $`[0,t_>]`$. There exists a $`t_{\mathrm{}}^{}(t_>)[0,t_{min}]`$  such that,

$$\frac{1}{tt^{}}G_{\mathrm{}}(t/\tau _{\mathrm{}})G_{\mathrm{}}(t^{}/\tau _{\mathrm{}})=\frac{1}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{t_{\mathrm{}}^{}(t_>)}{\tau _{\mathrm{}}}\right).$$
(34)
Using Eqs. (28) and (34) in (30), we find that

$$\kappa (t,t^{})=\frac{\mathrm{\Gamma }_0}{t_>}+\underset{\mathrm{}=1}{\overset{\mathrm{\Lambda }}{}}\frac{\mathrm{\Gamma }_{\mathrm{}}}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{t_{\mathrm{}}^{}(t_>)}{\tau _{\mathrm{}}}\right).$$
(35)
For all times shorter than $`\tau _1`$ the autocorrelation function is the sum of diffusive contributions (proportional to $`1/t`$) plus non-diffusive contributions. These contributions implicitly depend on $`t_>`$ through $`t_{\mathrm{}}^{}(t_>)`$. We assume that $`f_{\mathrm{}}`$ satisfies the conditions stated in Appendix A, so that the dependence of $`f_{\mathrm{}}`$ on $`t`$ is weaker than $`1/t`$ (for total time scales shorter than $`\tau _{\mathrm{}}`$; see Appendix A).

We next consider the behavior of Eq. (35) for time scales greater than $`\tau _1`$. Under the scale change $`tbt`$ such that $`\tau _1bt_>\tau _2`$, the contributions to the correlation function from the process with correlation time $`\tau _1`$ can be considered diffusive \[in other words, by virtue of Eqs. (10) and (12), $`f_1/\tau _1`$ has become a delta function\]. With $`bt_>\tau _2`$, the other processes preserve their old properties. Then, the autocorrelation function can be expressed

$$\kappa (bt,bt^{})=\frac{\mathrm{\Gamma }_0+\mathrm{\Gamma }_1}{bt_>}+\underset{\mathrm{}=2}{\overset{\mathrm{\Lambda }}{}}\frac{\mathrm{\Gamma }_{\mathrm{}}}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{t_b\mathrm{}^{}(t_>)}{\tau _b\mathrm{}}\right).$$
(36)
The complete derivation of Eq. (36) can be found in Appendix B. For a times larger than the correlation time $`\tau _\mathrm{\Lambda }`$, all contributions to the autocorrelation function are diffusive, the simulation can be considered ergodic, the sampling complete, and the temporal mean is equal to the ensemble mean within $`𝒪(1/t)`$ mean square fluctuations.

## IV Discrete time sequences and the MC metric

Monte Carlo simulations generate discrete sequences $`U_k`$ of values of the quantity under study. Additionally, in actual calculations the ensemble of sequences is represented by a finite rather than an infinite set. In this section, the model developed in the previous section is extended to finite sets of discrete sequences. We express the $`M`$ sequences $`\left\{U_k^{(m)}\right\}_{k=1}^K`$, where the label $`(m)`$ ranges from 1 to $`M`$. The exact ensemble mean value $`U`$ can be obtained in the limit that $`M`$ becomes infinite. In analogy with the model developed in Section III, each output is assumed to have the form

$$U_k^{(m)}=U+\sqrt{\mathrm{\Gamma }_0}\xi _k^{(m)}+\underset{\mathrm{}=1}{\overset{\mathrm{\Lambda }}{}}\sqrt{\mathrm{\Gamma }_{\mathrm{}}}g_{\mathrm{};k/\tau _{\mathrm{}}}^{(m)},$$
(37)
where

$`\xi _k^{(m)}`$ $`=`$ $`0`$ (38)
$`\xi _k^{(m)}\xi _k^{}^{(n)}`$ $`=`$ $`\delta _{m,n}\delta _{k,k^{}}`$ (39)
$`g_{\mathrm{};k/\tau _{\mathrm{}}}^{(m)}`$ $`=`$ $`0`$ (40)
$`g_{\mathrm{};k/\tau _{\mathrm{}}}^{(m)}g_{\mathrm{}^{};k^{}/\tau _{\mathrm{}^{}}}^{(n)}`$ $`=`$ $`\delta _{m,n}\delta _\mathrm{},\mathrm{}^{}f_{\mathrm{}}\left({\displaystyle \frac{|kk^{}|}{\tau _{\mathrm{}}}}\right).`$ (41)
The true ensemble average $`U`$ does not depend on the index $`m`$.

In the discrete case we define a metric

$$d_k=\frac{2}{M(M1)}\underset{i=2}{\overset{M}{}}\underset{j=1}{\overset{i1}{}}\left[\overline{U}_k^{(i)}\overline{U}_k^{(j)}\right]^2,$$
(42)
where the bars represent the temporal mean value

$`\overline{U}_k^{(m)}`$ $`=`$ $`{\displaystyle \frac{1}{k}}{\displaystyle \underset{k^{}=1}{\overset{k}{}}}U_k^{}^{(m)}`$ (43)
$`=`$ $`U+{\displaystyle \frac{\sqrt{\mathrm{\Gamma }_0}}{k}}W_k^{(m)}+{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{\Lambda }}{}}}{\displaystyle \frac{\sqrt{\mathrm{\Gamma }_{\mathrm{}}}}{k}}G_{\mathrm{};k/\tau _{\mathrm{}}}^{(m)},`$ (44)
with

$`W_k^{(m)}`$ $`=`$ $`{\displaystyle \underset{k^{}=1}{\overset{k}{}}}\xi _k^{}^{(m)}`$ (45)
$`G_{\mathrm{};k/\tau _{\mathrm{}}}^{(m)}`$ $`=`$ $`{\displaystyle \underset{k^{}=1}{\overset{k}{}}}g_{\mathrm{};k^{}/\tau _{\mathrm{}}}^{(m)}.`$ (46)
Observe that in the present case, our finite sample of the infinite ensemble is the set of outcomes from $`M`$ independent numerical experiments. The metric we have defined in Eq. (42) can be contrasted with alternative metrics  previously defined for molecular dynamics simulations. These alternative metrics examine the fluctuations of two simulations initialized from different components of configuration space averaged with respect to all the particles in the system. The metric we use in this work is determined using an average with respect to $`M`$ independent simulations that represent a subset of the full ensemble.

Using the model presented in Eq. (37), we now develop a way to predict the behavior of the MC simulation in the non-ergodic and the ergodic regimes. We first consider the case that the total simulation time $`K`$ is larger than the first correlation time $`\tau _1`$ but shorter than $`\tau _2`$, i.e. $`\tau _1K\tau _2`$. The expression for $`d_k`$ is given by

$`d_k`$ $`=`$ $`{\displaystyle \frac{2}{M(M1)}}{\displaystyle \underset{i=2}{\overset{M}{}}}{\displaystyle \underset{j=1}{\overset{i1}{}}}\left[\left(\overline{U}_k^{(i)}U\right)\left(\overline{U}_k^{(j)}U\right)\right]^2`$ (47)
$`=`$ $`{\displaystyle \frac{2}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}\left(\overline{U}_k^{(i)}U\right)^2{\displaystyle \frac{4}{M(M1)}}{\displaystyle \underset{i=2}{\overset{M}{}}}{\displaystyle \underset{j=1}{\overset{i1}{}}}\left(\overline{U}_k^{(i)}U\right)\left(\overline{U}_k^{(j)}U\right)`$ (48)
$`=`$ $`2\mathrm{\Gamma }_0{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}\left({\displaystyle \frac{W_k^{(i)}}{k}}\right)^2+2\mathrm{\Gamma }_1{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}\left({\displaystyle \frac{G_{1;k/\tau _1}^{(i)}}{k}}\right)^2+2{\displaystyle \underset{\mathrm{}=2}{\overset{\mathrm{\Lambda }}{}}}\mathrm{\Gamma }_{\mathrm{}}{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}\left({\displaystyle \frac{G_{\mathrm{};k/\tau _{\mathrm{}}}^{(i)}}{k}}\right)^2`$ (51)
$`+\mathrm{\hspace{0.17em}4}{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{\Lambda }}{}}}\sqrt{\mathrm{\Gamma }_0\mathrm{\Gamma }_{\mathrm{}}}{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \frac{W_k^{(i)}G_{\mathrm{};k/\tau _{\mathrm{}}}^{(i)}}{k^2}}+\mathrm{\hspace{0.17em}4}{\displaystyle \underset{\mathrm{}=2}{\overset{\mathrm{\Lambda }}{}}}{\displaystyle \underset{\mathrm{}^{}=1}{\overset{\mathrm{}1}{}}}\sqrt{\mathrm{\Gamma }_{\mathrm{}}\mathrm{\Gamma }_{\mathrm{}^{}}}{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \frac{G_{\mathrm{};k/\tau _{\mathrm{}}}^{(i)}G_{\mathrm{}^{};k/\tau _{\mathrm{}^{}}}^{(i)}}{k^2}}`$
$`{\displaystyle \frac{4}{M(M1)}}{\displaystyle \underset{i=2}{\overset{M}{}}}{\displaystyle \underset{j=1}{\overset{i1}{}}}\left(\overline{U}_k^{(i)}U\right)\left(\overline{U}_k^{(j)}U\right).`$
If the number of experiments $`M`$ is sufficiently large, we can neglect terms involving processes with different correlation times, and products of sequences belonging to different experiments. Under these assumptions we obtain

$$d_k=2\frac{\mathrm{\Gamma }_0}{k}\frac{1}{M}\underset{i=1}{\overset{M}{}}\frac{W_{k}^{(i)}{}_{}{}^{2}}{k}+2\frac{\mathrm{\Gamma }_1}{k}\frac{1}{M}\underset{i=1}{\overset{M}{}}\frac{G_{1;k/\tau _1}^{(i)}{}_{}{}^{2}}{k}+2\underset{\mathrm{}=2}{\overset{\mathrm{\Lambda }}{}}\mathrm{\Gamma }_{\mathrm{}}\frac{1}{M}\underset{i=1}{\overset{M}{}}\left(\frac{G_{\mathrm{};k/\tau _{\mathrm{}}}^{(i)}}{k}\right)^2.$$
(52)
Equation (52) preserves the form of Eq. (35). To make this statement explicit, let us rewrite Eq. (52) as

$$d_k=2\frac{\mathrm{\Gamma }_k}{k}+2\mathrm{{\rm Y}}_k,$$
(53)
where

$`\mathrm{\Gamma }_k`$ $`=`$ $`\mathrm{\Gamma }_0{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \frac{W_{k}^{(i)}{}_{}{}^{2}}{k}}+\mathrm{\Gamma }_1{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \frac{G_{1;k/\tau _1}^{(i)}{}_{}{}^{2}}{k}}`$ (54)
$`\mathrm{{\rm Y}}_k`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{\overset{\mathrm{\Lambda }}{}}}\mathrm{\Gamma }_{\mathrm{}}{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}\left({\displaystyle \frac{G_{\mathrm{};k/\tau _{\mathrm{}}}^{(i)}}{k}}\right)^2.`$ (55)
In Appendix B we present a study of the way non-diffusive contribution become diffusive under time scale changes. If $`M`$ is sufficiently large and $`\tau _1K\tau _2`$, by virtue of Appendix B, $`\mathrm{\Gamma }_k`$ must roughly be a constant. By roughly a constant we mean a constant $`C`$ plus some rapidly fluctuating function $`\zeta _k`$, with the following properties: a) $`\zeta _k=0`$ and b) $`|C|\mathrm{max}_{k=1,2,\mathrm{},K}(|\zeta _k|)`$. Then

$$\mathrm{\Gamma }_k\mathrm{\Gamma }_K+\zeta _k$$
(56)
If $`K`$ is enlarged, we expect to have a larger value of $`\mathrm{\Gamma }_K`$. $`\mathrm{{\rm Y}}_k`$ is a quantity related to the memory functions $`f_{\mathrm{}}`$ with correlation times $`\tau _{\mathrm{}}K`$. In the continuous time model, the colored noise processes contribute to the autocorrelation function with terms proportional to $`f_{\mathrm{}}(t_{\mathrm{}}^{}(t_>)/\tau _{\mathrm{}})`$, which are weakly dependent on $`t`$ (see Appendix A). We can expect $`\mathrm{{\rm Y}}_k`$ to be weakly dependent on $`k`$, and for sequences of length $`K`$ and for $`M`$ sufficiently large, we consider this quantity roughly to be a constant

$$\mathrm{{\rm Y}}_k\mathrm{{\rm Y}}_K+\beta _k.$$
(57)
where $`\beta _k`$ represents additional random noise. Then, for a given length $`kK`$, the MC metric $`d_k`$ can be approximated by

$$d_k=2\frac{\mathrm{\Gamma }_K}{k}+2\mathrm{{\rm Y}}_K+\gamma _k,$$
(58)
where $`\gamma _k=2(\zeta _k/k+\beta _k)`$ represents remaining stochastic noise from both contributions. In this approximation, $`\mathrm{\Gamma }_K`$ and $`\mathrm{{\rm Y}}_K`$ are the quantities that carry the long time dependence. Short time features appear in the $`1/k`$ dependence and in the remaining noise $`\gamma _k`$. If the sequences considered are increased in size by a factor of $`b`$, such that $`\tau _{\lambda 1}K\tau _\lambda bK`$ for a given $`1\lambda \mathrm{\Lambda }`$, $`\mathrm{\Gamma }_K`$ ($`\mathrm{{\rm Y}}_K`$) is increased (decreased) (see Appendix B). Then,

$$d_{bk}=2\frac{\mathrm{\Gamma }_{bK}}{bk}+2\mathrm{{\rm Y}}_{bK}+\gamma _{bk},$$
(59)
where $`\mathrm{{\rm Y}}_{bK}`$ must go to zero and $`\mathrm{\Gamma }_{bK}`$ must approach a constant when $`b`$ is increased. By virtue of the expected behavior of the non-diffusive contributions (see Appendix A), we propose the following expression for $`\mathrm{{\rm Y}}_{bK}`$

$$\mathrm{{\rm Y}}_{bK}=\mathrm{{\rm Y}}_K\varphi (b),$$
(60)
where $`\varphi (b)`$ is a decreasing function of $`b`$. Moreover, $`\mathrm{{\rm Y}}_{bK}`$ is a sum of non-diffusive contributions. As presented in Appendix A, each non-diffusive contribution to the autocorrelation function has a relative variation smaller than the relative variation of the diffusive contribution, namely $`11/b`$. If this inequality is applicable to the sum of non-diffusive contributions, we have that

$`1{\displaystyle \frac{1}{b}}`$ $`>`$ $`1{\displaystyle \frac{\mathrm{{\rm Y}}_{bK}}{\mathrm{{\rm Y}}_K}}`$ (61)
$`1{\displaystyle \frac{1}{b}}`$ $`>`$ $`1\varphi (b)`$ (62)
$`1`$ $`<`$ $`b\varphi (b),`$ (63)
for all $`b>1`$. Then, $`\varphi `$ must be either

$$\varphi (b)=b^\upsilon $$
(64)
or

$$\varphi (b)=\frac{1}{\eta \mathrm{ln}(b)+1},$$
(65)
with $`0<\upsilon <1`$ and $`0<\eta 1`$. Equation (65) can be thought as the limit of Eq. (64) when the exponent goes to zero. We know of no a priori argument to justify Eq. (60). However, as is discussed in Section V, our numerical experience has shown Eq. (60) to be obeyed in all cases we have examined.

Our goal is to develop a criterion to decide when the simulation can be considered ergodic. From the previous considerations it is clear that the ergodic limit is reached when $`\mathrm{{\rm Y}}_K`$ is indistinguishable from zero. The output from a MC simulation is usually noisy. Therefore, $`\gamma _k`$ can not be neglected. A useful way to separate diffusive and non-diffusive contributions and to eliminate the stochastic noise from Eq. (58), is to perform a Fourier analysis of the function $`kd_k`$. Let us define the frequencies $`\omega _n=(2\pi /K)n`$, with $`n=0,1,\mathrm{},K1`$. The discrete Fourier transform of the function $`kd_k`$ is the signal $`Y_K(\omega _n)`$

$`Y_K(\omega _n)`$ $`=`$ $`\widehat{kd_k}(\omega _n)={\displaystyle \frac{1}{K}}{\displaystyle \underset{k=1}{\overset{K}{}}}\mathrm{exp}(i\omega _nk)kd_k`$ (66)
$`=`$ $`{\displaystyle \frac{2}{K}}{\displaystyle \underset{k=1}{\overset{K}{}}}\mathrm{exp}(i\omega _nk)\mathrm{\Gamma }_K+{\displaystyle \frac{2}{K}}{\displaystyle \underset{k=1}{\overset{K}{}}}k\mathrm{exp}(i\omega _nk)\mathrm{{\rm Y}}_K+\widehat{k\gamma _k}(\omega _n)`$ (67)
$`=`$ $`2\delta _{n,0}\mathrm{\Gamma }_K+\left\{\delta _{n,0}(K+1)+\left(1\delta _{n,0}\right)\left(1+i\mathrm{cot}(\omega _n/2)\right)\right\}\mathrm{{\rm Y}}_K+\widehat{k\gamma _k}(\omega _n)`$ (68)
$`=`$ $`2\delta _{n,0}\mathrm{\Gamma }_K+\left(K\delta _{n,0}+1\right)\mathrm{{\rm Y}}_K+i\left(1\delta _{n,0}\right)\mathrm{cot}(\omega _n/2)\mathrm{{\rm Y}}_K+\widehat{k\gamma _k}(\omega _n)`$ (69)
In general, $`\widehat{k\gamma _k}(\omega _n)`$ is negligible except at high frequencies. For small positive values of the frequency we can make the approximation $`\mathrm{cot}(\omega _n/2)2/\omega _n`$. From this approximation we have

$$\mathrm{Im}\left(Y_K(\omega _n)\right)\frac{2}{\omega _n}\mathrm{{\rm Y}}_K.$$
(70)
The real part of Eq. (69) for positive frequencies is

$$\mathrm{Re}\left(Y_K(\omega _n)\right)=\mathrm{{\rm Y}}_K.$$
(71)
Even though simpler than Eq. (70), we have found Eq. (71) is more sensitive to the deviations of $`d_k`$ from the approximation Eq. (58). Therefore, the data obtained from the real part is of poorer quality than the data obtained from the imaginary part.

Equation (70) implies that for a given simulation length $`K`$, the contributions to the MC metric from the non-diffusive process can be determined from a simple relationship involving the Fourier transform of the function $`kd_k`$ at low frequencies. By increasing the length of the run $`K`$ by a factor of $`b`$, it is possible to observe the dependence of $`\mathrm{{\rm Y}}_{bK}`$ on $`bK`$.

## V Applications

The concepts developed in the previous sections are sufficiently general to be applied to any kind of MC simulation. We devote the present section to the application of the developments of this paper to the study of the Lennard-Jones 13-particle cluster in the canonical ensemble. This system has been introduced previously in Sec. II.

Some thermodynamic properties of clusters as a function of temperature exhibit rapid changes that are reminiscent of similar changes that occur for the same properties in bulk systems at phase transitions. In a bulk system a phase transition occurs at a single temperature. For clusters the rapid changes in thermodynamic properties occur over a finite temperature interval. To distinguish the temperature range where thermodynamic properties change rapidly in clusters from a true phase transition, we follow Berry et al.  and refer to such changes in physical properties as associated with a phase change. A common property that has been found to be useful in monitoring these phase change intervals of temperature is the heat capacity at constant volume

$$C_V(T)=\frac{1}{k_BT^2}(VV_T)^2_T+\frac{3}{2}Nk_B,$$
(72)
where $`_T`$ represents the classical canonical mean value.

In this work we consider the bare Metropolis (Met),  J-walking (Jw),  and parallel tempering (PT)  approaches to Monte Carlo simulations. The free variable of all these methods is the reduced temperature $`k_BT/\epsilon `$. In PT and Jw simulations, the highest temperature used ($`T_h`$) must be sufficiently large to ensure that Met is ergodic. From experience simulating a variety of systems, we have found that $`T_h`$ must also be lower than a temperature $`T_b`$ where cluster evaporation events become frequent. It is useful to think of $`T_b`$ as the cluster analogue of a boiling temperature. We have found that Met is unable to sample the boiling phase change region for clusters ergodically, using total time scales accessible to current simulations.

For the results that follow, $`U_k^{(m)}`$ is chosen to be represented by the potential energy of the system. In general $`U_k^{(m)}`$ can be any scalar property of the system. We define a pass to represent a set of single particle MC moves taken sequentially over the 13 particles in the cluster. We take $`U_k^{(m)}`$ to be the potential energy at the $`k^{}`$th pass, in the $`m^{}`$th experiment. Using Eq. (69) we can write

$$Y_K(0)=2\mathrm{\Gamma }_K+(K+1)\mathrm{{\rm Y}}_K.$$
(73)
In the non-ergodic regime, $`Y_K(0)`$ grows with $`K`$, while in the ergodic regime, the signal $`Y_K(0)`$ approaches a constant.

We begin by displaying results obtained for a calculation that has not attained ergodicity over the time scale of the simulation. We examine the 13-particle Lennard-Jones cluster with the Met algorithm setting $`R_c=4\sigma `$ at a temperature of $`k_BT_h/\epsilon =0.393`$. The temperature is chosen to be that typically used as the initial high temperature in Jw and PT studies of LJ<sub>13</sub>. By choosing a large constraining radius, the evaporation events are so frequent at the chosen temperature that attaining ergodicity proves to be quite difficult. We demonstrate the effect of reducing the constraining radius shortly.

Fig. 2: The upper panel is the signal $`Y_K(0)`$ (in units of $`\epsilon ^2`$) vs. $`K`$ for $`R_c=4\sigma `$ at $`k_BT_h/\epsilon =0.393`$. from $`M=40`$ independent experiments, of LJ<sub>13</sub>. The length of the simulation is $`10^4`$ MC passes. The lower panel shows the “time evolution” of $`\overline{U}_K`$ (in units of $`\epsilon `$) for 15 independent experiments. At least three basins with different energies are present. Clearly, the simulation at this scale of time, is not ergodic.

The number of replicas used in the calculation is $`M=40`$, and $`K=10^4`$. The upper panel of Fig. 2 shows the signal $`Y_K(0)`$ \[evaluated using Eq. (66)\], which grows along the entire simulation. This is the behavior expected in the non-ergodic regime. In the lower panel we can see the “time evolution” of the temporal mean values of 15 experiments.

Fig. 3: Upper panel: $`\mathrm{{\rm Y}}_{bK}`$ (in units of $`\epsilon ^2`$) as a function of $`\mathrm{log}_2(b)`$ for $`R_c=4\sigma `$, $`3\sigma `$, $`2.5\sigma `$, and $`2\sigma `$. For the two larger radii the full line is the best fit to the data points, according to Eq. (60) with $`\varphi `$ defined in Eq. (65). The lower panel shows the linear behavior of $`\mathrm{{\rm Y}}_K/\mathrm{{\rm Y}}_{bK}`$ vs. $`\mathrm{log}_2(b)`$, for $`R_c=4\sigma `$ and $`3\sigma `$. $`K`$ has been set to $`10^4`$.

There are three sets of curves, each of which is indicative of sampling of at least three different energy basins. At low values of $`K`$ the curves in the lower panel differ significantly. At $`K\mathrm{4\hspace{0.17em}000}`$ the high energy basin curves begin to decrease in energy. For a value of $`K`$ larger than the data displayed in Fig. 2, the curves can be expected to coalesce with the low energy basin curves. It is clear that for $`K\mathrm{10\hspace{0.17em}000}`$, the simulation is not ergodic.

In PT and Jw studies it is essential that the initial high temperature walk be ergodic. Ergodicity can be attained for LJ<sub>13</sub> by reducing the radius of the constraining potential so that evaporation events are rare. We now present a study of $`\mathrm{{\rm Y}}_K`$ as a function of $`K`$ for several values of $`R_c`$. To determine $`\mathrm{{\rm Y}}_K`$, we have calculated the Fourier transform function $`Y_K(\omega _n)`$ using Eq. (66) at a series of frequencies $`\omega _n=2\pi n/K`$ where $`n`$ has ranged from 1 to $`\mathrm{min}(\sqrt{12}bK/20\pi ,100)`$. This range of frequencies ensures the linear approximation used in Eq. (69) is valid while including sufficient numbers of points for accuracy. Using Eq. (70), we have calculated the slope of the imaginary part of $`1/Y_K(\omega _n)`$ as a function of $`\omega _n`$, for these frequencies. The data points appearing in Fig. 3 are the mean value over twenty independent calculations of the slope of $`1/Y_K(\omega _n)`$.

Fig. 4: $`\mathrm{{\rm Y}}_{bK}`$ (in units of $`\epsilon ^2`$) and its error vs. $`\mathrm{log}_2(b)`$ for $`R_c=2.5\sigma `$ and $`2\sigma `$, with $`K=10^4`$. When $`\mathrm{{\rm Y}}_{bK}`$ is on the order of its own error, the simulation can be considered ergodic. For $`R_c=2\sigma `$ the simulation becomes ergodic at $`\mathrm{log}_2(b)4`$ ($`bK16\times 10^4`$). For $`R_c=2.5\sigma `$ a longer simulation is needed to reach ergodicity.

Starting from random configurations, we have performed $`5\times 10^4`$ Met passes at $`k_BT_h/\epsilon =0.393`$. After this warmup process, we have created sequences of sizes $`bK=10^4`$, $`2\times 10^4`$, $`4\times 10^4`$, $`\mathrm{}`$, $`64\times 10^4`$. The results are presented in Fig. 3 for $`R_c=4\sigma `$, $`3\sigma `$, $`2.5\sigma `$, and $`2\sigma `$. The upper panel shows $`\mathrm{{\rm Y}}_{bK}`$ as a function of $`\mathrm{log}_2(b)`$, for fixed $`K=10^4`$. We have chosen to present the data using base 2 logarithms for clarity (each increase by 1 unit of $`\mathrm{log}_2(b)`$ represents a factor of 2 scale increase). All the data decrease with increasing $`b`$, but only $`R_c=2\sigma `$ and $`R_c=2.5\sigma `$ appear to vanish to within the error bars over the time scale of the current simulation. In the lower panel we present $`\mathrm{{\rm Y}}_K/\mathrm{{\rm Y}}_{bK}`$ as a function of $`\mathrm{log}_2(b)`$ for $`R_c=4`$ and $`3\sigma `$. The decay law suggested in Eq. (60) with $`\varphi `$ given by Eq. (65) is satisfied for both radii.

Fig. 5: The upper panel shows the decay behavior of $`\mathrm{{\rm Y}}_{bK}`$ (in units of $`\epsilon ^2`$) as a function of $`\mathrm{log}_2(b)`$ for PT and Met, at the temperature of the melting peak of the heat capacity, $`k_BT_m/\epsilon =0.282`$. From Eq. (64), we plot $`\mathrm{log}_2(\mathrm{{\rm Y}}_K/\mathrm{{\rm Y}}_{bK})`$ vs. $`\mathrm{log}_2(b)`$, to extract the value of the exponent $`\upsilon `$ (the slope of the linear fit). We have found $`\upsilon =0.93\pm 0.03`$ for PT, and $`\upsilon =0.94\pm 0.02`$ for Met. The straight lines are the best linear fits of the data points.

We have stated that the simulation can be considered effectively ergodic when $`\mathrm{{\rm Y}}_K`$ is indistinguishable from zero. In Fig. 4 we have plotted $`\mathrm{{\rm Y}}_{bK}`$ and its statistical error as a function of $`\mathrm{log}_2(b)`$ for $`R_c=2.5`$ and $`2\sigma `$. For $`R_c=2\sigma `$ the crossing point of $`\mathrm{{\rm Y}}_{bK}`$ and its error is at $`bK16\times 10^4`$. For $`R_c=2.5\sigma `$ the crossing point is at a $`bK>64\times 10^4`$. We can conclude that for $`k_BT_h/\epsilon =0.393`$ and $`R_c=2\sigma `$ the simulation can be considered effectively ergodic after $`16\times 10^4`$ Met passes.

Once a constraining radius is chosen, PT and Jw simulations require the highest temperature $`T_h`$ be chosen so that Met is ergodic. For a given $`R_c`$, the extent of ergodicity can be tested using the same metric that has been used for determining the optimum value of $`R_c`$, but by varying the temperature. For the parameters $`k_BT_h/\epsilon =0.393`$ and $`R_c=2\sigma `$ the simulation is ergodic even at very short sequence lengths. We have found that for $`k_BT_h/\epsilon <0.393`$ the simulations are not ergodic. To be sure that the parameters are appropriate, we have performed a short PT simulation ($`10^4`$ passes, ten PT passes consists of nine Met passes plus an exchange attempt) with 40 equally spaced temperatures in the range $`k_BT/\epsilon =`$\[0.028,0.393\] in order to obtain a first estimate of the position of the melting and boiling temperature regions. The boiling peak in the specific heat appears to be located at a higher temperature than $`k_BT/\epsilon =`$ 0.393. Moreover, the value of $`C_V`$ at $`k_BT/\epsilon =0.393`$ is about one-half the value of $`C_V`$ at the temperature of the melting peak $`k_BT_m/\epsilon =0.282`$. From these results we feel it is safe to choose $`R_c=2\sigma `$ and $`k_BT_h/\epsilon =0.393`$ for the calculations that follow.

We now illustrate the convergence characteristics of $`\mathrm{{\rm Y}}_K`$ when we increase the total time scale of the calculation by a factor $`b`$. We illustrate this behavior using a PT simulation of LJ<sub>13</sub>, and we focus on results at the temperature of the melting peak in the heat capacity ($`k_BT_m/\epsilon =0.282`$). We choose this temperature, because from experience  we know the statistical fluctuations are large at the melting heat capacity maximum. The large statistical fluctuations make it possible to emphasize the behavior of $`\mathrm{{\rm Y}}_K`$. We have run the PT simulation at 40 equally spaced temperatures in the range $`k_BT/\epsilon =`$ \[0.028,0.393\]. The initial warmup time has been set to $`10^4`$ Met passes, followed by $`2\times 10^4`$ PT passes. Following the warm-up period, we perform simulations of $`10^5`$, $`2\times 10^5`$, $`4\times 10^5`$, $`8\times 10^5`$, $`16\times 10^5`$, and $`32\times 10^5`$ PT passes. In each case the initial configuration has been taken to be the last configuration of the previous run. The output of the simulation are sequences of the potential energy. $`\mathrm{{\rm Y}}_K`$ has been determined in the same way as in the calculation of the high temperature parameters (presented in Fig. 3 and Fig. 4). The data points appearing in the upper panel of Fig. 5 are the mean value over twenty independent calculations of the slope of $`1/Y_K(\omega _n)`$. In the lower panel of Fig. 5 we have plotted $`\mathrm{log}_2(\mathrm{{\rm Y}}_K/\mathrm{{\rm Y}}_{bK})`$ as a function of $`\mathrm{log}_2(b)`$, where $`K=10^4`$ and $`b=1,2,4,\mathrm{},32`$. The slope of the linear fit is the exponent $`\upsilon `$, according to Eq. (64). At the temperature of the melting peak, $`\upsilon =0.93\pm 0.03`$.

It is of interest to perform a similar study of the behavior of $`\mathrm{{\rm Y}}_K`$ as a function of the time scaling for an Met calculation. We have taken the final configuration of the PT simulation at $`k_BT_m/\epsilon =0.282`$ as an initial configuration, and we have performed a simple Met simulation at that melting temperature. A graph of $`\mathrm{{\rm Y}}_{bK}`$ and $`\mathrm{log}_2(\mathrm{{\rm Y}}_K/\mathrm{{\rm Y}}_{bK})`$ as a function of $`\mathrm{log}_2(b)`$ for Met is also presented in Fig. 5. From the upper panel of Fig. 5, it is evident that Met results are not ergodic within the same scaled time as the PT results. It is also evident that the power law exponent for both Met and PT are not distinguishable. Similar studies of the power law using the Jw method also give the same exponent. Neither an increase in the number of temperatures nor changing the distribution of temperatures in both Jw and PT simulations has any effect on the calculated exponent.

Fig. 6: Comparison of the Met and Jw diffusion coefficients with the PT diffusion coefficient as a function of the reduced temperature. The dashed line represents equivalence between methods.

By using the results to compare the relative efficiencies of Met, Jw and PT simulations for the LJ<sub>13</sub> system. We have found that PT and Jw simulations can be considered ergodic if the run length is on the order of $`2\times 10^5`$ passes, while Met simulations that are initialized from configurations generated from an ergodic PT study are ergodic when the total run length consists of $`2\times 10^6`$ passes or more.

In order to compare approaches, we have calculated $`\mathrm{\Gamma }`$ as a function of the reduced temperature, for the three methods. The comparison of diffusion coefficients from different algorithms has also been used by Andricioaei and Straub . The comparison of Jw and Met with PT is presented in Fig. 6. The Jw and PT simulations are found to have comparable efficiencies using $`\mathrm{\Gamma }`$ as a measure for all calculated temperatures. At intermediate temperatures, Met is significantly less efficient. We have chosen to truncate the Jw study at $`k_BT/\epsilon =0.12`$. For temperatures below $`k_BT/\epsilon =0.12`$, Jw simulations require significant effort, because a large set of external distributions must be generated. Because at temperatures below $`k_BT/\epsilon =0.12`$ LJ<sub>13</sub> is dominated by structures close to the lowest energy icosahedral isomer, we expect the Jw and PT methods to have similar efficiencies (as measured by $`\mathrm{\Gamma }`$) for all temperatures.

## VI Conclusions

In this paper we have presented a study of the approach to the ergodic limit in MC simulations. In all the cases examined, the behavior of the MC metric $`d_k`$ can be approximated by Eq. (58), and the behavior of $`\mathrm{{\rm Y}}_{bK}`$ satisfies Eq. (60). Because the exponent $`\upsilon `$ is smaller than one for all the cases studied, the dependence of the non-diffusive contributions on $`d_k`$ is weaker (in the sense of Appendix A) than the diffusive contributions. The assumption on which we have built the stochastic model have been verified numerically for a system having a sufficiently complex potential surface to be viewed as prototypical of a large set of many-particle systems.

The MC metric used in this work appears to be a valuable tool to study the ergodicity properties of MC simulations. The non-ergodic components of the MC metric enable the prediction of the minimum length a MC simulation must have in order to be considered ergodic. The comparison of $`\mathrm{\Gamma }`$ from different algorithms gives a reasonable estimate of their relative efficiencies.

From the study of the melting region of 13 particle clusters, we have found that the exponent $`\upsilon `$ depends both on the method used and the nature of the potential energy function. We have performed calculations, not discussed in this work, where the functional form of the potential energy is modified. These studies have shown $`\upsilon `$ to be dependent on the details of the potential. We have not found the exponent $`\upsilon `$ to be a strong function of method. Although PT and Met have significantly different efficiencies as measured by their relative diffusion coefficients, $`\upsilon `$ is nearly the same in the two methods. The difference in the decay of $`\mathrm{{\rm Y}}_K`$ appears to be dominated by the coefficient in Eqs.(60) and (64) rather than the exponent.

As discussed in the text, parallel tempering and J-walking studies of many-particle systems must have an initial high temperature component that is chosen so that a Met simulation is known to be ergodic. For cluster simulations that require an external constraining potential to define the cluster, the radius of the constraining potential must be carefully chosen in order to achieve ergodic results. We have found the metric and associated decay laws developed in this work to be a particularly valuable method of choosing these initial parameters in both parallel tempering and J-walking simulations.

We also remark that the metric introduced here may be a more sensitive probe of ergodicity than may be required in some applications. For example in previous J-walking studies of the 13-particle Lennard-Jones cluster, the heat capacity curve determined with a constraining radius of 4$`\sigma `$ is nearly indistinguishable from the curve obtained with a constraining radius of 2$`\sigma `$. From the results of this work, we know the initial high temperature walk is not ergodic when a constraining radius of 4$`\sigma `$ is used. It is striking that the non-ergodicity as measured by the energy metric is not apparent in the heat capacity curve.

We have constructed a metric based on an ensemble of MC trajectories. By using an ensemble we attempt to cover sufficient portions of space so that all components are accessible. In practice only a finite subset of a full ensemble can be included, and it is always possible that components of space are missed. In such a case $`\mathrm{{\rm Y}}_K`$ may decay to zero numerically within the subspace, and the behavior may give misleading evidence that the simulation is ergodic. Because components of space may be missed in any finite simulation, it is impossible to guarantee ergodicity. It is hoped by using a sufficiently large ensemble of trajectories to define the metric, the possibility of missing components is minimized.

## Acknowledgments

We would like to thank Dr. O. Osenda for helpful comments. This work has been supported in part by the National Science Foundation under grant numbers CHE-9714970 and CDA-9724347. This research has been supported in part by the Phillips Laboratory, Air Force Material Command, USAF, through the use of the MHPCC under cooperative agreement number F29601-93-0001. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Phillips Laboratory or the U.S. Government.

## A Weak dependence of the non-diffusive contributions

We have considered two overall time scales for a MC simulation. Properties calculated at short times (labeled $`k`$ in the discrete case) provide information about each step of the MC process, and properties averaged over the total simulation time (labeled $`K`$ in the discrete case) give information about the approach to ergodic behavior. When $`K`$ is sufficiently short we have both diffusive and non-diffusive contributions as a function of $`k`$. In this Appendix we explain the relative time dependence of the diffusive and non-diffusive contributions to the autocorrelation function.

It has been assumed that the autocorrelation function Eq. (35) can be expressed as the sum of diffusive terms plus non-diffusive terms, i.e.

$$\kappa (t,t^{})=\kappa _d(t,t^{})+\underset{\mathrm{}=\lambda +1}{\overset{\mathrm{\Lambda }}{}}\kappa _{nd,\mathrm{}}(t,t^{}),$$
(A1)
where

$`\kappa _d(t,t^{})`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_0+\mathrm{\Gamma }_1+\mathrm{\Gamma }_2+\mathrm{}+\mathrm{\Gamma }_\lambda }{t_>}}`$ (A2)
$`\kappa _{nd,\mathrm{}}(t,t^{})`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_{\mathrm{}}}{\tau _{\mathrm{}}}}f_{\mathrm{}}\left({\displaystyle \frac{t_{\mathrm{}}^{}}{\tau _{\mathrm{}}}}\right).`$ (A3)
Increasing the time variables by a factor $`b>1`$, such that $`\tau _\lambda bt_>\tau _{\lambda +1}`$, with $`\lambda 1`$, we can study the relative variations of each contribution to the correlation function, diffusive and non-diffusive (labeled by $`\mathrm{}>\lambda `$). In this Appendix we only consider values of $`b`$ such that the transformation $`tbt`$ does not increase the time scale beyond the local correlation time. In Appendix B values of $`b`$ are considered that do cross such time scales.

By relative variations we mean

$`\mathrm{\Delta }_d(t,t^{};b)`$ $`=`$ $`\left|{\displaystyle \frac{\kappa _d(bt,bt^{})\kappa _d(t,t^{})}{\kappa _d(t,t^{})}}\right|`$ (A4)
$`\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)`$ $`=`$ $`\left|{\displaystyle \frac{\kappa _{nd,\mathrm{}}(bt,bt^{})\kappa _{nd,\mathrm{}}(t,t^{})}{\kappa _{nd,\mathrm{}}(t,t^{})}}\right|.`$ (A5)
The relative variation of each non-diffusive contribution is

$$\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)=\left|1\frac{1}{b^2}\frac{_0^{bt}𝑑t_1_0^{bt^{}}𝑑t_2f_{\mathrm{}}\left(\frac{|t_1t_2|}{\tau _{\mathrm{}}}\right)}{_0^t𝑑t_1_0^t^{}𝑑t_2f_{\mathrm{}}\left(\frac{|t_1t_2|}{\tau _{\mathrm{}}}\right)}\right|.$$
(A6)
whereas the relative variation of the diffusive contribution is

$$\mathrm{\Delta }_d(t,t^{};b)=1\frac{1}{b},$$
(A7)
If $`\mathrm{\Delta }_d(t,t^{};b)>\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)`$ for all pair of times $`t`$ and $`t^{}`$ and for all $`b>1`$ such that $`bt_>\tau _{\mathrm{}}`$, we say that the non-diffusive contributions are weaker than the diffusive contribution in their dependence on $`t`$. We explore, in the remainder of this appendix, the properties $`f_{\mathrm{}}`$ must have in order that the inequality $`\mathrm{\Delta }_d(t,t^{};b)>\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)`$ is satisfied.

Lemma: If the function $`H_{\mathrm{}}(t;\tau )`$

$$H_{\mathrm{}}(t;\tau )=_0^t𝑑t^{}f_{\mathrm{}}\left(\frac{t^{}}{\tau }\right)>0$$
(A8)
satisfies the inequality

$$H_{\mathrm{}}(t;\tau )>tf_{\mathrm{}}\left(\frac{t}{\tau }\right)t\mathrm{and}\tau ,$$
(A9)
then, $`H_{\mathrm{}}(t;\tau )`$ is an increasing function of $`\tau `$.

Demonstration: For $`\mathrm{}`$ and $`t`$ fixed, the function $`H_{\mathrm{}}(t,\tau )`$ evaluated in $`\tau ^{}`$ is

$`H_{\mathrm{}}(t;\tau ^{})`$ $`=`$ $`{\displaystyle _0^t}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau ^{}}}\right)`$ (A10)
$`=`$ $`{\displaystyle _0^t}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{\tau t^{}}{\tau ^{}\tau }}\right)`$ (A11)
$`=`$ $`{\displaystyle \frac{\tau ^{}}{\tau }}{\displaystyle _0^{\tau t/\tau ^{}}}𝑑uf_{\mathrm{}}\left({\displaystyle \frac{u}{\tau }}\right)`$ (A12)
$`=`$ $`{\displaystyle \frac{\tau ^{}}{\tau }}H_{\mathrm{}}(\tau t/\tau ^{};\tau ),`$ (A13)
then, for $`\mathrm{\Delta }\tau >0`$

$`{\displaystyle \frac{H_{\mathrm{}}(t;\tau +\mathrm{\Delta }\tau )H_{\mathrm{}}(t;\tau )}{\mathrm{\Delta }\tau }}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }\tau }}\left\{{\displaystyle \frac{\tau +\mathrm{\Delta }\tau }{\tau }}{\displaystyle _0^{\tau t/(\tau +\mathrm{\Delta }\tau )}}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau }}\right){\displaystyle _0^t}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau }}\right)\right\}`$ (A14)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }\tau }}\left\{{\displaystyle \frac{\mathrm{\Delta }\tau }{\tau }}{\displaystyle _0^{\tau t/(\tau +\mathrm{\Delta }\tau )}}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau }}\right){\displaystyle _{\tau t/(\tau +\mathrm{\Delta }\tau )}^t}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau }}\right)\right\}`$ (A15)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }\tau }}\left\{{\displaystyle \frac{\mathrm{\Delta }\tau }{\tau }}{\displaystyle _0^{\tau t/(\tau +\mathrm{\Delta }\tau )}}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau }}\right){\displaystyle \frac{t\mathrm{\Delta }\tau }{\tau +\mathrm{\Delta }\tau }}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau }}\right)\right\},`$ (A16)
where $`t^{}[t\tau /(\tau +\mathrm{\Delta }\tau ),t]`$. In the limit $`\mathrm{\Delta }\tau 0`$, and by virtue of the continuity of $`f_{\mathrm{}}`$, the derivative takes the form

$$\frac{H_{\mathrm{}}(t;\tau )}{\tau }=\frac{1}{\tau }\left\{H_{\mathrm{}}(t;\tau )tf_{\mathrm{}}\left(\frac{t}{\tau }\right)\right\}.$$
(A17)
Then, $`H_{\mathrm{}}(t;\tau )/\tau >0`$, and $`H_{\mathrm{}}(t;\tau )`$ is an increasing function of $`\tau .\mathrm{}`$

Here we have presented the two first conditions $`f_{\mathrm{}}`$ must have, namely Eqs. (A8) and (A9). From Eq. (23) $`f_{\mathrm{}}(0)`$ is a global maximum, and the memory functions must have a positive peak at zero. The area below that peak must be sufficiently large to satisfy Eq. (A8). Moreover, $`f_{\mathrm{}}(0)`$ must be sufficiently large to satisfy Eq. (A9), even at points where $`f_{\mathrm{}}(t/\tau )`$ is a local maximum. Then, to satisfy this Lemma, we need a memory function with a sufficiently large global maximum at $`t=0`$.

Corollary: Suppose $`H_{\mathrm{}}(t;\tau _{\mathrm{}})>tf_{\mathrm{}}(t/\tau _{\mathrm{}})`$. If $`b>1`$, then $`0<\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)<1`$ for all pair of times $`t`$ and $`t^{}`$.

Demonstration: Under the change of scale in time $`tbt`$, $`\kappa _{nd,\mathrm{}}(t,t^{})`$ can be written

$`\kappa _{nd,\mathrm{}}(bt,bt^{})`$ $`=`$ $`{\displaystyle \frac{1}{b^2tt^{}}}{\displaystyle _0^{bt}}𝑑t_1{\displaystyle _0^{bt^{}}}𝑑t_2{\displaystyle \frac{1}{\tau _{\mathrm{}}}}f_{\mathrm{}}\left({\displaystyle \frac{|t_1t_2|}{\tau _{\mathrm{}}}}\right)`$ (A18)
$`=`$ $`{\displaystyle \frac{1}{tt^{}}}{\displaystyle _0^t}𝑑t_1{\displaystyle _0^t^{}}𝑑t_2{\displaystyle \frac{1}{\tau _{\mathrm{}}}}f_{\mathrm{}}\left({\displaystyle \frac{b|t_1t_2|}{\tau _{\mathrm{}}}}\right),`$ (A19)
then, the quotient $`\kappa _{nd,\mathrm{}}(bt,bt^{})/\kappa _{nd,\mathrm{}}(t,t^{})`$ is

$`{\displaystyle \frac{\kappa _{nd,\mathrm{}}(bt,bt^{})}{\kappa _{nd,\mathrm{}}(t,t^{})}}`$ $`=`$ $`{\displaystyle \frac{_0^t𝑑t_1\left\{_0^{t_1}𝑑tf_{\mathrm{}}\left(\frac{t}{\tau _{\mathrm{}}/b}\right)+_0^{t_>t_1}𝑑tf_{\mathrm{}}\left(\frac{t}{\tau _{\mathrm{}}/b}\right)\right\}}{_0^t𝑑t_1\left\{_0^{t_1}𝑑tf_{\mathrm{}}\left(\frac{t}{\tau _{\mathrm{}}}\right)+_0^{t_>t_1}𝑑tf_{\mathrm{}}\left(\frac{t}{\tau _{\mathrm{}}}\right)\right\}}}`$ (A20)
$`=`$ $`{\displaystyle \frac{_0^t𝑑t_1\left\{H_{\mathrm{}}(t_1;\tau _{\mathrm{}}/b)+H_{\mathrm{}}(t_>t_1;\tau _{\mathrm{}}/b)\right\}}{_0^t𝑑t_1\left\{H_{\mathrm{}}(t_1;\tau _{\mathrm{}})+H_{\mathrm{}}(t_>t_1;\tau _{\mathrm{}})\right\}}}.`$ (A21)
By Eq. (A8), $`H_{\mathrm{}}(t;\tau )>0`$ $`t`$ and $`\tau `$. By the Lemma the numerator is smaller than the denominator. Then $`0<\kappa _{nd,\mathrm{}}(bt,bt^{})/\kappa _{nd,\mathrm{}}(t,t^{})<1`$ and then, $`0<\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)<1.\mathrm{}`$

Theorem: Suppose that $`b>1`$ is such that $`\tau _\mathrm{}1bt_>\tau _{\mathrm{}}`$, $`H_{\mathrm{}}(t;\tau _{\mathrm{}})>tf_{\mathrm{}}(t/\tau _{\mathrm{}})`$, and all $`f_{\mathrm{}}`$ satisfy the Lipschitz condition  (for all closed interval $`𝒜`$ exists a real positive number $`C_{\mathrm{}}`$ such that

$$\left|f_{\mathrm{}}(x)f_{\mathrm{}}(y)\right|C_{\mathrm{}}\left|xy\right|$$
(A22)
for all $`x`$ and $`y`$ in $`𝒜`$). Then $`\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)<\mathrm{\Delta }_d(t,t^{};b)`$ if and only if $`f_{\mathrm{}}`$ is non-negative in the interval $`[0,t_>)`$.

Demonstration: If $`\mathrm{\Delta }_{nd,\mathrm{}}(t,t^{};b)<\mathrm{\Delta }_d(t,t^{};b)`$, then

$`1{\displaystyle \frac{1}{b}}`$ $`>`$ $`1{\displaystyle \frac{1}{b^2}}{\displaystyle \frac{_0^{bt}𝑑t_1_0^{bt^{}}𝑑t_2f_{\mathrm{}}\left(\frac{|t_1t_2|}{\tau _{\mathrm{}}}\right)}{_0^t𝑑t_1_0^t^{}𝑑t_2f_{\mathrm{}}\left(\frac{|t_1t_2|}{\tau _{\mathrm{}}}\right)}}`$ (A23)
$`1`$ $`<`$ $`{\displaystyle \frac{1}{b}}{\displaystyle \frac{_0^{bt}𝑑t_1_0^{bt^{}}𝑑t_2f_{\mathrm{}}\left(\frac{|t_1t_2|}{\tau _{\mathrm{}}}\right)}{_0^t𝑑t_1_0^t^{}𝑑t_2f_{\mathrm{}}\left(\frac{|t_1t_2|}{\tau _{\mathrm{}}}\right)}}`$ (A24)
where the operations to reach Eq. (A24) are valid by using Corollary. Then

$`0`$ $`<`$ $`{\displaystyle _0^{bt}}𝑑t_1{\displaystyle _0^{bt^{}}}𝑑t_2{\displaystyle \frac{1}{b}}f_{\mathrm{}}\left({\displaystyle \frac{|t_1t_2|}{\tau _{\mathrm{}}}}\right){\displaystyle _0^t}𝑑t_1{\displaystyle _0^t^{}}𝑑t_2f_{\mathrm{}}\left({\displaystyle \frac{|t_1t_2|}{\tau _{\mathrm{}}}}\right)`$ (A25)
$`0`$ $`<`$ $`{\displaystyle _0^t}𝑑t_1{\displaystyle _0^t^{}}𝑑t_2\left\{bf_{\mathrm{}}\left({\displaystyle \frac{b|t_1t_2|}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{|t_1t_2|}{\tau _{\mathrm{}}}}\right)\right\}`$ (A26)
$`0`$ $`<`$ $`{\displaystyle _0^{t_<}}𝑑t_1\left\{{\displaystyle _0^{t_1}}𝑑t_2\left[bf_{\mathrm{}}\left({\displaystyle \frac{b(t_1t_2)}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{t_1t_2}{\tau _{\mathrm{}}}}\right)\right]+{\displaystyle _{t_1}^{t_>}}𝑑t_2\left[bf_{\mathrm{}}\left({\displaystyle \frac{b(t_2t_1)}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{t_2t_1}{\tau _{\mathrm{}}}}\right)\right]\right\}`$ (A27)
$`0`$ $`<`$ $`{\displaystyle _0^{t_<}}𝑑t_1\left\{{\displaystyle _0^{t_1}}𝑑t\left[bf_{\mathrm{}}\left({\displaystyle \frac{bt}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right]+{\displaystyle _0^{t_>t_1}}𝑑t\left[bf_{\mathrm{}}\left({\displaystyle \frac{bt}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right]\right\}`$ (A28)
$`0`$ $`<`$ $`{\displaystyle _0^{t_<}}𝑑t_1\left\{{\displaystyle _0^{bt_1}}𝑑tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right){\displaystyle _0^{t_1}}𝑑tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+{\displaystyle _0^{b(t_>t_1)}}𝑑tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right){\displaystyle _0^{t_>t_1}}𝑑tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right\}`$ (A29)
$`0`$ $`<`$ $`{\displaystyle _0^{t_<}}𝑑t_1\left\{{\displaystyle _{t_1}^{bt_1}}𝑑tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+{\displaystyle _{t_>t_1}^{b(t_>t_1)}}𝑑tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right\}.`$ (A30)
Using the intermediate value theorem,  we have

$`{\displaystyle _t^{bt}}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau _{\mathrm{}}}}\right)`$ $`=`$ $`(b1)tf_{\mathrm{}}\left({\displaystyle \frac{t^{}(t)}{\tau _{\mathrm{}}}}\right)`$ (A31)
$`=`$ $`(b1)tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+(b1)t\left[f_{\mathrm{}}\left({\displaystyle \frac{t^{}(t)}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right],`$ (A32)
where $`t^{}(t)[t,bt]`$. Let be $`t_\alpha ^{}(t)`$ and $`t_\beta ^{}(t)`$ the values at which the intermediate value theorem is satisfied, in the intervals $`[t,bt]`$ and $`[t_>t,b(t_>t)]`$ respectively

$`(b1)tf_{\mathrm{}}\left({\displaystyle \frac{t_\alpha ^{}(t)}{\tau _{\mathrm{}}}}\right)`$ $`=`$ $`{\displaystyle _t^{bt}}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau _{\mathrm{}}}}\right)`$ (A33)
$`(b1)(t_>t)f_{\mathrm{}}\left({\displaystyle \frac{t_\beta ^{}(t)}{\tau _{\mathrm{}}}}\right)`$ $`=`$ $`{\displaystyle _{t_>t}^{b(t_>t)}}𝑑t^{}f_{\mathrm{}}\left({\displaystyle \frac{t^{}}{\tau _{\mathrm{}}}}\right),`$ (A34)
then, the remainder can be written as

$$R_{\mathrm{}}(t_<,t_>;b)=_0^{t_<}𝑑t\left\{t\left[f_{\mathrm{}}\left(\frac{t_\alpha ^{}(t)}{\tau _{\mathrm{}}}\right)f_{\mathrm{}}\left(\frac{t}{\tau _{\mathrm{}}}\right)\right]+(t_>t)\left[f_{\mathrm{}}\left(\frac{t_\beta ^{}(t)}{\tau _{\mathrm{}}}\right)f_{\mathrm{}}\left(\frac{t_>t}{\tau _{\mathrm{}}}\right)\right]\right\}.$$
(A35)
By the Lipschitz condition, we have that

$`R_{\mathrm{}}(t_<,t_>;b)`$ $``$ $`{\displaystyle _0^{t_<}}𝑑t\left\{t\left|f_{\mathrm{}}\left({\displaystyle \frac{t_\alpha ^{}(t)}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right|+(t_>t)\left|f_{\mathrm{}}\left({\displaystyle \frac{t_\beta ^{}(t)}{\tau _{\mathrm{}}}}\right)f_{\mathrm{}}\left({\displaystyle \frac{t_>t}{\tau _{\mathrm{}}}}\right)\right|\right\}`$ (A36)
$`<`$ $`{\displaystyle _0^{t_<}}𝑑t\left\{tC_{\mathrm{}}\left|{\displaystyle \frac{t_\alpha ^{}(t)t}{\tau _{\mathrm{}}}}\right|+(t_>t)C_{\mathrm{}}\left|{\displaystyle \frac{t_\beta ^{}(t)(t_>t)}{\tau _{\mathrm{}}}}\right|\right\}`$ (A37)
$`<`$ $`{\displaystyle \frac{C_{\mathrm{}}}{\tau _{\mathrm{}}}}{\displaystyle _0^{t_<}}𝑑t\left\{t\left|btt\right|+(t_>t)\left|b(t_>t)(t_>t)\right|\right\}`$ (A38)
$`<`$ $`{\displaystyle \frac{C_{\mathrm{}}}{\tau _{\mathrm{}}}}(b1){\displaystyle _0^{t_<}}𝑑t\left[t^2+(t_>t)^2\right]`$ (A39)
$`<`$ $`{\displaystyle \frac{C_{\mathrm{}}}{\tau _{\mathrm{}}}}(b1)\left({\displaystyle \frac{2}{3}}t_<^3+t_<t_>(t_>t_<)\right)`$ (A40)
$`<`$ $`{\displaystyle \frac{2}{3}}t_>^3{\displaystyle \frac{C_{\mathrm{}}}{\tau _{\mathrm{}}}}(b1),`$ (A41)
where $`C_{\mathrm{}}`$ is a suitable positive real constant. Using Eqs. (A32) and (A35) in Eq. (A30) we have

$`0`$ $`<`$ $`{\displaystyle _0^{t_<}}𝑑t(b1)\left\{tf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+(t_>t)f_{\mathrm{}}\left({\displaystyle \frac{t_>t}{\tau _{\mathrm{}}}}\right)\right\}+(b1)R_{\mathrm{}}(t_<,t_>;b)`$ (A42)
$`0`$ $`<`$ $`{\displaystyle _0^{t_<}}𝑑ttf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+{\displaystyle _{t_>t_<}^{t_>}}𝑑ttf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+R_{\mathrm{}}(t_<,t_>;b)`$ (A43)
$`0`$ $`<`$ $`{\displaystyle _0^{t_<}}𝑑ttf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+{\displaystyle _0^{t_>}}𝑑ttf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right){\displaystyle _0^{t_>t_<}}𝑑ttf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+R_{\mathrm{}}(t_<,t_>;b)`$ (A44)
$`0`$ $`<`$ $`F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<)+{\displaystyle \frac{2}{3}}t_>^3{\displaystyle \frac{C_{\mathrm{}}}{\tau _{\mathrm{}}}}(b1),`$ (A45)
where

$$F_{\mathrm{}}(t)=_0^t𝑑t^{}t^{}f_{\mathrm{}}\left(\frac{t^{}}{\tau _{\mathrm{}}}\right),$$
(A46)
is a continuous and differentiable function of $`t`$. The inequality (A45) holds for any $`b>1`$. Suppose that $`F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<)<0`$. Then, if $`b`$ is such that

$$b=1+\frac{3}{2L}\frac{\tau _{\mathrm{}}}{t_>^3C_{\mathrm{}}}\left|F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<)\right|,$$
(A47)
where $`L>2`$, we have that

$`0`$ $`<`$ $`F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<)+{\displaystyle \frac{1}{L}}\left|F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<)\right|`$ (A48)
$`0`$ $`<`$ $`{\displaystyle \frac{L1}{L}}\left[F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<)\right]`$ (A49)
in contradiction with the hypothesis that $`F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<)`$ is negative. Then

$$0F_{\mathrm{}}(t_<)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>t_<).$$
(A50)
Let us define the function

$$\mathrm{\Delta }F_{\mathrm{}}(t)=F_{\mathrm{}}(t)F_{\mathrm{}}(t_>t),$$
(A51)
where $`t(0,t_>)`$. The right derivative at $`t=0`$ of $`\mathrm{\Delta }F_{\mathrm{}}(t)`$ is

$`\underset{\mathrm{\Delta }t0^+}{lim}{\displaystyle \frac{\mathrm{\Delta }F_{\mathrm{}}(\mathrm{\Delta }t)\mathrm{\Delta }F_{\mathrm{}}(0)}{\mathrm{\Delta }t}}`$ $`=`$ $`\underset{\mathrm{\Delta }t0^+}{lim}{\displaystyle \frac{F_{\mathrm{}}(\mathrm{\Delta }t)F_{\mathrm{}}(0)+F_{\mathrm{}}(t_>)F_{\mathrm{}}(t_>\mathrm{\Delta }t)}{\mathrm{\Delta }t}}`$ (A52)
$`=`$ $`\underset{\mathrm{\Delta }t0^+}{lim}{\displaystyle \frac{1}{\mathrm{\Delta }t}}\left\{{\displaystyle _0^{\mathrm{\Delta }t}}𝑑ttf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)+{\displaystyle _{t_>\mathrm{\Delta }t}^{t_>}}𝑑ttf_{\mathrm{}}\left({\displaystyle \frac{t}{\tau _{\mathrm{}}}}\right)\right\}`$ (A53)
$`=`$ $`\underset{\mathrm{\Delta }t0^+}{lim}{\displaystyle \frac{1}{\mathrm{\Delta }t}}\left\{\mathrm{\Delta }tt_1^{}f_{\mathrm{}}\left({\displaystyle \frac{t_1^{}}{\tau _{\mathrm{}}}}\right)+\mathrm{\Delta }tt_2^{}f_{\mathrm{}}\left({\displaystyle \frac{t_2^{}}{\tau _{\mathrm{}}}}\right)\right\}`$ (A54)
where $`t_1^{}[0,\mathrm{\Delta }t]`$ and $`t_2^{}[t_>\mathrm{\Delta }t,t_>]`$. Thus

$$\frac{\mathrm{\Delta }F_{\mathrm{}}(t)}{t}|_{t0^+}=t_>f_{\mathrm{}}\left(\frac{t_>}{\tau _{\mathrm{}}}\right).$$
(A55)
If the right derivative at 0 of $`\mathrm{\Delta }F_{\mathrm{}}(t)`$ is negative, $`\mathrm{\Delta }F_{\mathrm{}}(t)`$ approaches $`F_{\mathrm{}}(t_>)`$ from below, when $`t0`$. There exists a time $`0<\stackrel{~}{t}<t_>`$, such that $`0>F_{\mathrm{}}(t_>)+\mathrm{\Delta }F_{\mathrm{}}(\stackrel{~}{t})`$, in contradiction with Eq. (A50). Then, $`f_{\mathrm{}}`$ must be non-negative for $`t(0,t_>)`$. By the property Eq. (23) $`f_{\mathrm{}}(0)`$ must be positive. This proves that $`\mathrm{\Delta }_{nd,\mathrm{\hspace{0.17em}1}}(t,t^{};b)<\mathrm{\Delta }_d(t,t^{};b)f_{\mathrm{}}(t)0`$ for $`0t<t_>`$. To demonstrate that if $`f_{\mathrm{}}`$ is positive yields $`\mathrm{\Delta }_{nd,\mathrm{\hspace{0.17em}1}}(t,t^{};b)<\mathrm{\Delta }_d(t,t^{};b)`$ (i.e. the converse), follow the argument backwards, from Eq. (A30). $`\mathrm{}`$

In conclusion, if the memory functions are positive, satisfy the Lipschitz condition, and satisfy the condition Eqs. (A8) and (A9), the non-diffusive contributions are more weakly dependent on time than $`1/t`$.

The results of the present appendix are valid in the limit of a complete ensemble. In our numerical experiments only partial samples of the ensemble can be considered. The memory functions that appear in our numerical calculations come from partial mean values of the product of discontinuous functions (every noise process is a discontinuous function). These memory functions are discontinuous. The behavior of the non-diffusive contributions observed in our numerical experiments is in agreement with these analytic (infinite ensemble limit) results. We can infer that there might be a version of the theorem applied to discontinuous memory functions, but we have been unable to develop such a theorem.

## B Consequences of the time scale change in the non-diffusive contributions

In this appendix we show the behavior of the function $`f_1`$ when its correlation time is changed according to $`\tau _1\tau _{b1}=\tau _1/b`$, with $`b1`$; i.e. when the total simulation time is scaled to exceed the correlation time of the first colored noise process.

We multiply the time variables by a number $`b`$, such that $`\tau _1bt_>\tau _2`$. We have that the $`g_1`$ process contributes to the autocorrelation function with

$`{\displaystyle \frac{1}{b^2tt^{}}}G_1(bt/\tau _{\mathrm{}})G_1(bt^{}/\tau _{\mathrm{}})`$ $`=`$ $`{\displaystyle \frac{1}{b^2tt^{}}}{\displaystyle _0^{bt_<}}𝑑t_1{\displaystyle _0^{bt_>}}𝑑t_2{\displaystyle \frac{1}{\tau _1}}f_1\left({\displaystyle \frac{|t_1t_2|}{\tau _1}}\right)`$ (B1)
$`=`$ $`{\displaystyle \frac{1}{btt^{}}}{\displaystyle _0^{t_<}}𝑑t_1^{}{\displaystyle _0^{t_>}}𝑑t_2^{}{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t_1^{}t_2^{}|}{\tau _{b1}}}\right)`$ (B2)
where $`t^{}=t/b`$ and $`\tau _{b1}=\tau _1/b`$. We want to compute this contribution both within the neighborhood $`t_1=t_2`$ as well as outside such a region. To do so, we can split the integral in Eq. (B2) in three parts

$$\frac{1}{b^2tt^{}}G_1(bt/\tau _{\mathrm{}})G_1(bt^{}/\tau _{\mathrm{}})=I_1+I_2+I_3$$
(B3)
where

$`I_1`$ $`=`$ $`{\displaystyle \frac{1}{btt^{}}}{\displaystyle _0^{t_<}}𝑑t_1{\displaystyle _0^{\mathrm{max}(0,t_1ϵ/2)}}𝑑t_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{t_1t_2}{\tau _{b1}}}\right)`$ (B4)
$`I_2`$ $`=`$ $`{\displaystyle \frac{1}{btt^{}}}{\displaystyle _0^{t_<}}𝑑t_1{\displaystyle _{\mathrm{max}(0,t_1ϵ/2)}^{\mathrm{min}(t_>,t_1+ϵ/2)}}𝑑t_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t_1t_2|}{\tau _{b1}}}\right)`$ (B5)
$`I_3`$ $`=`$ $`{\displaystyle \frac{1}{btt^{}}}{\displaystyle _0^{t_<}}𝑑t_1{\displaystyle _{\mathrm{min}(t_>,t_1+ϵ/2)}^{t_>}}𝑑t_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{t_2t_1}{\tau _{b1}}}\right)`$ (B6)
with $`t_<>ϵ>0`$ (observe that the only integral involving $`t_1=t_2`$ is $`I_2`$). Consider $`I_1`$. If $`t_1<ϵ/2`$ the inner integral is zero. Therefore, $`t_1`$ must be bigger than $`ϵ/2`$ and

$$I_1=\frac{1}{bt_<t_>}_{ϵ/2}^{t_<}𝑑t_1_0^{t_1ϵ/2}𝑑t_2\frac{1}{\tau _{b1}}f_1\left(\frac{t_1t_2}{\tau _{b1}}\right),$$
(B7)
which, by virtue of the continuity of $`f_1`$, can be bounded as follows

$`{\displaystyle \frac{1}{bt_<t_>}}{\displaystyle _{ϵ/2}^{t_<}}𝑑t_1{\displaystyle \frac{b}{\tau _1}}\left(t_1{\displaystyle \frac{ϵ}{2}}\right)f_1\left({\displaystyle \frac{bt_{min}}{\tau _1}}\right)`$ $`I_1`$ $`{\displaystyle \frac{1}{bt_<t_>}}{\displaystyle _{ϵ/2}^{t_<}}𝑑t_1{\displaystyle \frac{b}{\tau _1}}\left(t_1{\displaystyle \frac{ϵ}{2}}\right)f_1\left({\displaystyle \frac{bt_{max}}{\tau _1}}\right)`$ (B8)
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{(t_<ϵ/2)^2}{t_<t_>}}{\displaystyle \frac{1}{\tau _1}}f_1\left({\displaystyle \frac{bt_{min}}{\tau _1}}\right)`$ $`I_1`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{(t_<ϵ/2)^2}{t_<t_>}}{\displaystyle \frac{1}{\tau _1}}f_1\left({\displaystyle \frac{bt_{max}}{\tau _1}}\right)`$ (B9)
where $`t_{max}`$ ($`t_{min}`$) is the time in the interval $`[ϵ/2,t_<]`$ at which the function $`f_1`$ reaches its maximum (minimum) value. Because $`f_1`$ is continuous, there exists $`t_1^{}[t_{min},t_{max}]`$ at which

$$I_1=\frac{1}{2}\frac{(t_<ϵ/2)^2}{t_<t_>}\frac{1}{\tau _1}f_1\left(\frac{bt_1^{}}{\tau _1}\right).$$
(B10)
Consider now $`I_3`$. If $`t_1+ϵ/2>t_>`$, the inner integral is zero. Therefore, $`0<t_1<\mathrm{min}(t_<,t_>ϵ/2)`$ and

$`I_3`$ $`=`$ $`{\displaystyle \frac{1}{btt^{}}}{\displaystyle _0^{\mathrm{min}(t_<,t_>ϵ/2)}}𝑑t_1{\displaystyle _{t_1+ϵ/2}^{t_>}}𝑑t_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{t_2t_1}{\tau _{b1}}}\right)`$ (B11)
$`=`$ $`{\displaystyle \frac{\mathrm{min}(t_<,t_>ϵ/2)}{t_<t_>}}\left[t_>{\displaystyle \frac{ϵ}{2}}{\displaystyle \frac{1}{2}}\mathrm{min}(t_<,t_>ϵ/2)\right]{\displaystyle \frac{1}{\tau _1}}f_1\left({\displaystyle \frac{bt_3^{}}{\tau _1}}\right),`$ (B12)
where $`t_3^{}[t_{min},t_{max}]`$, and now $`t_{max}`$ ($`t_{min}`$) is the time in $`[ϵ/2,t_>]`$ at which the function $`f_1`$ reaches its maximum (minimum) value.

Fig. 7: The area under the curve represents the first integral in Eq. (B13). The darker piece is half of the integral in the interval $`[t_1,t_1]`$, the lighter is half of the integral in $`[ϵ/2,ϵ/2]`$.

Let us consider now $`I_2`$. First observe that for the integral in $`t_1`$, if $`0t_1ϵ/2`$, $`\mathrm{max}(0,t_1ϵ/2)=0`$ and $`\mathrm{min}(t_>,t_1+ϵ/2)=t_1+ϵ/2`$. If $`ϵ/2t_1t_<`$ then $`\mathrm{max}(0,t_1ϵ/2)=t_1ϵ/2`$. Then

$$I_2=\frac{1}{bt_<t_>}\left\{_0^{ϵ/2}𝑑t_1_0^{t_1+ϵ/2}𝑑t_2\frac{1}{\tau _{b1}}f_1\left(\frac{|t_1t_2|}{\tau _{b1}}\right)+_{ϵ/2}^{t_<}𝑑t_1_{t_1ϵ/2}^{\mathrm{min}(t_>,t_1+ϵ/2)}𝑑t_2\frac{1}{\tau _{b1}}f_1\left(\frac{|t_1t_2|}{\tau _{b1}}\right)\right\}.$$
(B13)
The integral in $`t_2`$ between 0 and $`t_1+ϵ/2`$ can be evaluated with the help of Fig. 7

$$_0^{t_1+ϵ/2}𝑑t_2\frac{1}{\tau _{b1}}f_1\left(\frac{|t_1t_2|}{\tau _{b1}}\right)=\frac{1}{2}_{ϵ/2}^{ϵ/2}𝑑t\frac{1}{\tau _{b1}}f_1\left(\frac{|t|}{\tau _{b1}}\right)+\frac{1}{2}_{t_1}^{t_1}𝑑t\frac{1}{\tau _{b1}}f_1\left(\frac{|t|}{\tau _{b1}}\right).$$
(B14)
The second integral in $`t_1`$ can be separated in two parts; the first for $`ϵ/2t_1\mathrm{min}(t_<,t_>ϵ/2)`$ and the second for $`\mathrm{min}(t_<,t_>ϵ/2)t_1t_<`$. If $`t_>t_<<ϵ/2`$ the second term is zero. Then

$`{\displaystyle _{ϵ/2}^{t_<}}𝑑t_1{\displaystyle _{t_1ϵ/2}^{\mathrm{min}(t_>,t_1+ϵ/2)}}`$ $`dt_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t_1t_2|}{\tau _{b1}}}\right)={\displaystyle _{ϵ/2}^{\mathrm{min}(t_<,t_>ϵ/2)}}𝑑t_1{\displaystyle _{t_1ϵ/2}^{\mathrm{min}(t_>,t_1+ϵ/2)}}𝑑t_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t_1t_2|}{\tau _{b1}}}\right)+`$ (B16)
$`\mathrm{\Theta }\left({\displaystyle \frac{ϵ}{2}}+t_<t_>\right){\displaystyle _{t_>ϵ/2}^{t_<}}𝑑t_1{\displaystyle _{t_1ϵ/2}^{\mathrm{min}(t_>,t_1+ϵ/2)}}𝑑t_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t_1t_2|}{\tau _{b1}}}\right),`$
where $`\mathrm{\Theta }`$ is the step function. If $`t_1\mathrm{min}(t_<,t_>ϵ/2)`$ then $`\mathrm{min}(t_>,t_1+ϵ/2)=t_1+ϵ/2`$. The last integral in $`t_2`$ can be rearranged in the same way as Eq. (B14). Then

$`{\displaystyle _{ϵ/2}^{t_<}}𝑑t_1`$ $`{\displaystyle _{t_1ϵ/2}^{\mathrm{min}(t_>,t_1+ϵ/2)}}𝑑t_2{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t_1t_2|}{\tau _{b1}}}\right)={\displaystyle _{ϵ/2}^{\mathrm{min}(t_<,t_>ϵ/2)}}𝑑t_1{\displaystyle _{ϵ/2}^{ϵ/2}}𝑑t{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t|}{\tau _{b1}}}\right)+`$ (B18)
$`{\displaystyle \frac{1}{2}}\mathrm{\Theta }\left({\displaystyle \frac{ϵ}{2}}+t_<t_>\right){\displaystyle _{t_>ϵ/2}^{t_<}}𝑑t_1\left[{\displaystyle _{ϵ/2}^{ϵ/2}}𝑑t{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t|}{\tau _{b1}}}\right)+{\displaystyle _{t_1t_>}^{t_>t_1}}𝑑t{\displaystyle \frac{1}{\tau _{b1}}}f_1\left({\displaystyle \frac{|t|}{\tau _{b1}}}\right)\right],`$
We can observe that the correlation time $`\tau _{b1}`$ goes to zero when $`b`$ is increased. The function $`(1/\tau _1)f_1(bt/\tau _1)`$ becomes negligible outside a neighborhood of $`t=0`$ \[observe Eqs. (B10) and (B12)\]. Equation (18) holds, then, if $`b`$ is sufficiently large, $`(1/\tau _{b1})f_1(t/\tau _{b1})`$ can be considered a delta function. The integrals $`I_1`$ and $`I_3`$ become zero, and the integrals involving $`t=0`$ in the expression of $`I_2`$ converge to one. $`I_2`$ becomes

$$I_2=\frac{1}{bt_<t_>}\left\{\mathrm{min}(t_<,t_>\frac{ϵ}{2})+\mathrm{\Theta }\left(\frac{ϵ}{2}+t_<t_>\right)\left(\frac{ϵ}{2}+t_<t_>\right)\right\}=\frac{1}{bt_>},$$
(B19)
which is a diffusive contribution to the autocorrelation function. The autocorrelation function becomes then

$$\kappa (bt,bt^{})=\frac{\mathrm{\Gamma }_0+\mathrm{\Gamma }_1}{bt_>}+\underset{\mathrm{}=2}{\overset{\mathrm{\Lambda }}{}}\frac{\mathrm{\Gamma }_{\mathrm{}}}{\tau _{\mathrm{}}}f_{\mathrm{}}\left(\frac{t_b\mathrm{}^{}(t_>)}{\tau _b\mathrm{}}\right).$$
(B20)
The same argument can be used when $`b`$ is such that $`\tau _2bt_>\tau _3`$. After such changes in the time scale, the diffusion coefficient $`\mathrm{\Gamma }=\mathrm{\Gamma }_0+\mathrm{\Gamma }_1`$ is enlarged, and the non-diffusive contributions are reduced. There is an ultimate scale change, such that $`\tau _\mathrm{\Lambda }bt_>`$. Beyond this maximum time scale the process can be considered diffusive.


text:
# REFERENCES

Comment on “Confinement of Slave Particles in U(1) Gauge Theories of Strongly Interacting Electrons”

In a recent Letter , Nayak argued that slave partices are always confined in gauge theories of strongly-correlated electron systems like the t-J model. The argument mostly relies on Elitzur’s theorem and the fact that the effective gauge theory under discussion is at infinite coupling. On the other hand, in the previous papers  we studied dynamics of the gauge theory of the t-J model showing that the slave particles are in the Coulomb phase below certain critical temperature $`T_{\mathrm{CSS}}`$ that depends on hole doping. Thus, with a finite 3D coupling, the spin-charge separation occur, while in pure 2D, slave particles interact via logarithmic potential (in contrast to the linear-rising confining potential). We comment on the discrepancy between the results in Refs.1,2.

First, the existence of the deconfining Coulomb phase does not contradicts Elitzur’s theorem which says that expectation values of any gauge-variant quantities are vanishing. A well-known example is the weak-coupling phase of the U(1) gauge theory in $`(3+1)`$ dimensions which corresponds to the quantum electrodynamics. Mean field theory can also describe deconfining phases enforcing Elitzur’s theorem . Secondly, though the gauge coupling is infinite in the gauge theories of the t-J model as Nayak argued, this fact does not necessarily mean that the slave particles are confined. Actually, it is known at present that the pure SU(3) gauge theory in $`(3+1)`$ dimensions is always in the confinement phase at $`T=0`$ regardless of strength of the coupling. However, if light $`N_f`$-flavour fermions ($`N_f`$ is the number of flavours) are coupled to the SU(3) gauge theory, the system is in the deconfinemnt phase even at infinite coupling for $`N_f>7`$ . This means that matter fields which couple to the gauge field strongly influence the phase structure of the system. Similar phenomenon happens in the gauge theory of the t-J model as we showed in Ref.2. Spinons and holons generate there nontrivial gauge dynamics and the deconfinement phase appears at low $`T`$. We employed a non-perturbative method, not relying upon the perturbative calculations contrary to Nayak’s claim in Ref.1.

Let us go into details of the arguments in Ref.1. The effective low-energy Lagrangian (8) is obtained by simply ignoring high-energy modes in Eq.(6). However, in order to obtain an effective model, the high-energy (and high-momentum) modes must be integrated out. In some cases, this procedure does not generate any important terms that influence the dynamics of the low-energy modes. However, in the present case, it is conceivable that “kinetic terms” of the gauge field appear effectively from the high-energy modes. For 1D quantum spin models, this possibility is realized . The form of the continuum Lagrangian (8) is essential in order to solve the “constraints” on holons and spinons. In presence of a kinetic term of the gauge field, the constaints are not strictly satisfied by the low-energy modes, though they are strictly satisfied before integration over the high-energy modes. To judge whether the deconfinement of the slave particles occurs in low-energy quasiexcitations, careful investigation on the gauge dynamics of the low-energy theory that is faithful to the t-J model is necessary .

In Ref.1, the 1D t-J model at half filling is dicussed for illustrative purpose. However 1D gauge systems are quite different from those in higher dimensions because of the lack of the transverse direction. Study of the t-J model on a chain and a ladder in the framework of the slave-fermion gauge theory is made in Ref.6. Gauge field dynamics is nontrivial and it explains clearly and self consistently why quasiexciations on a chain and a ladder are so different, reflecting that the coefficient of the CP1 $`\theta `$-term is $`\pi `$ for chain whereas $`2\pi `$ for ladder.

For the 2D case, there are also problems to be solved before the conclusions in Ref.1 are accepted. Here again the possiblity of appearance of the gauge-field kinetic term is ignored. This is the essential point to discuss the confinement problem of the slave particles in the present context. Furthermore, restoration of the spin SU(2) symmetry must be shown in a convincing way from (25). This is a crucial test of the effective action (25).

To summarize, (1) Matter fields which couple to gauge field can change drastically the phase structure of the system, so study of the gauge dynamics is indispensable to identify quasiexcitations at low enegies. (2) To obtain the low-energy effective theory, possiblity of the gauge-field kinetic term should be addressed instead of taking the naive continuum limit. Appearance of the kinetic term prevents one from solving the constraint. (3) Rewriting the Lagrangian in terms of the gauge-invariant variables does not directly implies that particles with nonvanishing gauge charges are confined.

Ikuo Ichinose

Institute of Physics

University of Tokyo

Komaba, Tokyo, 153-8902, Japan

Tetsuo Matsui

Department of Physics

Kinki University

Higashi-Osaka, 577-8502, Japan