We consider a Markov chain with samples x t at time t = 0,1 , …   in a sample space Ω . In our case, its transition matrix P implements the heat-bath algorithm [17, 26, 27] (in other words, Glauber dynamics) for the Edwards–Anderson model or a version of the Metropolis algorithm [43] for hard spheres. We define the element P ( x , x ′ ) as the conditional probability to move from configuration x at time t to configuration x ′ at time t + 1 . With an initial configuration x 0 , the distribution π { t = 0 } is a delta function centered at x 0 . The distribution evolves over time as π { t + 1 } ( x ′ ) = ∑ x π { t } ( x ) P ( x , x ′ ) for each time step t . The approach to equilibrium is quantified by the mixing time, which is the time it takes for π { t } (which depends on the choice of x 0 ) to approach the stationary distribution π { t → ∞ } = π : τ mix ϵ = min t max x 0 ∈ Ω π { t } − π TV ⏟ d t < ϵ . (1) Here, ‖ ⋯ ‖ TV denotes the total variation distance [39], that is, one half of the absolute difference between π { t } and π over all the sample space, and ϵ is an arbitrary positive parameter that must be taken smaller than 1 2 . In Equation 1, the “ max ” refers to the worst initial choice for the approach of π { t } (which depends on x 0 ) to π , and this allows one to define the distance d ( t ) between π { t } and π , without explicit reference to the starting distribution π { t = 0 } . The mixing time is a non-asymptotic time scale [2] that describes the initial approach of π { t } toward the equilibrium distribution π on a finite distance scale ϵ . It comes with an exponential bound, valid from τ mix up to t → ∞ , while the asymptotic approach toward equilibrium, described by the (absolute) inverse gap of the transition matrix, can be much faster [39].

For a given transition matrix P of a Markov chain on a sample space Ω , a coupling is defined as a bivariate stochastic process with a configuration ( x t , y t ) at time t on the sample space Ω × Ω , such that P x t + 1 = x ′ | x t , y t = x , y = P x , x ′ , P y t + 1 = y ′ | x t , y t = x , y = P y , y ′ .

The bivariate process that updates the two copies x and y need not be Markovian [28] at a difference of its two projections. Non-Markovian couplings are theoretically important but have not been used yet in applications. Markovian couplings are described by a transition matrix P coup [ ( ⋅ ) , ( ⋅ ) ] on the sample space Ω × Ω that satisfies ∑ y ′ P coup [ ( x , y ) , ( x ′ , y ′ ) ] = P x , x ′ , ∑ x ′ P coup [ ( x , y ) , ( x ′ , y ′ ) ] = P y , y ′ , so that the transition matrix of the coupled Markov chain, which acts on two copies of the sample space Ω , when projected on either copy, returns the original transition matrix.

Couplings can take a variety of forms. The “classic” coupling performs two statistically independent Markov chains until, by accident, they couple, from when on they are glued together: P coup [ ( x , y ) , ( x ′ , y ′ ) ] = P x , x ′ P y , y ′ if  x ≠ y , P x , x ′ if  x = y , x ′ = y ′ , 0 if  x = y , x ′ ≠ y ′ , (2) (see the lhs of Figure 1). At the coupling time τ coup , the trajectories first meet: τ coup = min t x t = y t .

