Acceptance rate is a thermodynamic function in local Monte Carlo algorithms

TL;DR

This paper explores the acceptance rate in local Monte Carlo algorithms, deriving analytic expressions for the 1D Ising model and numerically analyzing other models, revealing a near-linear relationship with energy.

Contribution

It provides the first analytic expression for the acceptance rate in the 1D Ising model and extensive numerical analysis across various models and dimensions.

Findings

Acceptance rate is a linear function of energy in the 1D Ising model.

Acceptance rate remains nearly linear with energy in critical regions of other models.

Variance of acceptance rate stays finite near phase transitions despite singularities in specific heat.

Abstract

We study properties of Markov chain Monte Carlo simulations of classical spin models with local updates. We derive analytic expressions for the mean value of the acceptance rate of single-spin-flip algorithms for the one-dimensional Ising model. We find that for the Metropolis algorithm the average acceptance rate is a linear function of energy. We further provide numerical results for the energy dependence of the average acceptance rate for the 3- and 4-state Potts model, and the XY model in one and two spatial dimensions. In all cases, the acceptance rate is an almost linear function of the energy in the critical region. The variance of the acceptance rate is studied as a function of the specific heat. While the specific heat develops a singularity in the vicinity of a phase transition, the variance of the acceptance rate stays finite.

Equations

$$H=-J\sum _{⟨ij⟩}{S}_i{S}_j,$$

$$H=-J\sum _{⟨ij⟩}\delta (S_i,S_j),$$

$$H=-J\sum _{⟨ij⟩}\cos (S_i-S_j),$$

$$\cdots \to \mu \to \nu \to \cdots ,$$

$$p(\mu \to \nu )=\min (1,e^{-\beta \Delta E}),$$

$$p(\mu \to \nu )=\frac{e^{-\beta E_{\nu }}}{e^{-\beta E_{\nu }}+e^{-\beta E_{\mu }}}.$$

$$p(\mu \to \nu )$$

$$H=-J\sum _{i=1}^L{S}_i{S}_{i+1},$$

$$Q_i=\frac{1}{2}\left(S_i{S}_{i+1}+1\right),$$

$$H=-2J\sum _{i=1}^L{Q}_i+JL,$$

$$\sum _{i=1}^L{Q}_i\text{is even.}$$

$$Z=2x^{-L/2}\sum _{l=0}^{L/2}{C}_L^{2l}{x}^{2l},$$

$$Z=x^{-L/2}\Bigl{[}(x+1)^{L}+(x-1)^{L}\Bigr{]}\;,$$

$$1-R=\frac{x^{2}-1}{Z}\Bigl{[}(x+1)^{L-2}+(x-1)^{L-2}\Bigr{]}x^{-L/2}\;.$$

$$R=\frac{2}{x+1}.$$

$$\epsilon =-\frac{x-1}{x+1},$$

$$R=1+\epsilon ,$$

$$R=\frac{x}{1+x^2}.$$

$$R=\frac{1}{2}\frac{1-{\epsilon }^2}{1+{\epsilon }^2},$$

$$C=J\frac{d\epsilon }{dT}=\frac{⟨H^2⟩-⟨H{⟩}^2}{V{T}^2},$$

$$1-R=\sum _{l=0}^{L/2}\left(1-x^{-2}\right)\frac{2l}{L}\frac{2l-1}{L-1}\frac{C_L^{2l}{x}^{2l}2x^{-L/2}}{Z}.$$

$$(x+1{)}^L=\sum _{k=0}^L{C}_L^k{x}^k.$$

$$1-R=\frac{x^{2}-1}{Z}\Bigl{[}(x+1)^{L-2}+(x-1)^{L-2}\Bigr{]}x^{-L/2}\;,$$

$$-E/J$$

$$-\epsilon =1-R$$

$$R=1+\epsilon$$

$$R=\sum_{l=0}^{L/2}\Bigl{(}\frac{1}{1+x^{2}}\frac{2l}{L}\frac{2l-1}{L-1}\\ +\frac{x^{2}}{1+x^{2}}\frac{L-2l}{L}\frac{L-2l-1}{L-1}\\ +\frac{L-2l}{L}\frac{2l}{L-1}\Bigr{)}\frac{C_{L}^{2l}x^{2l}2x^{-L/2}}{Z}\;,$$

