Monte Carlo simulation of classical spin models with chaotic billiards

It has recently been shown that the computing abilities of Boltzmann machines, or Ising spin-glass models, can be implemented by chaotic billiard dynamics without any use of random numbers. In this paper, we further numerically investigate the capabilities of the chaotic billiard dynamics as a deterministic alternative to random Monte Carlo methods by applying it to classical spin models in statistical physics. First, we verify that the billiard dynamics can yield samples that converge to the true distribution of the Ising model on a small lattice, and we show that it appears to have the same convergence rate as random Monte Carlo sampling. Second, we apply the billiard dynamics to finite-size scaling analysis of the critical behavior of the Ising model and show that the phase transition point and the critical exponents are correctly obtained. Third, we extend the billiard dynamics to spins that take more than two states and show that it can be applied successfully to the Potts model. We also discuss the possibility of extensions to continuous-valued models such as the XY model.


I. INTRODUCTION
Many important classical spin models such as the Ising model and the Potts model are described as probability distributions of spin configurations. To investigate the behavior of such models numerically, we normally design Monte Carlo methods, and obtain samples from the model using random numbers. However, there is no reason in principle to use randomness for obtaining samples; it does not matter whether the samples are generated randomly or deterministically, if the samples properly represent the probabilistic models in question.
Chaotic Boltzmann machines recently proposed in Ref. [1] have chaotic billiard dynamics that yields samples from Ising spin-glass models without any use of random numbers. They have been numerically shown to have computing abilities comparable to conventional (stochastic) Boltzmann machines. In this paper, we further numerically investigate the capabilities of the chaotic billiard dynamics as a deterministic alternative to random Monte Carlo methods for classical spin models. Although there have been no studies that utilize billiard dynamics for deterministic simulations of spin models, the following three streams of studies can be considered as closely related to the present study.
Firstly, several deterministic cellular automaton models for the Ising model have been proposed. The Q2R automaton [2] is the simplest automaton model of the Ising model, and it evolves conserving the energy exactly. Creutz [3] proposed another automaton model in which demons are introduced as additional degrees of freedom. Demons absorb and release energy at each site conserving the total energy. Despite the simple and deterministic update rules, these models reproduce the probabilistic behavior of the Ising model. However, these automata cannot be directly regarded as deterministic samplers from the Ising model; in the Q2R automaton, spin configurations only move on a microcanonical ensemble, and in the Creutz automaton, the system temperature is internally * hideyuki@mist.i.u-tokyo.ac.jp determined. Although chaotic Boltzmann machines are not automata, there is a similarity to the Creutz automaton in the sense that the additional degrees of freedom introduced allows the system to deterministically generate samples from canonical ensembles.
Secondly, spin models based on coupled map lattices (CMLs) [4] have been proposed [5][6][7][8][9][10][11]. CMLs are typically composed of discrete-time chaotic elements on a lattice interacting with each other. They are known to exhibit rich spatiotemporal nonlinear dynamics, and form an important class of dynamical systems with a large number of degrees of freedom. By associating symbols to partitions in the state space of each element, CMLs can be regarded as deterministic Isinglike spin models. In the sense that the probabilistic behavior of Ising spins is realized by chaotic dynamics, CML-based spin models can be related to chaotic Boltzmann machines. However, each element constituting chaotic Boltzmann machines has continuous-time and non-chaotic dynamics.
Thirdly, it should be noted that random numbers used in conventional Monte Carlo simulations on ordinary computers, i.e., pseudo-random numbers, are deterministically generated. Therefore, in a sense, we have already been using deterministic Monte Carlo methods. Of course, such pseudo-random numbers are designed so that they can be regarded as truly random in many aspects. Since Monte Carlo methods rely on truly random numbers, it is crucial to use good random numbers in principle. However, the generation of good random numbers is costly, and this is one important issue in large-scale Monte Carlo simulations. Moreover, it has been pointed out [12] that low-quality pseudo-random numbers can actually be used in Monte Carlo methods and may even improve the performance. Thus, it is controversial as to the quality of randomness actually required. It is also noteworthy that recently a deterministic sampling algorithm, which is called herded Gibbs sampling [13], was proposed on the basis of the herding algorithm [14,15]. Although the algorithm is not practical for large-scale simulations, it is theoretically far more efficient than conventional Monte Carlo methods. Thus, it is intriguing to explore what is possible with deterministic Monte Carlo algorithms, and chaotic Boltzmann machines can be regarded as one of the approaches to the problem.
With these motivations in mind, in this paper, we numerically investigate the capabilities of the chaotic billiard dynamics as a deterministic alternative to random Monte Carlo methods by applying it to classical spin models in statistical physics. We confirm that the billiard dynamics can yield samples that converge to the true distribution of the Ising model, and we show that the phase transition point and the critical exponents of the Ising model are correctly obtained. Furthermore, we extend the billiard dynamics to spins that take more than two states, and we apply it to the Potts model and the XY model. We also point out that the billiard sampling dynamics is reversible and can be a good example for discussing the microscopic origins of irreversible macroscopic behavior on the Ising model.
It is to be noted that although the term "Monte Carlo" may imply that the algorithm is probabilistic, we intentionally keep using this term for deterministic algorithms also, thereby indicating that they can be used exactly in the place of conventional Monte Carlo methods.

