# Optimized energy calculation in lattice systems with long-range interactions

## I Introduction

The numerical study of systems with long-range interactions is notoriously difficult, due to the large number of interactions that has to be taken into account. Specifically, the number of operations to calculate the energy of a single particle scales as the total number of particles in the system, in contrast to the case of short-range interactions, where the corresponding number of operations is of order unity. This implies that Monte Carlo-based methods are generally restricted to very small system sizes, which are still hampered by strong finite-size effects. Some years ago, this problem was resolved for the case of $`\mathrm{O}(n)`$ spin systems with (ferromagnetic) long-range interactions, for which a dedicated cluster algorithm was developed . Since the efficiency of this algorithm is independent of the number of interactions per spin, speed improvements of several orders of magnitude could be obtained compared to a conventional cluster algorithm. This speed-up pertains to the generation of independent configurations, for which the calculation of the energy is not required. Indeed, a variety of interesting physical results could be obtained by means of this method, see, e.g., Refs. . Whenever one needs to sample the internal energy, however, the improvement is much less dramatic: the major remaining advantage is that one only has to calculate the energy for truly independent configurations, rather than in every Monte Carlo step. Whereas this still implies that one can study systems which are an order of magnitude larger than those that can be accessed via Metropolis-type simulations (cf. Ref. ), one is eventually limited by the fact that the total computing time scales quadratically with the system size. One major disadvantage of this inaccessibility of the energy is the fact that it is not possible to apply histogram interpolations in order to obtain information on thermodynamic quantities over a large parameter space . In this paper, we point out that, for systems with periodic boundary conditions, this problem can be circumvented by calculating the internal energy in momentum space. Thus, one can apply a Fast Fourier Transform (FFT), reducing the total computational effort to $`𝒪(N\mathrm{log}N)`$ for a system containing $`N`$ spins. Indeed, this is a natural choice if one recognizes that the total energy is just given by a (discrete) convolution, which is one of the major applications of the FFT. Remarkably, the computational overhead entailed by the FFT turns out not to be a limiting factor: already for very small systems it is more efficient than a direct calculation of the energy.

The remainder of this paper is organized as follows. First, we derive an expression for the energy in terms of the Fourier-transformed spin system and point out that several other observables can be obtained on the fly, at negligible additional cost. We also give a detailed comparison of our approach and the conventional method. Next, we illustrate our approach by means of several new physical results for one-dimensional systems with long-range interactions. We end with a summary of our results.

## II Energy calculation

We will now first illustrate our approach for a $`(d=1)`$-dimensional system with an $`n`$-component order parameter, i.e., a generalized $`\mathrm{O}(n)`$ spin chain. This system is described by the Hamiltonian

$$=\frac{1}{2}\underset{x=1}{\overset{N}{}}\underset{y=1}{\overset{N}{}}J(xy)𝐒(x)𝐒(y),$$
(1)   
where the spins $`𝐒(x)`$ are $`n`$-component unit vectors and $`N`$ is the system size. The interaction $`J(x)`$ is defined for all $`x`$. Under the condition that periodic boundary conditions are employed, the effective coupling $`\stackrel{~}{J}(x)`$ between two spins is given by the sum over all periodic copies,

$$\stackrel{~}{J}(x)\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}J(x+mN)$$
(2)   
and hence has a period $`N`$. We set the self energy $`\stackrel{~}{J}(mN)`$, which is just an additive constant in the total energy, equal to zero. Each component of the spin configuration $`𝐒(x)`$ and the interaction $`\stackrel{~}{J}(x)`$ can then be written as a Fourier sum

$$f(x)=\frac{1}{\sqrt{N}}\underset{k=0}{\overset{N1}{}}f_ke^{2\pi ikx/N},$$
(3)   
where the Fourier coefficients $`f_k`$ are obtained from the discrete Fourier transform of $`f(x)`$. By means of the discrete convolution theorem  it is then straightforward to show that Eq. (1) can be written as

$$=\frac{\sqrt{N}}{2}\underset{k=0}{\overset{N1}{}}\stackrel{~}{J}_k𝐒_k𝐒_k.$$
(4)   
The essential step is now, that application of the Fast Fourier Transform  reduces the computational effort for the calculation of the $`N`$ Fourier coefficients from $`𝒪(N^2)`$ to $`𝒪(N\mathrm{log}N)`$, thus, in principle, greatly speeding up the calculation of the total energy. The sum in (4) adds another $`𝒪(N)`$ operations, but this is compensated for by the fact that one typically also wants to calculate the magnetic susceptibility $`N^1|_{x=0}^{N1}𝐒(x)|^2`$, which in the momentum-space representation is immediately given by $`|𝐒_{k=0}|^2`$. For maximum efficiency, one has to restrict the system size to powers of 2. Naturally, the calculation of the coefficients $`\stackrel{~}{J}_k`$ has to be carried out only once.