$$R=\frac{x}{1+x^2}\frac{1-{\kappa }^L}{1+{\kappa }^L},$$

Full text

I Introduction

One of the most challenging problems in high-performance computing (HPC) is using the full strength of supercomputers. It is very demanding, however, to efficiently use all the power of a supercomputer in a single run. The main barrier is that most of the currently available algorithms do not scale well on the complex multi-node, multi-core, and multi-accelerator hybrid architectures which are the dominant today and will be dominant for the nearest future. Hence the computation must be divided into millions of tasks to be scheduled on individual cores. It is thus of crucial importance to develop new, fully scalable algorithms, new programming techniques, and new methods to build programs which can efficiently use the power of supercomputers.

Markov chain Monte Carlo methods  Berg ; BinderLandau are among the prime approaches of supercomputer simulations in physics, chemistry, and materials sciences. A Monte Carlo (MC) simulation is a sequence of local updates of microscopic degrees of freedom, and the overall performance of a simulation strongly depends on these local elementary updates BinderLandau . Two features of the family of local MC algorithms are worth noting: first, simulations based on these algorithms are naturally embarrassingly parallel, which makes them even more attractive for massively parallel computations; second, they are applicable in the presence of external fields and/or competing ferro- and antiferromagnetic interactions, where more sophisticated schemes (e.g., cluster updates) break down.

A program of the US DOE, Office of Science, entitled the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) INCITE awarded on average 60 projects per year, from which 30 percent are in physics, 28 percent in engineering, 15 percent in materials sciences, 9 percent in earth sciences, and 7 per cent in chemistry, etc. The average number of projects which used MC methods is 5 per year in the last ten years.

In the core of MC methods is the Metropolis algorithm Metropolis known for more than 60 years. What we find surprising is that a careful analysis of an “old” method could bring completely new knowledge on the computations. Namely, we note that the acceptance rate of an MC simulation is a thermodynamic function which displays a well-defined temperature dependence. These findings provide a deeper understanding of the algorithm, and we believe that our observation may influence further improvements of MC algorithms in general, and their parallel scalable versions.

A Markov process is defined by its set of transition probabilities between microscopic states, which must obey a balance condition Berg ; BinderLandau . Efficiency of an MC process is governed by several factors: the typical acceptance rate (i.e., the fraction of states in the Markov chain which differ from the previous state), the computational complexity of an elementary update, and the autocorrelation time of the process. Different choices for both trial moves and their transition probabilities are possible BinderLandau ; Janke2008 . A rigorous theorem claims that a specific choice of elementary updates — which is in fact a global update, attempting to update all degrees of freedom at once by drawing random increments from a Gaussian distribution — maximizes the efficiency if the mean value of the acceptance rate is tuned to the special value 0.234 Baesian ; Baesian2 . However, such processes have unfavorable correlation properties PotterSwendsen2013 , and are not suitable for MC simulations of physical systems with a large number of degrees of freedom.

We now make the following observation: the mean value of the acceptance rate of an MC simulation using a given set of transition probabilities has a well-defined temperature dependence. Therefore, it can be viewed as a thermodynamic function of the model under the MC dynamics, on par with other thermodynamic functions, e.g., the energy. A natural question is then: what is the relation between the mean acceptance rate and other thermodynamic quantities? We find that for the one-dimensional (1D) Ising model, the acceptance rate of the Metropolis algorithm Metropolis is a linear function of energy. An immediate question is then whether this linear relation is a one-off artifact of the Metropolis MC dynamics for the 1D Ising model, or whether it generalizes for other related models and MC algorithms.

The rest of the paper is organized as follows. In Sec. II we briefly describe the models for which we calculate analytically or compute numerically the acceptance rate. In Sec. III we present analytical results on the acceptance rate of the Metropolis and heat-bath algorithms for the 1D Ising model. Section IV contains computational results for the acceptance rate and its variance for the one-dimensional models described in Sec. II. Next Sec. V presents computational results for the acceptance rate and its variance for the two-dimensional models. In Sec. VI we summarize our findings. Two Appendixes contain more details on the analytical calculation of the acceptance rate when applying the Metropolis and heat-bath algorithms to the 1D Ising model.

II Models and Update Algorithms

We consider several well-known classical lattice models. The Ising model is defined by the Hamiltonian function

[TABLE]

