Event-chain Monte Carlo algorithms for hard spheres 

Get Complete Project Material File(s) Now! »

Non-universality of the transition

In higher spatial dimensions, renormalization theory [2] suggests that the nature of the transition in spins systems does not depend upon the microscopic model. Rather, it depends solely upon the dimension and the degree of symmetry of the system. This property is called strong universality. In 2D systems, the strong universality is broken, as can be seen in the different XY models.
As discussed in Section 1.2, it was not expected that the transition in the 2D rotor XY model was necessarily of the KT type, the transition could have been first order [51]. A phase transition is said to be first order (see Section 1.3.1) when the system changes from one phase to the another discontinuously. This transition has the particularity that the phase in the disordered region has a finite correlation length at the transition point. In a continuous transition (as the KT transition or a second-order transition), the correlation length becomes infinite at the transition. A first-order transition is therefore referred to as “weak” when the correlation lengths (in unit of lattice spacing) are large, and “strong” when they are comparable to the atomic dimensions.
In order to understand the validity of the KT theory, we examine its assumptions. The KT theory predicts a transition mechanism between a phase whose correlation function is algebraic and a disordered phase. The existence of the algebraic phase is not debated as the harmonic approximation becomes exact for any microscopic model in the limit of low temperatures. Under the condition that the harmonic approximation is still valid up to the vortices unbinding, the renormalization treatment of the transition is correct and does not introduce any other assumptions. The KT theory is thus built on the assumption that the low-temperature phase remains stable up to the vortex-unbinding transition. While the temperature TKT yields an upper limit for the stability of the algebraic low-temperature phase, a transition at a lower temperature can preempt the KT transition.

Scenarios for two-dimensional melting

We now examine the case of 2D particle systems. The melting transition in these systems is more complex that for lattice systems, and it is the subject of a long-standing debate [17, 18]. In this section, the principal theories are presented. One is a classic first-order transition, and the other one follows the ideas of Kosterlitz and Thouless. This section exposes elements needed to understand the results of the hard-disk model in Chapter 4.

First-order transition

A transition is said to be first order when the system changes from one phase to another discontinuously. In 3D and higher, solids generally melt through a first-order transition. In two dimensions, the question is open. In hard disks, the loop seen in the equation of state of Alder and Wainwright was first interpreted as a classic first-order transition. In disk, between a liquid and a “hexatic” phase (see Section 4.1).
A first-order transition has the particularity to show a phase coexistence given an appropriate thermodynamic ensemble. For example, a first-order transition in a XY model shows a phase separation in the microcanonical (N, E) ensemble. This is a result of the energy per particle of both phases being different. For 3D hard spheres [49], the system is phase-separated at the transition in the constant-density ensemble (N,V, T).
Indeed, the densities of both phases are different. This phenomenon is understood by the concavity of the free energy Fpure(V) that a scenario without phase separation would have (see Fig. 1.9).

The KTHNY theory

The KT theory of Kosterlitz and Thouless described for the XY model in Section 1.2.2 has also been originally developed for the melting of 2D solids [13]. In this case, the solidliquid transition is ruled by unbinding of the topological defect of the positional order: the dislocations. This theory was completed by Nelson and Halperin [32, 68], as well as Young [69]. They noticed that the unbinding of dislocations did not lead to the disordered liquid phase but to a new phase, called hexatic. In the hexatic phase the positional order is short ranged and the orientational order is quasi-long ranged. The transition from the hexatic to the liquid is again a KT transition, ruled by the unbinding of the topological defects of the orientation: the disclinations. This two-stage melting theory is referred to as the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory. Later, various first-order scenarios, based on unbinding of grain boundary, vacancies condensation, or simultaneous unbinding of dislocation and disclination, were considered [70, 71, 72, 73].
In Chapter 4, we show that the melting transition for hard disks follows a two-stage scenario with a hexatic phase. The solid-hexatic transition is of the KT type, and the liquid-hexatic transition is first order. The following reviews the principal results of the KTHNY theory, which are similar to the results of the KT transition developed in the XY model.

Hard-sphere simulations

Molecular-dynamics simulation is the direct integration of Newton’s laws of motion, which describe the time evolution of particle systems. Usually, the equations are solved by discretizing the time in time step δt. For each time step, interactions between particles are computed and velocities updated. This discretization induces an error that is negligible if δt is small enough.
The hard-disk system is composed of elastic disks of mass m, and the state of the system is fully determined by the disks positions and velocities (see Fig. 2.0). In this system, molecular dynamics can be achieved without discretization. The time of the next collision (the next “event”) is computed, and the velocities after a collision are updated by momentum and energy conservation: ~v′1 +~v′2 = ~v1 +~v2 and v′2 1 + v′2 2 = v21 + v22 .

The event-chain algorithm

The following introduces a Markov-chain Monte Carlo algorithm developed in collaboration with Werner Krauth and David B. Wilson: the event-chain algorithm [5] (see Section 7.1). The event-chain algorithm is a rejection-free Markov-chain Monte Carlo algorithm designed for hard-core interactions. In a single move, it displaces an arbitrary long chain of particle in a long-range coherent motion. This algorithm has the particularity that it can be set up “irreversible”, (that is, without the detailed-balance condition).
This property, which is usually seen in simple Markov chains [104], increases the efficiency of the thermalization. This algorithm is the first algorithm that clearly outperforms the local Metropolis algorithm.