Even more can be gained, if one also desires to calculate the spin–spin correlation function $`g(r)=𝐒(0)𝐒(r)=N^1_{x=0}^{N1}𝐒(x)𝐒(x+r)`$. The discrete correlation theorem  states that the Fourier transform $`g_k`$ of $`g(r)`$ is equal to $`N^{1/2}𝐒_k𝐒_k`$, so that $`g_k`$ is obtained by $`N`$ multiplications rather than another $`𝒪(N^2)`$ operations in the real-space representation.

All the above estimates are only measures for the complexity of the algorithm, which become valid for sufficiently large $`N`$. It remains to be seen whether the FFT-based approach is actually faster for the range of system sizes that can be accessed in present-day Monte Carlo simulations, which for lattice models go up to $`N10^510^6`$. Figure 1 compares the required CPU time per spin for the calculation of the internal energy via Eq. (1) and Eq. (4), respectively, and the susceptibility. As expected, the former method scales asymptotically linearly with $`N`$. For the latter method, two estimates are given, which only differ in the choice of the FFT implementation. The slower results (open squares) were obtained by means of the routines of Ref.  and the faster (triangles) by means of those of Ref. . Although these two estimates differ by as much as a factor of $`5`$, both of them outperform the conventional method already for $`N10`$; for $`N=2^{18}`$ the improvement amounts to roughly four orders of magnitude. The initial downward trend in Fig. 1 is due to overhead being distributed over an increasing number of spins. Likewise can the “irregularities” in the FFT estimates be attributed to computational aspects. The slight deviations from linearity in the conventional estimates, however, are due to statistical inaccuracies in the timing: for $`N𝒪(10^4)`$ this method becomes already prohibitively slow. Thus, we conclude that the method presented here provides a very efficient approach to energy calculations in lattice systems with long-range interactions; while there is still a weak system-size dependence in the required computational effort per spin, this no longer constitutes a bottleneck for practical applications. Note that higher-dimensional models can also be treated in this fashion.

## III Applications

### A The Ising chain

As a first example, we consider the Ising chain with algebraically decaying interactions, $`J(x)J|x|^{1\sigma }`$. The critical behavior for this system is essentially classical for $`0<\sigma \frac{1}{2}`$ and nonclassical for $`\frac{1}{2}<\sigma 1`$, see Ref. . Numerical results for the thermal exponent $`y_t`$  have indicated that the latter regime can be subdivided into two parts: $`y_t>\frac{1}{2}`$ for $`\frac{1}{2}<\sigma 0.65`$ and $`y_t<\frac{1}{2}`$ for interactions that decay faster. This implies that the specific heat only diverges in a part of the nonclassical regime and should display a cusp-like singularity in the remaining part. By means of illustration, we have calculated the specific heat for $`\sigma =0.25`$ and $`\sigma =0.90`$. In both cases, we expect to find a function that does not diverge at the critical point, although the behavior should be qualitatively different. Simulations were carried out for $`N=2^p`$, with $`3p16`$, at a number of different couplings, for several times $`10^6`$ independent samples per system size. The full curves were determined by means of the multiple-histogram method , where great care was taken to minimize systematic errors due to the histogram interpolation.

Figure 2 shows the specific heat $`C`$ for $`\sigma =0.25`$, as a function of the reduced coupling $`KJ/(k_\mathrm{B}T)`$. It displays several close similarities to the specific heat of the mean-field model, including the build-up of a jump discontinuity at the critical point, the crossing of the finite-size curves in a single point at $`K_c`$ (up to corrections to scaling), and (not visible on this scale) an excess peak in the curves for finite systems, i.e., $`lim_N\mathrm{}C_{\mathrm{max}}(N)lim_{KK_c}lim_N\mathrm{}C(K,N)`$ . As shown in Ref. , the location of the specific-heat maximum shifts as a function of system size, according to