where the coupling constant $J>0$, $\langle ij\rangle$ denotes nearest-neighbor pairs, and $S_{i}=\pm 1$ are Ising spins, located at the sites of a $d$-dimensional lattice of linear size $L$ (and volume $V=L^{d}$) with periodic boundary conditions.

In the $q$-state Potts model, spins can take $q$ possible values, $S_{i}\in{1,\ldots,q}$ Wu :

[TABLE]

where the coupling constant $J>0$ and $\delta(S_{i},S_{j})$ is the Kronecker delta symbol, which equals one whenever $S_{i}=S_{j}$, and zero otherwise.

Finally, we consider the XY model, defined by Kosterlitz :

[TABLE]

where the coupling constant $J>0$ and $S_{i}$ are continuous variables, $S_{i}\in[0,2\pi)$.

MC simulations provide a way of studying models (1)-(3) in thermodynamic equilibrium. An MC simulation constructs an ergodic random walk in the configuration space of a model,

[TABLE]

which generates the equilibrium Gibbs distribution of a model as its stationary distribution BinderLandau . For local updating schemes, successive configurations $\mu$ and $\nu$ only differ by the value of a single spin.

An elementary update of the local Metropolis algorithm Metropolis for the Ising model proceeds in two steps: (i) select a random site and (ii) flip its spin, $S_{i}\to-S_{i}$, with the probability

[TABLE]

where $\beta$ is the inverse temperature and $\Delta E=E_{\nu}-E_{\mu}$ the energy difference between the updated and original states Janke2008 . The generalization to models with more than two states per spin is straightforward: (ii) is simply replaced by $S_{i}\to\widetilde{S}_{i}$, where $\widetilde{S}_{i}$ is any admissible spin value.

The heat-bath algorithm for the Ising model differs from the Metropolis algorithm only in that a spin-flip update is accepted with the probability Janke2008

[TABLE]

This can be recast into the form

[TABLE]

which is the general Glauber update rule Glauber that can also be applied to models with more than two states per spin and a generalized update proposal $S_{i}\to\widetilde{S}_{i}$. Note that only for the Ising model with two states per spin, the heat-bath process coincides with the Glauber dynamics Glauber . In the general case, the heat-bath process in a strict sense (there is some confusion with the notation in the literature) involves all possible spin values and is hence more complicated.

III Acceptance rates of MC simulations of the 1D Ising model

To calculate the expected value of the acceptance probability of an MC simulation of the 1D Ising model we first convert Eq. (1) to bond variables Suzuki1972 ; Mueller2017 .

III.1 Bond representation

In one dimension Eq. (1) takes the form

[TABLE]

where the indices in (7) are taken modulo $L$, i.e., the term $S_{L}S_{L+1}$ is understood as $S_{L}S_{1}$.

To calculate the expected values of the acceptance probability of a MC simulation of the 1D Ising model (7) for the Metropolis and heat-bath updates, we first convert the model (7) to a bond representation. We define for a bond connecting sites $i$ and $i+1$ the ’charge’ Suzuki1972 ; Mueller2017 ,

[TABLE]

which takes values of 0 (for $S_{i}\neq S_{i+1}$) and 1 (for $S_{i}=S_{i+1}$ ). In this representation, Eq. (7) takes the form

[TABLE]

where the sum is taken over the bonds of the lattice. With periodic boundary conditions, the number of bonds equals the number of sites of the lattice. This way, the state space of the model (7) is spanned by a collection of $L$ integers $Q_{i}=\left\{0,1\right\}$, subject to the constraint: the parity of $\sum_{i=1}^{L}Q_{i}$ is the same as the parity of the number of bonds. We take $L$ to be even throughout, so that

[TABLE]

The partition function corresponding to Eq. (9) then reads

[TABLE]

where $x\equiv e^{2\beta J}$ and $\beta$ is the inverse temperature. The summation runs over the values of $\sum_{i}Q_{i}=2l$, and the binomial coefficient, $C_{L}^{2l}$, counts the number of ways of distributing $2l$ values of $Q=1$ over $L$ bonds. The factor of 2 accounts for a double-counting of the representation (8): each value of $Q_{i}$ can be realized by two possible combinations of $S_{i}$ and $S_{i+1}$ (e.g., $Q_{i}=0$ means that either $S_{i}=-1$ and $S_{i+1}=1$ or vice versa).

Performing the summation in (10) we obtain

[TABLE]

which agrees with the standard result Baxter .