Transition matrices, as the ones in Equation 2, are implemented in Monte Carlo algorithms with the use of random elements, that is, one or several random numbers for selecting a particle or a spin, for choosing a move, and for accepting or rejecting it, etc. For example, the move from x at time t may produce an outcome x ′ that depends on the realization of the random element, but when this element is specified, as ϒ t ( x ) , it becomes a function called a random map { t } × Ω → Ω : x → x ′ = ϕ x , ϒ t ( x ) . The random map ϕ x , ϒ t ( x ) implementing this move must satisfy P ϕ x , ϒ t x = x ′ = P x , x ′ , as it must reproduce the transition matrix P . A random map ϕ (also called a “grand” coupling [39]) specifies a coupling, and it automatically implements a “gluing” operation, as two Markov chains that meet at a position x at time t encounter the same random element. For the classic coupling of Equation 2, the randomness at time t is a vector ϒ t = { ϒ t ( x ) : x ∈ Ω } of i.i.d random variables, that is of random numbers drawn from the same distribution (see center of Figure 1). For the “random-share” coupling, one uses, at time t , the same random element for all configurations x ∈ Ω : ϒ t = { ϒ t , … , ϒ t } . In other words, all configurations are updated with the same random numbers. Many other couplings exist, and it is only important that the projection onto a single copy produces a valid Markov chain. While every random map corresponds to a coupling, it appears that not all couplings (for example, the path couplings in Ref. [32]) can be expressed as random maps.

The connection between mixing times and coupling times is as follows ([39], corollary 5.3): d t ≤ max x 0 , y 0 ∈ Ω P x 0 , y 0 τ coup > t , (3) where d ( t ) is the distance entering the definition of the mixing time in Equation 1. From our previous discussion, it is evident that for random walks on large graphs, the classic coupling time can be much larger than the mixing time simply because the two Markov chains must hit the same configuration at the same time. In contrast, the random-share coupling time is of the same order as the mixing time for many random walks. In the problems at the focus of this article, we will witness different regimes, as a function of external parameters, that are separated by a phase transition. In this context, it is of great interest that an optimal coupling [28] realizes the coupling at time t and at position x t of two Markov chains that have started at time t = 0 at configurations x 0 and y 0 with the minimum of the probabilities to go from x 0 or from y 0 to x t . The optimal coupling is non-Markovian and virtually impossible to construct in practice, but it demonstrates that the bound of Equation 3 can be saturated.

A special class of couplings for which the inequality of Equation 3 can be tight (up to logarithms) requires the concept of monotonicity. In monotone couplings, there exists a partial ordering “ ⪯ ” of configurations so that x t ⪯ y t implies x t + 1 ⪯ y t + 1 . In terms of the random map, x ⪯ y implies ϕ ( x , ϒ ) ⪯ ϕ ( y , ϒ ) . No partial ordering exists for the random walk on the path graph with a classic coupling, and trajectories of Markov chains may cross (see the lhs of Figure 1). In contrast, for the random-share coupling of the one-dimensional random walk (which is a grand coupling), the ordering is complete. For a monotone grand coupling, with l the length of the longest “chain” in the partially ordered subset, the mean coupling time τ coup satisfies τ coup < 2 τ mix 1 / e 1 + log l . (4)

With Equation 3, there are thus upper and lower bounds for the monotone coupling time in terms of the mixing time, and the two agree up to a logarithm. For a monotone coupling with extremal elements, one must only survey their evolution, which will bracket all other configurations (see the rhs of Figure 1). Full surveys are possible in other cases [15], but the upper bound in Equation 4 is then often lost.

We can consider families of Markov chains that correspond to physical systems with size N , which may represent the number of sites, spins, or particles. As N increases and approaches infinity, under suitable conditions, such as constant temperature for spin systems or constant density for particle systems, the behavior of these systems can be studied. We may refer to “fast” coupling if the mean coupling time τ coup scales not slower than a power of the system size N (in later sections, we will use an N log N scaling).

As mentioned in the introduction, we may imagine the worst-case initial configurations x 0 and y 0 as the end points of a path of configurations, with adjacent elements on the path being neighbors, with respect to some metric. Under some conditions, it is often possible to show that any pair of neighboring configurations come in expectation even closer after one step of the Markov chain, and this establishes that the distance between x 1 and y 1 contracts, and similarly for later times, leading to a proof of fast coupling [13].