II. BILLIARD DYNAMICS FOR THE ISING MODEL
In this section, we briefly introduce the implementation of probabilistic spin models by billiard dynamics. Although we limit the description to the Ising model, it can be extended to arbitrary Ising spin-glass models or Boltzmann machines [1].
Let us consider the Ising model composed of N sites on a lattice. The Hamiltonian for a spin configuration σ = (σ 1 , . . . , σ N ) ∈ {−1, +1} N of the Ising model is given by where the summation is taken for all the adjacent pairs in a lattice. The probability distribution of the spin configurations of the Ising model is given by the Gibbs distribution where T denotes the temperature and Z denotes the partition function given by For large spin systems, it is difficult to evaluate the probability and directly obtain the samples. To obtain samples from the probability distribution (2), we normally use Monte Carlo methods. Here we consider the heat-bath algorithm, which is also known as Gibbs sampling in the field of machine learning. We choose a spin i from 1, . . . , N randomly or sequentially. For each chosen spin, we update the state according to the probability where σ \i denotes the configuration of the spins in the system except for the ith spin. This process defines a Markov chain having the Gibbs distribution (2) as the stationary distribution. Therefore, we can eventually obtain a sample sequence from the Gibbs distribution.
Here, we consider using billiard dynamics instead of the heat-bath algorithm for sampling from the Gibbs distribution (2). We introduce an internal state x i ∈ [−1, +1] for each node i. The internal state of the ith node evolves according to the following differential equation: The state of the ith node σ i changes when x i reaches +1 or −1 as follows: Note that x i decreases when σ i = +1 and increases when σ i = −1. Therefore, the internal state x i continues oscillating between +1 and −1. The continuous-time dynamics defined by Eqs. (5) and (6) is a hybrid dynamical system [16] on the state space {−1, +1} N × [−1, +1] N , because it has both discrete and continuous state variables. The differential equation (5) is designed to be consistent with Eq. (4) in the following sense. We assume here that the states of the neighboring nodes of the ith node are fixed. Then, x i oscillates between +1 and −1 periodically. According to Eq. (5), it takes 2 exp( σ j /T ) unit time for x i to move from +1 to −1, and it takes 2 exp(− σ j /T ) unit time for x i to move from −1 to +1. Hence, the period is 2(exp( σ j /T ) + exp(− σ j /T )) unit time, in which the state σ i takes on +1 for 2 exp( σ j /T ) unit time. Therefore, the probability that we observe σ i = +1 at a random instant is consistent with Eq. (4). Note that this consistency is derived under the assumption that the states of the neighboring nodes are fixed. Since the states in the system actually change, this is merely an explanation that justifies Eq. (5) only intuitively. It is expected but not assured that σ follows the Gibbs distribution (2).
In the explanation above, it is only essential that the relative time duration for which σ i takes on +1 in a period coincides with the probability P[σ i = +1 | σ \i ] defined in Eq. (4). In other words, what is required for the consistency is that the speed |dx i /dt| is proportional to P[σ i | σ \i ] −1 . In fact, there are two other natural implementations in place of Eq. (5). One possibility is to define the speed as P[σ i | σ \i ] −1 as follows: Another possibility is to define the speed as P[−σ i | σ \i ] as follows: Considering the interactions with the neighboring nodes, these three definitions have different dynamics. As shown in the next section ( Fig. 2(a)), all the definitions appear to work similarly, so that we mainly employ only Eq. (5) in this paper. We also note that implementations are not limited only to these three types. For instance, arbitrary constants different from site to site can be multiplied to the speeds.
There are two views on the dynamics of this system. If we focus on each node in the system, the internal state x i continues oscillating between −1 and +1, interacting with the neighboring nodes. In this sense, the system can be regarded as a type of coupled oscillator system. On the other hand, if we observe all the internal states simultaneously, the internal state (x 1 , . . . , x N ) travels in a straight line in the hypercube [−1, +1] N according to Eq. (5), and changes its direction when it collides with the boundary of the hypercube according to Eq. (6). In this sense, the dynamics of this system can be regarded as pseudo-billiard [17]. Unlike deterministic spin models implemented by CMLs, each node does not have chaotic dynamics. However, it has been shown [1,18] that chaotic behavior naturally emerges from the interactions in the system, and it is considered to work as a heat-bath to realize probabilistic behavior of the spin configurations Statistics of the Ising spin model can be obtained from the continuous-time billiard system in the following manner. Let t 0 , t 1 , . . . , t k , . . . be the sequence of times at which state switchings occur according to Eq. (6). Then, the expectation value of a statistic Φ(σ) can be estimated from the samples until time t K as follows: (9) where τ k = t k − t k−1 and σ k denotes the state that the system is taking on in the time interval from t k−1 to t k . Thus, in a sense, we obtain a sample sequence σ 1 , σ 2 , . . . weighted by time intervals τ 1 , τ 2 , . . . . The statistics are calculated in this manner in the present study. Another method to obtain a (unweighted) sample sequence is to observe σ = (σ 1 , . . . , σ N ) uniformly, at every one unit time for example, or randomly. Note that sampling from a continuous-time billiard system is analogous to sampling with continuous-time (kinetic) Monte Carlo algorithms such as the Bortz-Kalos-Lebowitz (BKL) algorithm [19] and the Gillespie algorithm [20].
Although the billiard dynamics is described as a continuous-time system, numerical calculation can be carried out by iterating the Poincaré map x(t k ) → x(t k+1 ) induced on the boundary of the hypercube [0, 1] N as explained in Refs. [1,18]. Hence, when the system temperature is low and the spins seldom flip, the simulation performs efficiently, in a manner similar to continuous-time Monte Carlo algorithms. However, it is generally less efficient on ordinary computers compared with the ordinary heat-bath algorithm [1].