III.2 Acceptance rate of the Metropolis algorithm

We use the bond representation (9) to calculate the acceptance rates for the Metropolis and the heat-bath algorithms. In this section, we only state the main results and relegate the details of calculations to the Appendix.

We start with the Metropolis update (4). Denoting the expected value of the acceptance probability by $R$, the expected value of the rejection probability is

[TABLE]

In the thermodynamic limit, $L\rightarrow\infty$, the second term in brackets is negligible, and Eq. (12) simplifies to

[TABLE]

We now compare Eq. (13) to the thermodynamic mean value of the internal energy of the system, $E$. Using the partition function (11), we obtain in the thermodynamic limit the standard result Baxter

[TABLE]

where the reduced energy density $\varepsilon=E/JL$.

Comparing Eqs. (13) and (14), we find

[TABLE]

i.e., the expected value of the acceptance probability is a linear function of the energy. In fact, relation (15) holds for all values of $L$, see Appendix A.

III.3 Acceptance rate of the heat-bath algorithm

For the expected value of the acceptance probability $R$ of the heat-bath update (5) we find in the thermodynamic limit $L\rightarrow\infty$

[TABLE]

Comparing to (14), we have

[TABLE]

which approaches the linear behavior $(1+\varepsilon)/2$ for $1+\varepsilon\ll 1$ (cf. Fig. 1). Details can be found in Appendix B.

IV Simulation results in 1D

In this section we first verify the analytical results for the Metropolis and heat-bath acceptance rates of the Ising model in one dimension and then test the observed qualitative features for the other models defined in Sec. II. Results for the generalization to two dimensions will be presented in the next section.

IV.1 First moments of the energy and acceptance rate

We performed MC simulations of the 1D Ising model (7) using Metropolis updates (4) and heat-bath updates (5) for temperatures ranging from $T/J=0.2$ to $10$. We used $N_{T}=10^{6}$ MC steps (MCS) for thermalization and collected statistics over $N_{A}=10^{7}$ MCS. Here an MCS is defined as $L$ elementary update attempts for a chain with $L$ spins. We first focus on the collected statistics for the total energy of the system, $E$, and the acceptance rate, $R$, which we specifically define as the ratio of the numbers of accepted and attempted elementary updates.

Figure 1 shows the relation between the mean values of the acceptance rates of the MC process and the reduced energy density, $\varepsilon$, for a chain of $L=512$ spins and periodic boundary conditions. The results of the MC simulations agree with Eqs. (15) and (17) in the whole range of energies (hence, temperatures). We note that for the heat-bath algorithm the dependence of $R$ on the reduced energy density is approximately linear in a wide range of temperatures: the relative difference between $R$ and its linear approximation is below 10% for $T/J<1.1$.

It is instructive to compare the behavior of local MC algorithms for related classical spin models. We performed MC simulations of the 3- and 4-state Potts models (2), and the classical XY model (3) in one dimension using the Metropolis and Glauber algorithms. There are several ways of organizing the local updates. We take the simplest possible prescription: we select a spin $S_{i}$ at random, and then draw a proposed value $\widetilde{S}_{i}$ for the update $S_{i}\to\widetilde{S}_{i}$ from a uniform discrete distribution of $q$ values for the Potts model, and from a uniform distribution on $[0,2\pi)$ for the XY model. In these simulations we use $N_{T}=10^{4}$ MCS for thermalization and $N_{A}=10^{5}$ MCS for the averaging.

The energy dependence of the acceptance rate for the Metropolis algorithm is summarized in Fig. 2(a). In general, the dependence turns out to be a non-linear featureless curve. The maximum difference is observed to be between the 3-state Potts and Ising models. The 4-state Potts model is closer to the Ising result, and for the XY model, the acceptance rate approaches that for the Ising model at very large temperatures, $T\gg J$. Results for the Glauber updates turn out to be qualitatively similar, and we present them for one spatial dimension in Fig. 3(a).

IV.2 Second moments of energy and acceptance rate

Mean values of energy and acceptance rate are computed as first moments of samples generated by the MC process. Given a one-to-one, monotonic relation between the mean values, it is instructive to compare the second moments. In general, the second moment of the reduced energy density $\varepsilon=E/JV=\langle H\rangle/JV$ is related to the specific-heat capacity,

[TABLE]

where $\langle\cdots\rangle$ stands for the average over the states generated by the MC process.