The path-coupling analysis that is local in sample space and in time yet valid uniformly for any pair of neighboring configurations yields a rigorous global fast-coupling bound. We will discuss the limiting temperature T path for spin glasses and limiting density η path for hard-sphere systems for which the uniform contraction allows one to prove fast coupling. However, the path-coupling approach is quite conservative. Numerical evidence [4] indicates fast coupling down to a temperature T coup that is lower than T path , and up to a density η coup that is higher than η path . However, only T path and η path are known analytically. In the models that we study, the coupling is either exponential (and thus “slow”) or “fast.”

The path-coupling analysis provides a justification for “damage spreading,” which has been extensively studied for spin systems in the physics literature, with the random-share coupling. As in path coupling, two neighboring initial configurations x 0 and y 0 were chosen and were followed for very large times. The explicit relationship between the time to couple and the time to mix is lost, but the mean coupling time starting from neighboring initial configurations is again exponential below T coup and ∼ N log N or faster above. The connection between coupling and damage spreading was made in [4].

The coupling time in Equation 3 that allows bounding the mixing time follows the worst-case pair of starting configurations, x 0 and y 0 . For monotone coupling, these configurations are given by the two extremal elements, but in general, this requires a survey of the entire sample space. For the Glauber dynamics of spin glasses with the random-share coupling, the patch algorithm [15] rigorously surveys the | Ω | ∼ 1 0 600 configurations on a 64 × 64 lattice, and the same algorithm also applies to hard-sphere models, where it allows one to establish the grand-coupling time [4, 16]. It was found, however, that a few hundred random initial configurations contained worst-case pairs with high probability. Such a partial-survey approximation is easy to set up in practice.

We use the partial-survey approximation to evaluate the mean coupling time τ coup for spin-glass and hard-sphere systems. Here, a systematic numerical approach, inspired by the finite-size scaling analysis of second-order phase transitions, is discussed for distinguishing between fast and slow couplings. In this context, fixing the system size N corresponds to limiting the worst-case pair distance between initial configurations, and the scaling behavior is analyzed as N grows by varying N . Suppose we obtain τ coup ( N , β ) numerically as a function of the system size N and the model parameter β , which represents the inverse temperature in the case of spin-glass systems. For hard-sphere systems, this parameter may also be the density η . In the fast-coupling regime, the size dependence of τ coup exhibits N ⁡ log ⁡ N behavior at high temperatures, while in the slow-coupling regime, it increases exponentially at low temperatures. This phenomenon can be viewed as a dynamical phase transition, with the two behaviors changing at a certain critical temperature β coup .

Assuming that, as β approaches β coup , N * ( β ) provides a diverging scale that controls the coupling behavior, the scaling form is postulated to hold in the vicinity of β coup , expressed as τ coup N , β = N ϕ f N / N * β w i t h N * β = | β coup − β | − ω , (5) where ϕ and ω are positive parameters associated with the dynamical transition, and f is a universal scaling function. The two behaviors of fast and slow couplings are represented in the asymptotic form of this scaling function f ( x ) , with x = N | β coup − β | ω : f x = x 1 − ϕ log x as  x → ∞  for     β coup > β , exp a x as  x → ∞  for     β coup < β , with a positive constant a . The value of the scaling function f ( 0 ) at β = β coup is constant, and the parameter ϕ can be identified as the exponent of the power-law divergence of τ coup at β coup . In the case of a ferromagnetic Ising model with monotone coupling, where the coupling time and the mixing time coincide, these parameters characterize the universality class of the corresponding ferromagnetic phase transition and are related to the dynamical exponent z and the correlation length exponent ν through the dimensionality d . For example, in the case of the mean-field ferromagnetic Ising model, it has been rigorously shown that ϕ = 3 / 2 [19], which is consistent with z = 2 . However, in general, the singularity at β coup in this coupling time is not directly associated with an order parameter of the physical system.