$$K_{\mathrm{max}}=K_c+a_1L^{y_t^{}}+a_2L^{2y_t^{}}+b_1L^{\sigma 1}+\mathrm{},$$
(5)   
where $`y_t^{}=\frac{1}{2}`$ and the coefficients $`a_i,b_i`$ are nonuniversal. A fit to this expression yielded $`y_t^{}=0.51(6)`$ and $`K_c=0.1147(5)`$, in good agreement with $`K_c=0.114142(2)`$ . The inset shows the peak height as a function of system size, strongly suggesting that the maximum is indeed finite in the thermodynamic limit.

The case $`\sigma =0.90`$, shown in Fig. 3, clearly exhibits a distinctly different behavior. The specific heat is now nonzero in the thermodynamic limit, on either side of the critical point, and indeed displays the expected cusp-like singularity. The inset confirms that the maximum is convergent for $`N\mathrm{}`$. Since $`y_t`$ is still sufficiently close to $`\frac{1}{2}`$, i.e., the absolute value of the exponent $`\alpha `$ is sufficiently small, the location of the maximum cannot be distinguished from the critical point, unlike the case $`\sigma =1`$, where it is expected to occur at a coupling $`K<K_c`$.

### B The three-state Potts chain

The ferromagnetic Potts model provides a particular generalization of the Ising model with respect to the number $`q`$ of possible coexisting ordered phases. For $`q=2`$ the Ising model is recovered; for $`q>2`$ the ferromagnetic Potts model defines a genuine universality class distinct from the Ising and the more general $`\mathrm{O}(n)`$ universality class. The Potts model is of particular theoretical interest, because the phase transition it describes may be of first or second order depending on $`q`$ and the spatial dimension $`d`$, even in the absence of symmetry-breaking fields. For nearest-neighbor interactions in $`d=2`$ many properties of the Potts model are exactly known . In particular, if the model is in the first-order regime and $`q`$ is sufficiently large, the asymptotic finite-size properties of the nearest-neighbor ferromagnetic Potts model have been established in a rigorous fashion . For long-range interactions, however, much less is known and currently available numerical data are limited to rather small systems . We demonstrate in the following, that also for the Potts model the cluster algorithm introduced in Ref.  can be combined with the FFT, allowing the numerical treatment of much larger systems.

Again, we concentrate on the case of algebraically decaying interactions $`J(x)|x|^{1\sigma }`$. The Hamiltonian of the ferromagnetic Potts chain with periodic boundary conditions can then be written in the same form as Eq. (1), where the Potts spins $`𝐒(x)`$ are unit vectors which mark the corners of a (hyper)tetrahedron in $`q1`$ dimensions. For the present case $`q=3`$ we employ the complex notation

$$𝐒(x)𝒮(x)\{1,e^{2\pi i/3},e^{4\pi i/3}\},$$
(6)   
i.e., $`𝐒(x)𝐒(y)=\mathrm{}[𝒮(x)𝒮(y)^{}]`$, where the asterisk denotes the complex conjugate. The spin representation of the Potts model given by Eq. (1) is equivalent to the standard Kronecker representation, but it has the advantage that the configurational energy is directly accessible by means of the FFT. According to mean-field theory the ferromagnetic Potts model should always show a first-order phase transition for $`q>2`$. For our case of $`d=1`$ and algebraically decaying interactions, one therefore expects the mean-field prediction to be correct for sufficiently small values $`\sigma >0`$ of the decay exponent of the interaction, i.e., there should be a critical value $`\sigma _c`$ separating first- and second-order behavior. Mean-field theory provides an important guideline for the interpretation of our Monte Carlo data, so we briefly summarize the basic mean-field predictions. Following Ref.  we introduce the probability $`p_\kappa (x)`$ that lattice site $`x`$ is occupied by the Potts state $`\kappa `$, $`1\kappa q`$, and we define a homogeneous scalar order parameter $`s`$ indicating a broken symmetry with respect to Potts state $`\kappa =1`$:

$$p_1(x)m_1=\frac{1+(q1)s}{q},p_\kappa (x)m_\kappa =\frac{1s}{q},\mathrm{\hspace{0.33em}2}\kappa q.$$
(7)   
For a given value of $`s`$ the mean-field free-energy density $`f_{\mathrm{MF}}(s)`$ (in units of $`k_\mathrm{B}T`$) is then obtained as

