Two-size Probability-Changing Cluster Algorithm

We propose a self-adapted Monte Carlo approach to automatically determine the critical temperature by simulating two systems with different sizes at the same temperature. The temperature is increased or decreased by checking the short-time average of the correlation ratios of the two system sizes. The critical temperature is achieved using the negative feedback mechanism, and the thermal average near the critical temperature can be calculated precisely. The proposed approach is a general method to treat second-order phase transition, first-order phase transition, and Berezinskii-Kosterlitz-Thouless transition on the equal footing.


Introduction
Finite-size scaling (FSS) [1] is a basic concept in the study of phase transitions and critical phenomena. The Binder ratio [2], essentially the moment ratio, is widely used in the analysis of the numerical data. The moment ratios of the magnetization m, for different sizes scale as that the same procedure is also effective for studying the Berezinskii-Kosterlitz-Thouless (BKT) transition, where a fixed line instead of a fixed point exists. We select the 2D 5-state clock model that has two BKT transitions with higher and lower transition temperatures.
The remaining part of the paper is organized as follows: We describe the twosize PCC algorithm for the second-order transition in Section II. Two-size PCC studies of the first-order transition and BKT transition are discussed in Sections III and IV, respectively. Section V is devoted to summarizing the study and discussing results.

Second-order transition
Let us start with the second-order phase transition. As an example, we consider a 2D Ising model on the square lattice, whose Hamiltonian is given by The summation is taken over the nearest-neighbor pairs, and periodic boundary conditions are imposed in numerical simulations. Using the Wolff single-cluster flip algorithm for spin update, we simulate the two system sizes simultaneously. In determining the critical point (line), we use the ratio of the correlation functions with different distances, i.e., the correlation ratio [4], Here, g(r) is a correlation function with the distance r. For the values of r and r ′ , we choose r = L/2 and r ′ = L/4 for numerical calculations. The actual procedure for the two-size PCC algorithm is as follows. We use two systems of different sizes, say L and L/2. After simulating some steps at the same temperature, we measure the correlation ratios of both the systems, R(L) and R(L/2). We increase or decrease the inverse temperature β (= 1/T in units of the coupling J) according to the following rule: where ∆β > 0. There are two parameters to choose, the number of Monte Carlo steps (MCS) for taking a short-time average, N av , and the difference of β, ∆β. Note that the cluster flip algorithm is effective because the rapid equilibration is required after a change in temperature. The plot of the time evolution of β is shown in Fig. 1. For this plot, we chose N av = 4000 and ∆β = 0.00005; that is, after every 4000 MCS, β is changed by ±0.00005. We will discuss the choice of N av and ∆β later. In the figure, the time steps are given in units of 1000 MCS, and we show the data up to 4 × 10 6 MCS. The system sizes are L = 64 and L = 512; that is, the set of system sizes are (64,32) and (512,256). We see that the temperature oscillates around the average value. The width of fluctuation decreases as the system size increases because of the effect of self-averaging. For convenience, we denote the exact value of β c (=ln(1 + √ 2)/2 = 0.4406868) for an infinite system by a dotted line.
Next, we examine the histogram of β, h(β), for the two-size PCC algorithm. We plot h(β) of the 2D Ising model for L = 64, 96, 128, 192, 256, 384, and 512 in Fig. 2. Measurement is performed for 4 × 10 6 MCS after equilibration of 10,000 MCS. We made 32 runs for each system size in order to estimate statistical errors. For smaller sizes (L = 64, 96, and 128), we made 64 runs. The parameters N av and ∆β were chosen as 4000 and 0.00005, respectively. In the plot, the histogram h(β) is normalized by The obtained histogram is similar to normal distribution, and we see that the histogram becomes sharper with an increase in system size. Furthermore, the peak position approaches the exact value of β c for the infinite system.
The transition (inverse) temperature for each size, β c (L), was estimated using the averaged value of β. The plot of the size dependence of β c (L) as a function of 1/L is shown in Fig. 4, where the statistical errors were estimated by 2σ of the distribution of β c . We can see from Fig. 4 that β c (L) rapidly approaches the exact value ln(1 + √ 2)/2 = 0.4406868, which is denoted by a dotted line, even for small sizes. The rapid convergence of the present algorithm is apparent when we compare the size-convergence rate of β c (L) with that of the original PCC ( Fig. 1 of Ref. [10]). In the original version of the PCC algorithm, although the size-dependent β c (L) is automatically tuned, we still have to consider the size dependence of β c (L) based on the FSS. Instead, with the two-size PCC algorithm, the infinite-size critical temperature is There is a single peak, which will be compared with the case of the first-order transition later. We see that the distribution becomes sharper as the system size increases. The peak position approaches the exact energy at the critical temperature, that is, The exact critical value is denoted by a dotted line in the figure.
Because each system wanders around temperature, we can take a thermal average of physical quantities at each temperature. This is the same situation encountered in the replica exchange method [11]. There, the temperatures of replicas are exchanged following the transition probabilities based on the Boltzmann weight; the thermal average at a fixed β is obtained by averaging over different replicas. Before showing the data of the thermal average of physical quantities, we present the energy distribution for a fixed value of β. The energy distribution is decomposed as The data of p(E/N; β) at two typical temperatures, β = 0.44085 and 0.4405, together with the whole distribution of p(E/N), is shown in Fig. 6. The system size is fixed at L = 128. The value of the β-decomposed distribution is magnified twenty times for clarity. Two energy distributions with different values of β are related to each other through the equation It is a reweighting of the Boltzmann factor, which is the basis of the histogram method by Ferrenberg and Swendsen [12]. The thermal average of a physical quantity A at β ′ is obtained by the measurement at β through the relation where [· · ·] stands for the Monte Carlo average at β. We note that the β-decomposed energy distribution, p(E/N; β), does not depend on h(β). Thus, we obtain the thermal average of a physical quantity A, A β , without considering N av and ∆β. As an example of the physical quantities, we show the correlation ratio R(β) as a function of temperature for various system sizes in Fig. 7. We see that the data of different sizes intersect at the critical point within statistical errors. The critical value of the correlation ratio, R c , is calculated as follows: using the Jacobi θ-functions [13] (see also Refs. [14,4]). This value is denoted by a dotted line in Fig. 7. Our simulation reproduces the exact value of R c with an accuracy up to four digits. Now let us consider the FSS. Because the critical temperature and the critical exponents are known for the 2D Ising model, R(β) are plotted as a function of (β − β c )L 1/ν , as shown in Fig. 8, where β c = ln(1 + √ 2)/2 and 1/ν = 1. We see that the FSS works quite well.