III. CONVERGENCE
In this section, we numerically verify that the billiard dynamics yields samples that converge to the true distribution of the Ising model on a small lattice, for which probability distributions can be computed. Figure 1 shows the absolute error between the empirical distribution observed from the billiard system and the true distribution of the Ising model on a two-dimensional lattice of size L = 4. The absolute error is defined as where r σ (t) denotes the total time for which the system takes on the state σ until time t. The solid line in Fig. 1 indicates the absolute error averaged for 96 different initial values, and the dashed lines indicate the minimum and maximum errors. Note that the error is calculated for each orbit. Therefore, Fig. 1 shows that for all the 96 initial values randomly chosen from the uniform distribution on    To summarize, all the results show convergence to the true distribution with the convergence rate almost in the order of O(1/ √ t). Namely, the convergence rate is almost the same as that of random Monte Carlo sampling, and slower than the order O(1/t) of the herded Gibbs sampling [13].

IV. FINITE-SIZE SCALING ANALYSIS
To examine the capabilities of the billiard dynamics for larger lattices, we apply it to finite-size scaling analysis (see, e.g., Ref. [21]) of the critical behavior of the Ising model. In the context of CMLs, it is known that the universality class depends on the updating rules [7]. Therefore, it is intriguing to determine as to which universality class the continuous-time billiard dynamics belongs to, although it should belong to the Ising universality class if the billiard dynamics truly yields samples from the Ising model precisely.
We use the Ising model on a two-dimensional lattice of size L = 32, 40, . . . , 80 with the periodic boundary condition. To find the phase transition point, we calculated the Binder cumulant [22]    The exponent ν can be obtained from the derivative of the Binder cumulant with respect to the temperature at the critical temperature, because according to the finite-size scaling theory, the following relation holds: Figure 4(a) shows the derivatives of the polynomial fits at the theoretical critical temperature. The gradient of the log-log plot is estimated as 0.993 ± 0.016, which is consistent with ν = 1 of the Ising model. The exponents β and γ can be obtained similarly from the absolute magnetization |m| and magnetic susceptibility χ = N( m 2 − |m| 2 ) at the critical temperature, using the following relation: From the log-log plots shown in Fig. 4(b) and (c), the critical exponents β/ν and γ/ν are estimated as 0.1248 ± 0.0004 and 1.751 ± 0.004, respectively. These values are consistent with those of the Ising universality class: ν = 1, β = 1/8, and γ = 7/4.