$$f_{\mathrm{MF}}(s)=K\zeta (1+\sigma )s^2+\left\{[1+(q1)s]\mathrm{log}[1+(q1)s]+(q1)(1s)\mathrm{log}(1s)\right\}/q,$$
(8)   
where $`K=J/(k_\mathrm{B}T)`$ denotes the reduced coupling and $`\zeta (\alpha )`$ is the Riemann zeta function. Note that the replacement $`K(q1)q^1K`$ transforms Eq. (8) from the spin representation into the Kronecker representation. The transition point $`K=K_{\mathrm{MF}}^t`$ from the disordered phase $`s=0`$ to the ordered $`(\kappa =1)`$ phase $`s=s_{\mathrm{MF}}^t`$ follows from standard mean-field arguments :

$$K_{\mathrm{MF}}^t\zeta (1+\sigma )=\frac{(q1)^2}{q(q2)}\mathrm{log}(q1),s_{\mathrm{MF}}^t=\frac{q2}{q1}.$$
(9)   
According to Eqs. (7) and (9) the distribution function $`P(m_1)`$ for a finite system displays three maxima near the transition temperature $`T_{\mathrm{MF}}^tJ/(k_\mathrm{B}K_{\mathrm{MF}}^t)`$: one at $`m_1=1/q`$ for the disordered phase, one at $`m_1=(q1)/q`$ for the ordered $`(\kappa =1)`$ phase, and one at $`m_1=1/[q(q1)]`$ for the ordered phases with respect to the remaining Potts states $`(\kappa 2)`$. Note that all ordered phases appear with equal probability in the course of the simulation. For $`\sigma 0.4`$ and in our case of $`q=3`$ these three peaks in $`P(m_1)`$ are indeed located very close to their mean-field positions. For $`\sigma =0.6`$ the peaks are still clearly separated, but they occur at positions shifted with respect to the mean-field predictions and for $`\sigma 0.7`$ the peaks start to overlap strongly and can only be identified for very large systems (see below).

Although the algorithm introduced in Ref.  is by far the most efficient one for the simulation with spin systems with long-range interactions, it is not able to deal with first-order phase transitions beyond a certain system size. The reason is that like the Metropolis algorithm the Wolff cluster algorithm encounters an activation barrier between states with and without long-range order, which is set by the energy-density gap between the disordered and the ordered phase. For a given size of the gap the tunneling time between disordered and ordered phases, and therefore the required sampling time, increases exponentially with the system size so that the attainable system size $`N`$ is severely limited. This tunneling problem can be solved by employing the well-established ideas of multicanonical sampling , however, the generalization of the cluster algorithm  to an efficient multicanonical algorithm is beyond the scope of the present article.

The data we present in the following have been obtained from histograms of the energy taken at several temperatures. The data are again conveniently analyzed by the optimized multiple-histogram method . For the values $`\sigma =0.2`$ and $`\sigma =0.4`$ the Potts chain undergoes a strong first-order phase transition  which limits the chain length to $`N=2^{13}`$ in the former case and to $`N=2^{14}`$ in the latter case. We reinvestigate chains from $`N=2^{10}`$ spins to the respective maximum chain length by taking a few times $`10^6`$ independent samples for each system size and temperature, where a comparison with the finite-size theory of Borgs et al.  turns out to be very instructive. We restrict the detailed presentation to the case $`\sigma =0.4`$; the data for $`\sigma =0.2`$ are qualitatively very similar. Near the transition temperature the energy distribution function $`P(E)`$ displays two peaks characterizing the ordered and the disordered phase, respectively . In Ref.  the temperature of equal peak height is taken as an estimate for the transition temperature on the finite system, where the leading finite-size corrections are of the order $`𝒪(1/N)`$. If the systems are large enough, i.e., the peaks in the energy distribution are well separated, the ratio $`W_o/W_d`$ of the weight of the ordered phase $`W_o`$ and the weight of the disordered phase $`W_d`$ provides a far more convenient indicator of the transition temperature, because the associated finite-size corrections decay exponentially with $`N`$ . Our result for long-range interactions is shown in Fig. 4. The transition temperature is marked by the intersection of $`W_o/W_d`$ as a function of temperature for the three largest systems. For $`N=2^{11}`$ the peaks in the energy distribution are not well separated so that $`W_o/W_d`$ is not well defined in this case. For $`N2^{12}`$ the curves meet at $`W_o/W_d=1.67(2)`$ as shown by the solid line. For $`\sigma =0.2`$ we find a corresponding intersection at $`W_o/W_d=1.25(2)`$. For nearest-neighbor interactions in $`d2`$ and sufficiently large $`q`$ the value $`W_o/W_d=q`$ is expected to indicate the transition temperature . Surprisingly, we find a much smaller value here which appears to increase with $`\sigma `$. From Ref.  one furthermore expects that the curves displaying the energy density for different system sizes as a function of temperature exhibit an intersection close to the transition temperature, where the deviations are predicted to be exponentially small in $`N`$. In Fig. 5 this situation is shown for long-range interactions with $`\sigma =0.4`$. A corresponding result has been obtained for $`\sigma =0.2`$. The energy densities intersect near the transition temperature found in Fig. 4, where the shifts between mutual intersections seem to be compatible with exponentially small finite-size effects. Still too few data are available for a quantitative analysis of these shifts, but finite-size effects of the order $`1/N`$ can be ruled out. The fourth-order energy cumulant $`U_4`$ defined by