Straight event-chain algorithm

The simplest version of the event-chain algorithm is the “straight event-chain” (SEC) algorithm. A move of this algorithm starts by the sampling of an initial disk k and an angle θ. The disk k is displaced straight in the direction θ until it collides with another disk k′. The latter is then displaced in the same direction until it collides with yet another disk and so on (see Fig. 2.4). The move stops when all displacements add up to a total length parameter ℓ.

Stability and verifications

For long runs (up to 1014 displacements), rare events, such as rounding errors due to finite precision, can occur and disrupt the simulations. In order to prevent such an event, a simple addition was made to the algorithm: the collision length was taken as max(Lcoll , 0). In that way, if the collision length is Lcoll < 0 (meaning that the moving disk and collided disk are overlapping), the (backward) displacement is not accepted. The result is that the incoming disk does not move, and the collided overlapping disk is displaced forward: this reduces the overlapped region. As a global consequence, a starting configuration with many overlapping disks quickly “heals” by itself and becomes a valid configuration: the algorithm is stable and runs with more than 1014 displacements can be performed without trouble.
The implementation was carefully tested by comparing the distribution of a given observable obtained with the event-chain algorithm and with an independent algorithm such as the local Metropolis algorithm. To that purpose the orientational order parameter Y of Eq. (2.24) is used. In order to increase the efficiency of this test, the simulations are performed in the transition region (η ∼ 0.7), where a small violation of the balance condition has large consequences on the stationary distribution. The distribution for all the algorithms were found to be identical, confirming the validity of their stationary distribution (see Fig. 2.12).

Low-density phase: liquid

The nature of the low-density phase is determined by analyzing the pure phase at η ≃ 0.700. The system is ordered in clusters of maximum size ∼ 100σ −200σ and disordered above this scale; the phase is therefore liquid (see Fig. 4.9). The correlation length is more precisely computed through the orientational correlation function which behaves as Co(r) µ exp−x/ξo at large distances. This gives a correlation length of ξo ∼ 100σ (see Fig. 4.9 again). This large value partly explains the difficulty for small systems to discriminate between a first-order and a continuous transition. The correlation length is one order of magnitude larger than the disks, yet much shorter than the size of the system for 5122 and 10242 disks. This point is crucial for the validity of the results, and explains the collapse of P(v) for 5122 and 10242 disks in the pure-liquid branch (see Fig. 4.4).
Therefore, the behavior of an infinite system is obtained with these system sizes. The behavior of the correlation length at lower densities shows an important increase for η < 0.700 (see Fig. 4.10). This behavior can be wrongly interpreted as evidence for a KT transition. This behavior is consistent with the view that a KT transition is preempted by a first-order transition, and that for other microscopic models the transition could be continuous, as is seen in XY models (see Section 1.2.4). For “small” systems, the decay of the correlation function in the coexistence region can be wrongly interpreted as algebraic (see Fig. 4.10). It is the result of a two-phase region where the interface shows large fluctuations. At η = 0.708, it necessary to reach 10242 disk in order to see the two-phase behavior (see Fig. 4.10 again).

High-density phase: hexatic

The nature of the high-density phase is determined by analyzing the pure phase at η ∼ 0.718, which is clearly out of the transition region (see Fig. 4.4). The orientation at this density is conserved along the whole system (see Fig. 4.11 and Fig. 4.12). This does not indicate that the phase is solid as a cut of g(~r) along an axis shows that the positional order at this density is short ranged (see Fig. 4.11 again). The short-range positional order imposes that the shear modulus is zero [68], this phase can therefore not be a solid.
Moreover, in the absence of a shear modulus, the next term in the small-deformation development of the Hamiltonian is [32, 68] Hhex = 1 2 KA Z |~∇ θ(~r)|2 d2r , (4.13) KA being the Frank constant (see Chapter 1, Section 1.3.2). This shows that a phase with a short-range positional order cannot possesses a long-range orientational order. The orientational order is therefore quasi-long ranged, and the phase at η = 0.718 is hexatic.
The high-density phase in the coexistence region (η ∼ 0.716) shows the same property as for η = 0.718, it is therefore also hexatic. The analysis of η = 0.718 instead of η = 0.716 is motivated by the fact that η = 0.718 is clearly above the coexistence region.

Table of contents :

General introduction
I Two-dimensional melting 
1 Transitions in two dimensions 
1.1 Two-dimensional solids
1.2 Transitions in 2D XY models
1.3 Scenarios for two-dimensional melting
2 Event-chain Monte Carlo algorithms for hard spheres 
2.1 Hard-sphere simulations
2.2 The event-chain algorithm
2.3 Implementation of the algorithm
3 Simulation method 
3.1 Choice of ensemble
3.2 Statistical errors
3.3 Computation of observables
4 Hard-disk melting transition 
4.1 Phase coexistence
4.2 Nature of the phases
4.3 Thermalization and finite-size effects
4.4 Other simulations
II Perfect sampling 
5 The coupling-from-the-past method 
5.1 The Monte Carlo method
5.2 Coupling from the past
5.3 Survey problem
6 Damage spreading and coupling in Markov chains 
6.1 Chaos in Markov chains
6.2 Spin glass
6.3 Hard spheres
6.4 Griffeath’s coupling
7 Publications 
7.1 Publication I
7.2 Publication II
7.3 Publication III


Related Posts