V. POTTS MODEL
In this section, we extend the billiard dynamics to spins that take more than two states, and we apply it to the Potts model.
For simplicity, we describe the dynamics for the Potts model here; it is straightforward to extend the dynamics to more general spin systems.
The Hamiltonian of the q-state Potts model is given by where σ i ∈ {0, 1, . . . , q − 1} denotes the state of the ith site and δ(·, ·) represents the Kronecker delta function.
We implement the Potts model on the basis of the billiard sampling dynamics for the Ising model as follows. We define an internal state x i ∈ R/qZ = [0, q) for each ith node, which evolves according to the following differential equation: Note that the speed dx i /dt is always positive and the two end points of [0, q) are regarded as identical. Therefore, x i continues to increase in the interval [0, q), and when it reaches x i = q, it instantaneously jumps to x i = 0. The state σ i of the ith node is determined from x i as σ i = ⌊x i ⌋, where ⌊x i ⌋ denotes the largest integer not greater than x i . Equation (16) is designed according to essentially the same idea as Eq. (5) of the billiard dynamics of the Ising model. Namely, the speed is determined so that it is proportional to P[σ i | σ \i ] −1 . Therefore, when the probability is higher, the internal state moves more slowly, and the system remains in such a state for a longer duration.
This system can be regarded as a coupled oscillator system, because each internal state x i oscillates on a circle R/qZ = [0, q) interacting with the neighboring nodes. When q = 2, this system is essentially equivalent to the billiard dynamics for the Ising model. However, for q > 2, the internal state (x 1 , . . . , x N ) travels on a N-dimensional torus [0, q) N , and its dynamics cannot be reduced to a billiard system. Note that, as mentioned in Ref. [1], another implementation by pseudobilliard dynamics for the Potts model is possible by using switched arrival systems [23]. As a statistic that characterizes the behavior of the Potts model, we define the order parameter as follows: where N s denotes the number of the sites taking on the state s. The order parameter m takes 1 when all the spins take one state (max s N s = N) in the completely ordered phase, and it takes 0 when the spins take all the states equally (N s = N/q) in the completely disordered phase. Figure 5 shows statistics of the Potts model on a twodimensional lattice of size L = 24 with the periodic boundary condition for q = 4 and q = 6 calculated by the heat-bath algorithm and the oscillator sampling dynamics. The Potts model on a two-dimensional lattice is known to undergo a second order phase transition if q ≤ 4 and a first order phase transition if q > 4. In Fig. 5, the average order parameter m and the average energy per site ǫ = E /N as well as their variances per site, the magnetic susceptibility χ = N( m 2 − m 2 ) and the specific heat c = N( ǫ 2 − ǫ 2 ), are shown. In all the graphs, the lines for the two methods coincide with each other.

A. Extensions to continuous-valued spin models
The implementation for the Potts model naturally leads us to consider extensions to continuous-valued models such as the XY model. Here, we discuss the possibilities of such extensions.
The XY model is composed of spins that take continuousvalued states in [0, 2π). The Hamiltonian of the XY model is given by where θ i ∈ [0, 2π) denotes the state of the ith spin. The XY model can be regarded as the limit q → ∞ of the Potts model with the interaction term replaced by the cosine function. One way to design sampling dynamics for the XY model on the basis of the oscillator dynamics for the Potts model is as follows. We define an internal state x i ∈ R/Z = [0, 1) for each ith node. It evolves according to the following differential equation: Since dx i /dt is always positive, the internal state x i always increases from 0 to 1. When it reaches x i = 1, it instantaneously jumps to x i = 0 and the state θ i is updated as follows: where φ is an arbitrary irrational number. We use the golden mean φ = ( √ 5 − 1)/2 in the numerical simulations. As another possibility, it is natural to consider the limit q → ∞ in the implementation for the Potts model (with the interaction term replaced by cos). Then, the internal state x i can be viewed as the continuous-valued spin state θ i of the XY model. More specifically, the state θ i evolves according to the following differential equation: This equation defines a coupled oscillator system, although the interaction is completely different from ordinary coupled oscillator systems such as the Kuramoto model. Since the state changes continuously, the numerical simulation of this system is very different from other systems we have considered. In the following simulations, we integrate the differential equation by the fourth-order Runge-Kutta method with a time step of 0.01.
We evaluate these methods by examining energy distributions constructed from sample sequences for the XY model at T = 1.2 and 0.6, where energy values per site are discretized with a bin width of 0.01. Since we do not have the true distribution, we regard an empirical distribution constructed from samples until t = 10 8 generated by the Metropolis algorithm as the "true" distribution, where N Metropolis steps are regarded as 1 unit time. Figure 6(a), (b), and (c) shows the absolute errors from the true distribution of the empirical distributions obtained by the Metropolis algorithm, the irrationalrotation sampling dynamics, and the coupled-oscillator sampling dynamics, respectively. For all the algorithms, sample sequences are obtained by sampling uniformly at every 1 unit time. As regards the irrational-rotation sampling ( Fig. 6(b)), the error decreases almost in the order of O( √ t) even in the worst case. While the coupled-oscillator sampling dynamics works well for T = 1.2, it does not work at all for T = 0.6 ( Fig. 6(c)). We have not succeeded in tracking down the cause of this result yet; it is possible that Eq. (21) intrinsically does not offer proper sampling dynamics even if the equation can be integrated without numerical errors. Although it is also possible that this may be only due to numerical errors, this result at least shows that naive numerical integration does not work well.
To summarize, we have examined two types of sampling dynamics designed for continuous-valued spin models. While the irrational-rotation sampling dynamics (Eq. (19)) updates the state θ i discretely, the coupled-oscillator sampling dynamics (Eq. (21)) updates the state θ i continuously. Although the irrational-rotation sampling seems promising, further studies, especially on the differences between these two methods, are necessary to understand the capabilities and the limitations of the proposed dynamics.