From a mathematical perspective, the fact that the interactions { J i j } are quenched random variables complicates the analysis with respect to uniform interactions. The Sherrington–Kirkpatrick model [52], in other words, the Edwards–Anderson model on a complete graph corresponding to its infinite-dimensional limit, has been at the forefront of theoretical developments in spin-glass research. This model undergoes a thermodynamic phase transition separating a high-temperature paramagnetic phase from a low-temperature spin-glass phase at an exactly known temperature. The existence of this phase transition and the low-temperature properties were first established using the replica method [44] and later proven rigorously [54]. The study on the domain-wall free energy [51], which incorporates the fluctuation effects at the mean-field level, has indicated that the lower critical dimension is 2.5, which lies between the dimensions of 2 and 3.

Mathematically rigorous results for the Edwards–Anderson model in finite dimensions are very few. In systems with random interactions, local regions may exhibit low probabilities but strong correlations, leading to anomalous singularities in the free energy and divergences in high-temperature expansions. In a specific random system, the existence of this type of singularity has been mathematically proven and is known as the Griffiths singularity [29]. This singularity emerges at the phase transition temperature when the random interactions are assumed to be uniform. In the Edwards–Anderson model, the Curie temperature of the ferromagnetic Ising model (with all J i j equal to + J ) constitutes this Griffiths temperature. Despite these difficulties, it has been proven that, at sufficiently high temperature, the Edwards–Anderson order parameter vanishes identically, and the spin-glass susceptibility remains finite in short-range spin-glass models [6, 25]. This means that the high-temperature phase is paramagnetic, although rigorous temperature bounds seem to be absent. These temperature regions are far from the spin-glass transition temperature T SG suggested by the numerical simulations mentioned below. One expects that a spin-glass phase cannot exist at temperatures higher than the Griffiths temperature, so the Griffiths temperature likely serves as an upper bound for T SG . However, this seems not to be a rigorous statement.

Early numerical studies [12, 42] on domain-wall energies at zero temperature, though limited to small system sizes, were the first to propose the existence of a finite-temperature spin-glass transition in three dimensions and the absence of such a transition in two dimensions. These findings were subsequently strengthened by exact algorithms in two dimensions and more sophisticated heuristic algorithms [10], which allowed for larger system sizes and more accurate results. Following them, local Monte Carlo methods, particularly those using the heat-bath algorithm, played a crucial role in confirming these conclusions. These Monte Carlo studies provided direct evidence for a finite-temperature transition in three dimensions [7, 8, 46, 47] and the absence of such a transition in two dimensions [7, 34]. While neither has been proven rigorously, the fact that the ground state of the two-dimensional Edwards–Anderson model can be computed in a time polynomial in N [9, 55] is compatible with the hypothesis that complex phase transitions are unlikely to occur in systems where the ground state can be easily obtained. However, it should be noted that certain systems, such as the random-field Ising model [45], allow for efficient ground-state calculations yet still exhibit complex phase transitions at finite temperatures. These conclusions, both for three and higher dimensions as well as for two dimensions, were based on estimates of spin-glass order parameters. These order parameters examine the degree to which the equilibrium running averages of a given observable, such as the spin overlap between replicated systems, become independent of two independent starting configurations in Monte Carlo simulations ([7], Equation 4). Another route to studying spin glasses has consisted in analyzing the autocorrelation functions of observables (e.g., the value of σ k ( t ) ). Early results already pointed to a difference in the scaling behavior at late times ([46], Figure 7), [47], from which a finite spin-glass transition temperature in the range T SG ≃ 1.10 − 1.14 was inferred. Although no consensus has been reached on the nature of the spin-glass phase, more recent studies have refined estimates of the spin-glass transition temperature T SG in three dimensions, with different estimates such as T SG = 1.1019 ( 29 ) [3] and T SG = 1.109 ( 10 ) [30], which combine simulations for rather small system sizes with empirical extrapolations to the thermodynamic limit.