The variance of the acceptance rate is readily computed as the variance of a Bernoulli process of binary decisions ($1$ if an elementary update is accepted, and [math] otherwise). For a Bernoulli process, the variance, $\operatorname{var}R$, is related to the mean value, $R$, via $\operatorname{var}R=R(1-R)$. In the case of the 1D Ising model, given the exact results for $R$ derived above, hence also $\operatorname{var}R$ is known exactly.

Figure 2(b) displays the relation between the heat capacity and the variance of the acceptance rate for the Metropolis algorithm. We scale the variance of $R$ by $T^{2}$ in accordance with Eq. (18). At lowest temperatures, $T\ll J$, the heat capacity as a function of temperature has a maximum due to the well-known Schottky anomaly Kubo . Because of this, the curves in Fig. 2 form arcs for $C\sim 1$. Outside of this range, for $T>J$, the relation between second moments divided by $T^{2}$ is close to linear on the log-log scale for all one-dimensional models. Results of simulations using the heat-bath respectively Glauber updates, shown in Fig. 3(b), look qualitatively similar.

V Simulation results in 2D

It is instructive to compare the behavior of 1D models to higher dimensions, if only to see whether our observations are specific to 1D or have broader applicability. Specifically, 2D models with Hamiltonians (1)–(3) undergo a phase transition at a certain temperature $T_{c}$, between a high-temperature paramagnetic behavior and a low-temperature phase. For the Potts models, critical parameters are known analytically Wu , $T_{c}/J=1/\ln\left(1+\sqrt{q}\right)$, and for the Kosterlitz-Thouless transition of the XY model, MC simulations of Ref. Weber quote $T_{c}/J=0.887(2)$. Here the specific-heat capacity stays finite everywhere, with a smooth peak located about $20\%$ above $T_{c}/J$. The behavior of the MC process with local updates varies significantly between the paramagnetic phase ($T>T_{c}$), the low-temperature phase ($T<T_{c}$), and the critical region ($T\approx T_{c}$) BinderLandau .

Figure 4 shows results of simulations with the Metropolis algorithm of the two-dimensional models with $64^{2}$ spins for temperatures between $T/J=0.5$ to $10$. Here we use $N_{T}=10^{5}$ MCS for thermalization and $N_{A}=10^{6}$ MCS for averaging.

The first moments of the acceptance rate and energy shown in Fig. 4(a) are both smooth across the phase transition, with a relation which is close to linear in the critical region $T\sim T_{c}$. The second moments in Fig. 4(b) show two clearly separate branches, for $T>T_{c}$ and $T<T_{c}$, which join at the critical point. Note that while the heat capacity develops a singularity as $T\to T_{c}$, the variance of the acceptance rate remains smooth and does not show any signs of divergence. This can be explained by recalling that the heat capacity (18) is by definition proportional to the variance of the nonlocal total energy, while the acceptance rate and its variance refer to local measurements. It is also worth noting that the relative position of the low- and high-temperature branches for the Ising model differs from both Potts models and the XY model.

For heat-bath respectively Glauber updates, the results are qualitatively similar to those using Metropolis updates, see Fig. 5.

We have also verified that simulations of the three- and four-dimensional models behave qualitatively similar to the two-dimensional models. These results will be detailed elsewhere.

VI Conclusion

Concluding, we start with the observation that the acceptance rate of an update proposal of a Monte Carlo simulation can be regarded on par with thermodynamic functions of a model, and thus can itself be considered a thermodynamic function of a model under a given Monte Carlo dynamics.

For the one-dimensional Ising model we derive analytically that the mean value of the acceptance rate for local Metropolis updates is a linear function of energy. This linear dependence turns out to be specific for this combination of the updating scheme and the model: changing the updating algorithm to heat-bath updates changes the functional form of the relation between the mean values of the acceptance rate and energy, so that the relation is only linear away from the high-temperature region $T\gg J$.

We simulate several classical models in one and two spatial dimensions, the Ising model, the 3- and 4-state Potts models, and the XY model, and compute the dependence of the first and second moments of the acceptance rate on the mean value and the second moment of energy. We find that in general the relation is not linear in a wide range of temperatures, but is close to linear in the critical region around the transition temperature $T_{c}$.