$$U_4^4/^2^2$$
(10)   
is shown in Fig. 6 for different system sizes as a function of temperature, where $``$ is the Hamiltonian given by Eq. (1). These cumulants should also intersect at the transition temperature in the limit $`N\mathrm{}`$. The data displayed in Fig. 6 show the expected tendency, but the finite-size corrections are much larger than those for the weight ratio or the energy density (see Figs. 4 and 5, respectively). For $`\sigma =0.2`$ a corresponding result has been found. The systematic shift of the intersections of $`U_4`$ for different system sizes is compatible with a $`1/N`$ behavior as anticipated in Ref.  for nearest-neighbor interactions, but the present amount of data is too limited to give reliable quantitative evidence for this behavior.

For nearest-neighbor interactions and periodic boundary conditions the energy density asymptotically obeys the scaling law (cf. Eq. (1) of Ref. )

$$E(\beta ,L)\frac{E_d+E_o}{2}\frac{E_dE_o}{2}\mathrm{tanh}\left[\frac{E_dE_o}{2}(\beta \beta _t)N+\frac{\mathrm{log}q}{2}\right],$$
(11)   
where $`\beta =1/(k_\mathrm{B}T)`$ is the inverse temperature and $`\beta _t`$ is the transition point. It is instructive to compare the scaling form given by Eq. (11) with the data displayed in Fig. 5. The energies of the disordered phase $`E_d`$ and the ordered phase $`E_o`$ can be read off from the positions of the two maxima of the energy distribution function. It turns out that the data in Fig. 5 and their counterpart for $`\sigma =0.2`$ are in fact consistent with Eq. (11) within the error bars, provided the number of states $`q`$ on the right-hand side of Eq. (11) is replaced by the effective value $`q_{\mathrm{eff}}(\sigma =0.4)W_o/W_d|_{\beta =\beta _t}=1.67`$ measured in Fig. 4 at the transition temperature or its counterpart $`q_{\mathrm{eff}}(\sigma =0.2)=1.25`$, respectively. The quantitative comparison of our data for $`\sigma =0.2`$, $`N=2^{13}`$ and $`\sigma =0.4`$, $`N=2^{14}`$ with Eq. (11) is shown in Fig. 7. Within the statistical errors, the agreement is excellent except for larger values of the scaling variable $`(\beta \beta _t)N`$. These deviations are due to the fact that Eq. (11) only holds asymptotically for sufficiently large systems. For finite systems additional finite-size corrections enter through the residual $`N`$-dependence of $`E_d`$, $`E_o`$, and $`q_{\mathrm{eff}}(\sigma )`$ which appear as parameters in Eq. (11). Figure 7 demonstrates that the finite-size effects in the three-state Potts chain with periodic boundary conditions and long-range interactions can be interpreted in terms of the Borgs–Kotecký theory  for the nearest-neighbor Potts model in higher dimensions for an effective number of states $`q_{\mathrm{eff}}(\sigma )`$. The physical meaning of $`q_{\mathrm{eff}}(\sigma )`$, however, remains unclear. The proof of Eq. (11) also requires the assumption that $`q`$ is sufficiently large , so $`q=3`$ may not be sufficient for a quantitative comparison. On the other hand, numerical investigations have shown that the $`q5`$ nearest-neighbor Potts model in $`d=2`$  and the $`q=3`$ nearest-neighbor Potts model in $`d=3`$  follow the theory of Ref.  very closely despite the small values of $`q`$. Further analytical and numerical studies are required to settle this question.