Damage spreading in spin-glass systems was found as a dynamical anomaly in early numerical simulations [14, 18], which showed that it occurs at temperatures higher than the spin-glass transition temperature suggested by other studies. However, it remained unclear whether the anomaly was related to the spin glass transition itself or to the Griffiths singularity. The connection between damage spreading and coupling, which is the focus of this article, was recognized in Ref. [4].

In the finite-dimensional Edwards–Anderson model, we now consider the random-share coupling for the heat-bath algorithm of Equation 7. To establish coupling, we consider two arbitrary spin configurations as initial states of the two Markov chains and apply the path-coupling argument of Section 2.2. The two configurations differ in at most N sites so that we can connect them by a path of at most N neighboring configurations that differ by one spin only.

Let σ A and σ B be two such neighboring configurations (see the lhs of Figure 2) that differ by the spin j . The common random element { k , ϒ } of Equation 6 contains the spin k to be updated and the random number ϒ required for the heat-bath step of Equation 7. With probability p 1 → 0 = 1 / N , the spin j is updated. The field h j is the same for σ A and σ B , and so is ϒ in Equation 7. It follows that the distance decreases from 1 to 0 with p 1 → 0 .

With probability 2 d / N , spin l , one of the 2 d neighboring spins of j , is updated. The local fields h l ( σ A ) and h l ( σ B ) differ by exactly 2. The probability p 1 → 2 of making different decisions, which corresponds at most to the red region on the rhs of Figure 2, is at most equal to p 1 → 2 = 2 d N ⁢ max h ⁢ π + ⁢ h − π + ⁢ h ± 2 = 2 d N ⁢ π + ⁢ 0 − π + ⁢ − 2 = 2 d N ⁢ 1 2 − 1 1 + exp ( 4 β ) .

If p 1 → 0 > p 1 → 2 , the expected distance between σ A and σ B decreases after one step, for any choice of spin configuration and any choice of the couplings { J i j } , which is the case at high temperature. This is also a condition where the damage caused by a single spin difference does not spread in the initial stage of the damage spreading under random-share coupling. It provides the upper bound of the damage spreading temperature. More details are discussed in Section 3.3. The limiting temperature for the application of the path-coupling argument is when p 1 → 0 = p 1 → 2 , which translates into β path = 1 4 log 2 d d − 1 − 1 = 1 2 d + 1 6 d 3 + 1 10 d 5 + ⋯   , and equivalently, T path = 1 β path = 2 d − 2 3 d − 8 45 d 3 + ⋯   . (8)

For T > T path , we are assured of fast coupling in the Edwards–Anderson model. The argument also holds for sublattice parallel updates. As discussed, T path is obtained for any choice of interactions and any spin configuration. Consequently, T path is also the path-coupling bound for the ferromagnetic Ising model, although we know from monotonicity that fast coupling will take place down to the Curie temperature.

We now numerically evaluate the mean coupling time of the finite-dimensional Edwards–Anderson model in both two and three dimensions in view of the scaling analysis discussed in Section 2.3. The mean coupling time of the two-dimensional model was already evaluated under a random update rule, and it has been demonstrated that a dynamical phase transition occurs in which the size dependence of the coupling time qualitatively changes [4], confirming earlier results [14]. The mean coupling time results presented below are evaluated using the partial-survey approximation with the number N 0 of randomly chosen initial conditions. The results obtained with different values of N 0 are plotted at each data point, but they are completely contained within the size of the markers, thereby confirming that they are independent of N 0 . A dendrogram representation explains the independence of the mean coupling time of N 0 (see Figure 3).

All the figures shown below represent results averaged over 4,096 realizations of interactions, independent of N , with error bars indicating sample fluctuations from these realizations. The first results for the three-dimensional Edwards–Anderson model are presented in the two panels of Figure 4, which show the estimated mean coupling time for the partial-survey approximation under the parallel and random updates. Although the two updates differ in the high-temperature limit, both exhibit a N ⁡ log ⁡ N behavior for system size N at sufficiently high but finite temperatures. As the temperature decreases, the behavior of the N dependence of the mean coupling time changes from slow to fast increase at a certain temperature. There is a slight, yet significant, difference in the transition temperature between the two updates, with a lower transition temperature observed for the parallel updates. This illustrates that coupling has no direct thermodynamic significance.