Our result for the acceptance rate of the heat-bath algorithm for the one-dimensional Ising model can be viewed as an addition to the Glauber paper Glauber on the dynamics of the one-dimensional Ising model – since in that case, the heat-bath algorithm exactly reproduces the Glauber dynamics. The acceptance rate in the Glauber dynamics Glauber is the frequency of the spin flips, and it is given by Eqs. (16) and (17).

The acceptance rate can also be calculated analytically for any exactly solvable model (e.g., one-dimensional $q$-state Potts models and the two-dimensional Ising model), although this is not straightforward in all cases. One can compute the acceptance rate for any local Monte Carlo algorithm. In fact, we checked that for all models and for all dimensions for which we have code on hand, the acceptance rate is linear in the energy close to a second-order phase transition, as demonstrated for some two-dimensional models in Section V. Moreover, additional simulations show that in the vicinity of a first-order phase transition, the acceptance rate is a linear function of energy as well.

The source codes for our Monte Carlo simulations and data analysis are available from Ref. repo .

Acknowledgements.

E.B. and M.G. acknowledge the support of the Academic Fund Program at the National Research University Higher School of Economics (HSE), grant No. 18-05-0024, and by the Russian Academic Excellence Project “5-100”. W.J. thanks the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for support through project No. 189 853 844 – SFB/TRR 102 (project B04). L.S. acknowledges the support from the program 0033-2018-0010 of Landau Institute for Theoretical Physics and thanks Travis S. Humble (Oak Ridge National Laboratory) for the valuable discussions of the high-performance computing and INCITE program.

Appendix: Acceptance rates of local Monte Carlo updates for the

one-dimensional Ising model

In this appendix we turn our attention to the mathematical details of calculating the expectations of the acceptance rates of local MC updates.

VI.1 Acceptance rate of the Metropolis algorithm

Using the bond representation (9), we note that flipping a spin $S_{i}$ flips the values of two bond charges, $Q_{i}$ and $Q_{i-1}$. The acceptance probabilities depend on the sum $Q\equiv Q_{i}+Q_{i-1}$: for $Q=0$ or $1$, the Metropolis update is always accepted, since $\Delta E\leqslant 0$. For $Q=2$, the update is accepted with probability $e^{-4\beta J}=x^{-2}$.

Denoting the expected value of the acceptance probability by $R$, the expected value of the rejection probability is then

[TABLE]

Here the factor $2l(2l-1)/L(L-1)$ counts the probability that, in a configuration with $\sum_{i}Q_{i}=2l$, for a randomly chosen site $i$ we have $Q=2$, i.e., $Q_{i}=Q_{i-1}=1$.

The sum entering Eq. (19) is readily computed by differentiating twice the binomial formula

[TABLE]

The result is

[TABLE]

which is Eq. (12).

We now compare Eq. (12) to the internal energy of the system. The energy is given in general by $E=-d\ln Z/d\beta=(1/Z)dZ/d\beta$. Using (11) and $dx/d\beta=2Jx$ (since $x=e^{2\beta J}$), we obtain from the product rule

[TABLE]

which simplifies in the thermodynamic limit $L\rightarrow\infty$ to Eq. (14).

Comparing (VI.1) with (12), one readily sees that

[TABLE]

or

[TABLE]

is true for all lattice sizes $L$, i.e., the relation between the acceptance rate $R$ for the Metropolis update and the reduced energy density $\varepsilon$ of the 1D Ising model does not depend on the length $L$ of the one-dimensional chain with periodic boundary conditions.

VI.2 Acceptance rate of the heat-bath algorithm

The expected value of the acceptance probability $R$ of the heat-bath update can be calculated similarly to Eqs. (19) and (12). Here we directly compute the acceptance probability: acceptance probability of an elementary update of the spin $S_{i}$ is again defined by $Q\equiv Q_{i}+Q_{i-1}$, and the analog of Eq. (19) is

[TABLE]

where the terms in brackets correspond to $Q=2$, $Q=0$ and $Q=1$, respectively. Differentiating the binomial formula, we obtain

[TABLE]

where $\kappa=(x-1)/(x+1)=e^{-1/\xi}<1$ with $\xi$ denoting the correlation length. In the thermodynamic limit $L\rightarrow\infty$, the second factor $\frac{1-\kappa^{L}}{1+\kappa^{L}}=\frac{1-e^{-L/\xi}}{1+e^{-L/\xi}}$ in (24) approaches unity exponentially fast, and comparing to (14), we obtain Eq. (17).