We close our discussion of the three-state Potts chain with a brief summary of our results for $`\sigma =0.6`$, 0.7, and 0.75, which have been studied with reduced statistics ($`10^5`$ independent samples for each system size and temperature). The value $`\sigma =0.6`$ is still located in the first-order regime , but in order to obtain a well-defined weight ratio $`W_o/W_d`$ system sizes of $`N2^{16}`$ Potts spins are required, although the maxima in the energy distribution are well separated already for $`N2^{14}`$. We have performed simulations for $`N=2^{14}`$ up to $`N=2^{17}`$ at four to six temperatures for each system size. Even for $`N=2^{17}`$ the finite-size corrections are too large to identify an intersection of the energy densities as accurately as displayed in Fig. 5. Data acquisition for $`N>2^{17}`$ is strongly hampered by the energy gap so that we refrain from discussing our data for $`\sigma =0.6`$ in any more detail. In contrast to Ref.  the value $`\sigma =0.7`$ of the decay exponent can undoubtedly be identified as a member of the first-order regime. For system sizes $`N2^{16}`$ the energy distribution function displays the typical double-peak structure which becomes sharper as the system size is increased at fixed temperature. We illustrate this in Fig. 8, where the data for $`P(E)`$ are shown at $`T/T_{\mathrm{MF}}^t=0.8095`$, which is close to the transition point. The same analysis has been repeated for $`\sigma =0.75`$ and system sizes up to $`N=2^{19}`$ spins. Although $`P(E)`$ also develops a plateau similar to the one displayed in Fig. 8 for $`N=2^{15}`$, no double-peak structure could be resolved up to $`N=2^{19}`$ so that $`\sigma =0.75`$ may already belong to the second-order regime of the three-state Potts chain with long-range interactions. However, the detection of a double-peak structure in $`P(E)`$ for a given value of $`\sigma `$ is essentially a matter of attainable system size (see Fig. 8), so $`\sigma _c>0.7`$ is the only safe conclusion here.

## IV Conclusions

The combination of the recently developed cluster algorithm  for systems with long-range interactions with the Fast Fourier Transform for the calculation of the configurational energy leads to a Monte-Carlo algorithm with a very high efficiency. In particular, the FFT allows to extend the attainable system sizes by two orders of magnitude in comparison with other approaches (cf. Ref. ). Histogram interpolation methods then allow the investigation of thermodynamic properties of these systems with unprecedented resolution. By construction the algorithm can only deal with first-order phase transitions up to a limited system size. In order to avoid this limitation the algorithm must be generalized to include multicanonical sampling. Here, we have demonstrated the potential of the algorithm for the Ising chain and the three-state Potts chain with algebraically decaying interactions. For completeness, it is mentioned that also for system sizes that are not integer powers of two a considerable gain can be obtained by performing the discrete Fourier Transform via, e.g., a prime-factor algorithm.

For the Ising chain we have investigated the finite-size behavior of the specific heat in the classical regime for $`\sigma =0.25`$ and in the nonclassical regime for $`\sigma =0.9`$. In the former case the specific heat behaves essentially mean-field like, i.e., the expected discontinuity in the specific heat at the critical temperature in the thermodynamic limit builds up as the system size is increased. On the other hand, the choice $`\sigma =0.9`$ is expected to yield a negative specific-heat exponent, i.e., a cusp singularity should appear with increasing system size. Our numerical data confirm this behavior as well and clearly show the different shapes of the specific-heat curves in the two cases.

The three-state Potts chain is expected to show a first-order phase transition for $`\sigma <\sigma _c`$, where our results indicate that $`\sigma _c>0.7`$. For $`\sigma =0.2`$ and $`\sigma =0.4`$, for which the $`q=3`$ Potts chain displays a strong first-order phase transition, our data confirm the Borgs–Kotecký scenario of the first-order transition in Potts models with nearest-neighbor interactions in higher dimensions, provided the number $`q`$ of states is replaced by the effective number of states $`q_{\mathrm{eff}}(\sigma )=W_o/W_d|_{\beta =\beta _t}<q`$ which also enters the finite-size scaling form of the energy density near the transition temperature. For $`\sigma =0.6`$ the same behavior can be confirmed only on a semi-quantitative level, because much larger systems must be investigated in order to obtain sufficient resolution. The mechanism that leads to the reduction of the effective number of states and the physical interpretation of $`q_{\mathrm{eff}}(\sigma )`$ are not known.

###### Acknowledgements.

We gratefully acknowledge helpful discussions with T. Neuhaus and W. Janke, and stimulating comments by K. Binder. Furthermore, the authors wish to thank the organizers of the Twelfth Annual Workshop on Recent Developments in Computer Simulation Studies (Athens, Georgia), where the collaboration leading to this work was initiated. M. Krech also gratefully acknowledges financial support through the Heisenberg program of the Deutsche Forschungsgemeinschaft.