First-order phase transition
We now consider the case of the first-order phase transition. The 2D ferromagnetic q-state Potts model [15,16,17] is taken into account. The Hamiltonian is given by where δ ab is the Kronecker delta. This model is known to show the second-order phase transition for q ≤ 4 and first-order phase transition for q ≥ 5. Here, we provide the data for a two-size PCC calculation of the 2D 6-state Potts model. Hysteresis in the first-order transition systems should be considered, which is different from the conditions of the second-order transition. It is more feasible to employ the multi-cluster update of the Swendsen-Wang type [8] because a spin configuration changes extensively with such an update. For the systems with the second-order transition, there is no appreciable difference in the choice of the cluster update. The number of the steps for calculating the short-time average, N av , and the difference in β, ∆β, were chosen as 4000 and 0.00001, respectively. We used smaller values of ∆β as this could help reduce the effect of hysteresis. We conducted measurements for 4 × 10 6 steps after equilibration of 10,000 steps; such measurements were repeated 64 times for L = 64 and 96, and 32 times for L = 128, 192, 256, 384, and 512. We note that the correlation function of the q-state Potts model is given as The histogram of β, h(β), for the system sizes L = 64, 96, 128, 192, 256, 384, and 512 is shown in Fig. 9. The histogram exhibits sharp peaks, and the peak position gradually approaches the exact value. The exact value of the first-order transition inverse temperature for the infinite system is given by ln(1 + √ q) = 1.23823 (for q = 6), in units of J, which is shown by the dotted line.
The average value of β, β c (L), is plotted as a function of 1/L in Fig. 10. The exact value of the first-order transition inverse temperature (1.23823) is given by the dotted line. We can observe that the calculated estimate of the transition temperature converges to the exact value with five-digit accuracy.
The distribution of E, p(E/N), is shown in Fig. 11 for various sizes. We observe double peaks, which are specific to the first-order transition. In Fig. 12, we plot the β-decomposed energy distribution, p(E/N; β). The system size is set to be L = 128.  Here, we show the data of two typical temperatures, β = 1.23809 and 1.23794, which are on both sides of the peak value of h(β) shown in Fig. 9, together with the entire distribution of p(E/N). The value of the β-decomposed distribution is magnified twenty times for clarity. We can observe that the weight of high energies increases for the hightemperature (low-β) energy distribution. We now examine the peak positions of energy. Baxter [18] (see also [19]) calculated the exact difference in the higher energy peak E 2 and the lower energy peak E 1 , i.e., the latent heat. The exact result is where Θ = arcosh( √ q/2). The middle point (E 1 + E 2 )/2N is also given as Thus, for q = 6, E 1 /N and E 2 /N are calculated as 0.49102 and 0.69248, respectively. These values are given in Fig. 11. We can observe that the positions of the energy peaks approach the exact infinite values as the system size increases. The size dependences of the numerical estimates of E 1 /N and E 2 /N are plotted as a function of 1/L in Fig. 13. The statistical errors are within the size of marks. They converge to the exact values [18].