Figure 5 presents finite-size scaling plots of the mean coupling time for the three-dimensional Edwards–Anderson model, comparing both the parallel and random updates. The plot demonstrates that the scaling works well when the appropriate scaling parameters are chosen. This is consistent with the above argument that the transition temperatures, T coup or β coup , are significantly different for the two update rules. In contrast, the precision of the scaling exponents, ϕ and ω , is not as precise as that of the transition temperature, and it can be considered that these two rules yield almost the same values for these exponents. It remains unclear whether these exponents have a meaning analogous to the critical exponents of a second-order transition. Of particular interest is the exponent ω , which represents the divergence of the characteristic scale as it approaches the transition temperature. Our results suggest that this exponent has the same value on both the high- and low-temperature sides of the transition temperature. This is comparable to the correlation length exponent.

An analogous scaling analysis for the two-dimensional Edwards–Anderson model is shown in Figure 6. The left panel is the analysis result of our own numerical simulations using the sublattice parallel update, while the right panel presents the scaling analysis based on numerical data using the random update from [4]. In both cases, the scaling is consistent with a phase transition in the mean coupling time. As observed in the three-dimensional model, T coup depends on the underlying Markov chain, with a lower transition temperature for the parallel update. The scaling exponents depend on the dimensionality. However, the proper scaling variable may not be the number of spins, N , used here, but rather the linear dimension L . This suggests that the value of the exponents may depend on the dimensionality through the relationship N = L d .

Table 1 summarizes the key temperatures discussed in previous sections, including T path and T coup , as well as previously estimated results for T SG and T Griffiths . This table demonstrates the differences in transition temperatures for both two- and three-dimensional Edwards–Anderson models, providing a detailed overview of the coupling and spin-glass transitions.

On the one hand, path coupling demonstrates that above T path , the uniform contraction between neighboring configurations leads to fast coupling. Below T path , there are spins k (for example, those with h k = 0 ) for which, at least initially, there is no such contraction. Nevertheless, as our numerical simulations show, fast N log N coupling also takes place in the window T coup < T < T path . The absence of a regime change at T path can be illustrated, in the language of damage spreading, by following the mean damage as a function of time for two configurations that initially, at time t = 0 , are neighboring. Above T path , the mean damage decreases exponentially for all times (see inset of Figure 7), whereas for T < T coup , it increases rapidly. In the window T coup ≲ T ≲ T path , the mean damage initially increases, as expected, but then turns around and again vanishes exponentially. This turning point seems to occur when the damage reaches a certain size, which grows as the temperature approaches T coup . This behavior can be understood in analogy with the characteristic diverging scale N * ( β ) in the finite-size scaling analysis of Equation 5, which suggests a picture similar to a critical phase transition, where the threshold damage size corresponds to the diverging scale near T coup .

Rigorous results on hard-disk (and hard-sphere) models are very few. It is known that the close-packing density η = π / ( 2 3 ) in two dimensions is characterized by the hexagonal packing [24]. It thus corresponds to an essentially unique configuration that has long-range orientational and positional order. For densities below the close-packing density, the absence of long-range positional order was established rigorously [50] so that there is no crystal (with long-range orientational and positional order) below close packing. Indications for a phase transition were first found in the 1960s [1, 40]. The existence of two phase transitions and of three phases (liquid, hexatic, and solid) as a function of density is now well accepted [5, 40]. As in the Edwards–Anderson model (where the temperature T replaces the inverse of the density η as a control parameter), a rigorous proof of a transition away from close packing is still lacking. At low finite densities, the convergence of the virial expansion was proven early on [38], establishing the existence of the liquid phase. It extends up to a density η = 0.70 and is followed by a window of coexisting liquid and hexatic regions (see Table 2 below).