B. Spin echoes in the Ising model
Thus far, we have mainly described the billiard and oscillator sampling dynamics as a sampling method for classical spin systems. However, we note that the dynamics itself is also interesting as an abstract model for physical systems with a large number of degrees of freedom that can be related to both the CMLs and the coupled oscillator systems.
One important aspect of the billiard sampling dynamics is its reversibility, while its macroscopic behavior as spin models is irreversible. The microscopic origins of irreversible macroscopic behavior have been an important topic in statistical physics since more than a century ago (see, e.g., Ref. [24]). This topic has been discussed using models such as the Lorentz gas model (or Sinai's billiard) [25], multibaker maps [26,27], the Q2R automaton [28], and Nosé-Hoover thermostats [29][30][31]. The billiard sampling dynamics also ex- hibits both reversible microscopic behavior and irreversible macroscopic behavior, and therefore, it can be a good example for discussing this topic, particularly if its chaotic dynamics is investigated in a more theoretical manner in the future. Note that the bakermap lattice [9] for the Ising model can be a good example as well, because its dynamics is invertible and can be modified to be reversible.
As an example that demonstrates the reversibility, we show in Fig. 7 that spin echoes [32] can be observed in the Ising model with the billiard sampling dynamics. We use the Ising model in the paramagnetic (disordered) phase. At time t = 0, all the spins are aligned to +1, so that the magnetization is initially equal to 1. The initial internal state x is drawn from the uniform distribution on [−1, +1] N . As the system equilibrates, the magnetization decreases to values around zero. At time t = 50, we flip the signs of the internal states x instantaneously as x i → −x i , while the states σ are kept unchanged. After the flip, the simulation continues as if time is reversed. The magnetization remains in the equilibrium state for a while; however, it suddenly increases to 1 at time t = 100. Namely, all the spins are aligned again at the moment, contradictory to the second law of thermodynamics. This result is a natural consequence of the reversible dynamics, and theoretically, the memory of the initial state can be restored after an arbitrarily long time. However, it will practically become difficult to restore the initial state after a long equilibration time on a large lattice, due to the sensitive dependence on the initial conditions of the chaotic dynamics. This indicates that it will become almost impossible to find a microscopic state that evolves contrary to the second law of thermodynamics.

VII. SUMMARY
In this paper, we have numerically verified that the billiard dynamics can generate samples from the Ising model by examining the convergence and applying it to finite-size scaling analysis. We also have extended the billiard dynamics to multi-valued and continuous-valued spin models. In all the simulations (except for the oscillator dynamics for the XY model), the proposed dynamics works well as a deterministic alternative to random Monte Carlo sampling.
It is considered important to examine more spin models under various conditions on various network structures numerically. However, because there are infinitely many possible models to examine, we cannot clarify the capabilities and the limitations of this approach only by the means of numerical simulations. Therefore, we consider future studies on the theoretical foundations of the sampling dynamics to be also important on the basis of the promising results presented in this paper.