BKT transition
The 2D spin systems with continuous XY symmetry exhibit a unique phase transition called the Berezinskii-Kosterlitz-Thouless (BKT) transition [20,21,22]. There exists a BKT phase of a quasi long-range order (QLRO), where the correlation function decays as a power law. Here, we consider the q-state clock model, which is a discrete version of the classical XY model. The Hamiltonian is given by The 2D q-state clock model experiences the BKT transition for q ≥ 5, whereas the q = 4 clock model is two sets of the Ising model and the 3-state clock model is equivalent to the 3-state Potts model. The q = 2 clock model is simply the Ising model. For q ≥ 5, there is an interplay between the plane-rotator symmetry, which attempts to preserve the BKT phase, and the discreteness, which tends to create a long-range order (LRO) at low temperatures. Two transition temperatures, T 1 < T 2 , are observed; each corresponds to the transition between LRO and QLRO and between QLRO and a disordered phase.
We conducted a simulation of the two-size PCC algorithm for the 2D q = 5 clock model. As the system has a wide temperature range in the critical state, we selected a larger ∆β, 0.002. Again, we selected N av =4000. A histogram of β, h(β), of the q = 5 clock model is shown in Fig. 14. The data in the histogram are widely distributed, which is contrast to the case of strong transitions, i.e., the second-order and first-order transitions. This is related to the fact that there is a fixed line instead of a fixed point for the system with the BKT transition. In the figure, the numerical estimates of β 1 (1/T 1 = 1/0.911(5)) and β 2 (1/T 2 = 1/0.940(5)) [7] are shown by dotted lines for convenience.
The distribution of E, p(E/N), is shown in Fig. 15 for various sizes. In Fig. 16,   we plot the β-decomposed energy distribution p(E/N; β). The system size is L = 128. Data for three β's are shown; from left, β=1.10, 1.08, and 1.06. The value of the βdecomposed distribution is magnified ten times for clarity. Energy peaks are observed at certain values depending on the temperature.
We calculate the temperature dependence of physical quantities using the same procedure as the one followed by the Ising model. Although the histogram h(β) shown in Fig. 14 is not very smooth, fairly accurate estimates of the thermal average of physical quantities at a fixed β can be obtained as was discussed in the case of the Ising model.
The temperature dependence of the correlation ratio for various sizes is plotted in Fig. 17. The numerical estimates of β 1 and β 2 reported in [7] are shown by dotted lines for convenience. In the intermediate temperature range, the correlation ratios of different sizes take the same value, whereas they start to exhibit variations below β 2 and above β 1 . To locate the BKT transition temperatures precisely, a careful FSS treatment with exponential divergence behavior is required [7]. When the present method is used directly, the systems remain in the intermediate state for a long time. We may set windows for the allowed temperature range. In the case of the 2D q = 5 clock model, the temperature range may, for example, be restricted as β < 1.07 for the β 2 transition, and β > 1.09 for the β 1 transition.

Summary and discussion
In this paper, we described the two-size PCC algorithm. We simultaneously simulate two systems of different sizes at the same temperature. Comparing the short-time average of the correlation ratios of the two sizes, we increase or decrease the temperature based on the negative feedback mechanism given by Eq. (5). A temperature near the critical temperature is automatically selected. It is simply an Ehrenfest model for diffusion with a central force [23,24].
For the strong transitions including second-order or first-order transitions, the temperature peaks sharply at the critical temperature. Thus, we can locate the critical temperature in a self-adapted manner. The energy distribution is singly peaked in the case of the second-order transition, whereas it is doubly peaked in the case of the firstorder transition. As the system wanders around the temperature, we can calculate the thermal average of physical quantities for each temperature. We showed the results of the correlation ratios of the 2D Ising model, which demonstrated a satisfactory FSS behavior. For the first-order transition, we determined the double-peak positions E 1 and E 2 for the 2D q = 6 Potts model. The results were compared with the exact values obtained by Baxter [18].
In the case of the systems with the BKT transition, the temperature is widely distributed in the two-size PCC algorithm, which is owing to the existence of a fixed line in the BKT transition. An investigation of the temperature dependence of the correlation ratio for the 2D q = 5 clock model showed that correlation ratios of different sizes take the same value in the intermediate BKT state, whereas they start to vary below β 2 and above β 1 . We can easily determine the specific behavior of the BKT transition compared to the second-order transition or the first-order transition.
The advantage of the present method is that a temperature range can be automatically selected near the critical temperature. As the sampling of the temperature peaks at the critical temperature, the critical phenomena can be studied efficiently.
To summarize, we have proposed a unified method of numerical simulation that can treat the second-order phase transition, the first-order phase transition, and the BKT transition with equal footing. By simultaneously simulating two systems of different sizes, say L and L/2, we could measure the correlation functions, which are essential when investigating phase transition. Thus, we could easily determine the type of the phase transition. The proposed algorithm is general. We can apply this algorithm to various problems of any dimension. For example, the 2D ferromagnetic q-state Potts model with r invisible (redundant) states exhibits a change in the phase transition from the second order to the first order owing to the entropy effect of invisible states [25]. A study on the two-size PCC algorithm is now in progress for such a transition change problem.