We now consider path coupling for hard disks, using the random map based on Equation 9 and a Hamming metric that counts the number of different disk positions in any two configurations. Let X A and X B be two neighboring hard-disk configurations that differ in the position of disk j only (see Figure 8). Simplifying a coupling from Ref. [35], we use as the common random element { k , ϒ } the disk k to be updated and its new position, both identical for X A and X B . With probability 1 / N , the disk j is moved (that is, k = j ). The move is accepted in both configurations if it stays away (by 2 σ ) from the “halo” of all remaining disks in both configurations. This yields the probability of decreasing the Hamming distance from 1 to 0: p 1 → 0 ≥ 1 N 1 − N − 1 N 4 η . (10)

On the other hand, the Hamming distance can be increased from 1 to 2 if a disk different from j is moved less than 2 σ away (that is, into the halo), of disk j in one configuration but not in the other. The probability of increasing the Hamming distance from one to two can thus be bounded as: p 1 → 2 ≤ N − 1 N 8 N η , (11) where the factor 8 η / N on the rhs arises from the difference between two “halos” of area 4 π σ 2 for each of the disks j in the two configurations.

Again, for p 1 → 0 > p 1 → 2 , the expected Hamming distance between A and B decreases after one step, for any two neighboring disk configurations, which can be assured for η < η path = 1 12 . (13)

It follows [13] that the Hamming distance between configurations A and B that differ in the position of only one disk decreases in expectation at each step if the density is smaller than 1 / 12 .

As with the Edwards–Anderson model, we now analyze the mean coupling time of the two-dimensional hard-disk model under the global Metropolis algorithm with the random-share coupling of Equation 9. In this case, we reanalyze the data obtained in Ref. [4], which we replot on the lhs of Figure 9. The analogous scaling ansatz again provides an excellent fit of the data. The critical exponents do not differ significantly from those found in the Edwards–Anderson model, suggesting the possibility of some underlying universality. However, uncovering the intricate physical picture behind this similarity remains an open question for future research. It should be noted that these critical exponents are not directly related to the critical phenomena of physical systems in the conventional sense. Rather, they characterize the “phase transition” in computational algorithms associated with the coupling of Markov chains. From an algorithmic perspective, these exponents are of significant interest as they provide insight into the inherent challenges in achieving fast coupling.

The coupling approach to the hard-disk system has been intensely studied in recent years, and the random-share coupling of Equation 9 only provides the simplest possible choice. A number of refined couplings have been proposed. The one proposed in Ref. [35] moves disks differently for the configuration X A and X B and reaches a path-coupling density of 1 / 8 (see Table 2 for an overview). Building on this coupling, optimizing the metric reaches a limiting density of 0.154, which was later improved for a different algorithm to 1 / 6 . In addition to these rigorous bounds, numerical evidence for the birth–death algorithm [56] points to a coupling density of ∼ 0.3 [4]. These densities, and especially the rigorously proven ones, are still quite far from the “empirical” transition density η ∼ 0.70 of the liquid phase, which was only in recent years understood to be toward a hexatic, and which bounds on a region η ∈ [ 0.7 , 0.76 ] without a homogeneous solution, and then giving rise to a mixture of the hexatic and the liquid.

The crucial connection between fast coupling (thus, fast mixing) and physical ordering was made for the hard-sphere case in Ref. [32], where it was proven that O ( N log N ) random steps of the global Metropolis algorithm are insufficient to construct configurations with any kind of long-range order. Fast mixing of a single-particle algorithm, even a non-local one, thus implies that the resulting configuration (which is practically in equilibrium) has exponential spatial correlation functions. This, to all intents and purposes, shows the extension of the liquid phase. We believe that it does not, however, prove the convergence of the virial expansion [38] because of the possibility of a liquid–liquid phase transition, which cannot be captured in a mixing-time argument.