Event-chain Monte Carlo: foundations, applications, and prospects

This review treats the mathematical and algorithmic foundations of non-reversible Markov chains in the context of event-chain Monte Carlo (ECMC), a continuous-time lifted Markov chain that employs the factorized Metropolis algorithm. It analyzes a number of model applications, and then reviews the formulation as well as the performance of ECMC in key models in statistical physics. Finally, the review reports on an ongoing initiative to apply the method to the sampling problem in molecular simulation, that is, to real-world models of peptides, proteins, and polymers in aqueous solution.


INTRODUCTION
Markov-chain Monte Carlo (MCMC) is an essential tool for the natural sciences. It is also the subject of a research discipline in mathematics. MCMC has from the beginning [67], in 1953, focussed on reversible MCMC algorithms, those that satisfy a detailed-balance condition. These algorithms are particularly powerful when the Monte Carlo moves (from one configuration to the next) can be customized for each configuration of a given system. Such a priori choices allow for big moves to be accepted, and for sample space to be explored rapidly. Prominent examples for reversible methods with custom-built, large-scale, moves are path-integral Monte Carlo and the cluster algorithms for spin systems [49].
In many important problems, insightful a priori choices for moves are yet unavailable. MCMC then often consists of a sequence of unbiased local moves, such as tiny displacements of one out of N particles, or flips of one spin out of many. Local reversible Monte Carlo schemes are easy to set up for these problems. The Metropolis or the heatbath (Gibbs-sampling) algorithms are popular choices. They generally compute acceptance probabilities from the changes in the total potential (the system's energy) and thus mimic the behavior of physical systems in thermodynamic equilibrium. However, such algorithms are often too slow to be useful. Examples are the hard-disk model [1,7] where local reversible MCMC methods failed for several decades to obtain independent samples in large systems, as well as the vast field of molecular simulation [25]. That field considers classical models of polymers and proteins, etc., in aqueous solution. Local reversible MCMC algorithms were for a long time without alternatives in molecular simulation, but they remain of little use because of their extreme inefficiency.
Event-chain Monte Carlo (ECMC) [8,70] is a family of local, non-reversible MCMC algorithms that has been developed over the last decade. Its foundations, applications, and prospects are at the heart of the present review. At a difference with its local reversible counterparts (that are all essentially equivalent to each other and to the physical overdamped equilibrium dynamics), different representatives of ECMC present a wide spread of behaviors for the same problem. Two issues explain this spread. First, any decision to accept or reject a given move is made on the basis of a consensus [70] between statistically independent factors, which each regroup parts of the total potential. Usually, there is a choice between inequivalent factorizations. In the Metropolis algorithm, in contrast, all decisions to accept a move are made on the basis of changes in the total energy. Second, ECMC is fundamentally a non-reversible, "lifted", version of an underlying reversible algorithm. In a lifted Markov chain, some of the randomly sampled moves of the original ("collapsed") Markov chain are rearranged. Again, there are many options for rearranging moves, and each of these can give rise to a specific dynamic behavior. The development of ECMC is still in its infancy. Its inherent variability may however bring to life the use of local MCMC in the same way as reversible energy-based Monte Carlo has been empowered through the a priori probabilities.
In the following, Section 2 discusses the mathematical foundations underlying ECMC, namely global balance, factorization, lifting, and thinning (making statistically correct decisions with minimal effort). Section 3 reviews exact results for various Markov chains in the simplest setting of a path graph. Section 4 studies one-dimensional N-body systems Section 5 provides an overview of findings on N-particle systems in higher than one dimension. Section 6 reviews a recent proposal to apply ECMC to molecular simulations. Section 7, finally discusses prospects for event-chain Monte Carlo and other non-reversible Markov chains in the context of molecular simulation.

MARKOV CHAINS, FACTORIZATION, LIFTING, AND THINNING
The present section reviews fundamentals of ECMC. Non-reversibility modifies the positioning of MCMC from an analogue of equilibrium physics towards the realization of a non-equilibrium system with steadystate flows. Subsection 2.1 discusses transition matrices, the objects that express MCMC algorithms in mathematical terms. The balance conditions for the steady state of Markov chains are reviewed, as well as the performance limits of MCMC in terms of basic characteristics of the sample space and the transition matrix. Subsection 2.2 discusses factorization, the break-up of the interaction into terms that, via the factorized Metropolis algorithm, make independent decisions on the acceptance or the rejection of a proposed move. Subsection 2.3 reviews the current understanding of lifting, the mathematical concept that allows non-reversible Markov chains including ECMC to replace random moves by deterministic ones without modifying the stationary distribution and while remaining within the framework of (memoryless) Markov chains. Subsection 2.4 reviews thinning, the way used in ECMC to make complex decisions with minimum computations.

Transition matrices, balance conditions
For a given sample space Ω, a Markov chain consists in a time-independent transition matrix P . Its nonnegative element P ij gives the conditional probability for the configuration i to move to configuration j in one time step. The transition matrix P is stochastic-all its rows i sum up to j∈Ω P ij = 1. Commonly, the probability P ij is composed of two parts P ij = A ij P ij , where A is the a priori probability to propose a move (mentioned in Section 1), and P ij the socalled "filter" to accept the proposed move. Any term P ii in the above transition matrix is the probability to move from i to i. If non-vanishing, it may result from the probability of all proposed moves from i to other configurations j = i (by the a priori probability) to be rejected (by the filter P). In the lifted sample spaceΩ that will be introduced in Subsection 2.3, ECMC is without rejections.

Event-chain Monte Carlo
2.1.1 Irreducibility and convergence, basic properties of the transition matrix A (finite) Markov chain is irreducible if any configuration j can be reached from any other configuration i in a finite number of steps [56,Sect. 1.3]. As the matrix P t = (P t ) ij gives the conditional probability to move from i to j in exactly t time steps, an irreducible matrix has (P t ) ij > 0 ∀i, j, but the time t may depend on i and j.
The transition matrix P connects not only configurations but also probability distributions π {t−1} and π {t} at subsequent time steps t − 1 and t. By extension, the matrix P t connects the distribution π {t} with the (user-provided) distribution of initial configurations π {0} at time t = 0: The initial distribution π {(} 0) can be concentrated on a single initial configuration. In that case, π {0} is a discrete Kronecker δ-function for a finite Markov chain, and a Dirac δ-function otherwise.
An irreducible Markov chain has a unique stationary distribution π that satisfies π i = j∈Ω π j P ji ∀i ∈ Ω. Eq.
(2) allows one to define the flow F ji from j to i as the stationary probability π j to be at j multiplied with the conditional probability to move from j to i: where the left-hand side of Eq. (2) was multiplied with the stochasticity condition k∈Ω P ik = 1.
Although any irreducible transition matrix has a unique π, this distribution is not guaranteed to be the limit π {t} for t → ∞ for all initial distributions π {0} . But even in the absence of convergence, ergodicity follows from irreducibility alone, and ensemble averages i∈Ω O(i)π i of an observable O agree with the time averages lim t→∞ [56,Theorem 4.16]).
Convergence towards π of an irreducible Markov chain requires that it is aperiodic, that is, that the return times from a configuration i back to itself {t ≥ 1 : (P t ) ii > 0} are not all multiples of a period larger than one. For example, if the set of return times is {2, 4, 6, . . . }, then the period is 2, but if it is {1000, 1001, 1002, . . . }, then it is 1. These periods do not depend on the configuration i. For irreducible, aperiodic transition matrices, P t = (P t ) ij is a positive matrix for some fixed t, and MCMC converges towards π from any starting distribution π {0} . The existence and uniqueness of the stationary distribution π follows from the irreducibility of the transition matrix, but if its value is imposed (for example to be the Boltzmann distribution or the diagonal density matrix [49]), then the steady-state flow rate of Eq. (2) becomes a necessary condition for the transition matrix called "global balance". For ECMC, this global-balance condition of Eq. (2) must be checked for all members of the lifted sample spaceΩ.

Event-chain Monte Carlo
A reversible transition matrix is one that satisfies the detailed-balance condition Detailed balance implies global balance (Eq. (4) yields Eq. (2) by summing over j, considering that j∈Ω P ij = 1), and the flow into a configuration i coming from a site j goes right back to j. In ECMC, the reality of the global-balance condition is quite the opposite, because the entering flow F ji is compensated by flows to other sites than j (F ji > 0 usually implies F ij = 0). Checking the global balance condition is more complicated than checking detailed balance, as it requires monitoring all the configurations j that contribute to the flow into i, in an augmented sample space.
For a reversible Markov chain, the matrix A ij = π 1/2 i P ij π −1/2 j is symmetric, as trivially follows from the detailed-balance condition. The spectral theorem then assures that A has only real eigenvalues and that its eigenvectors form an orthonormal basis. The transition matrix P has the same eigenvalues as A, and closely related eigenvectors: The eigenvectorsx of P must be multiplied with √ π to be mutually orthogonal. They provide a basis on which any initial probability distribution π {0} can be expanded. An irreducible and aperiodic transition matrix (reversible or not) has one eigenvalue λ 1 = 1, and all others satisfy |λ k | < 1 ∀k = 1.
A non-reversible transition matrix may belong to a mix of three different classes, and this variety greatly complicates their mathematical analysis. P may have only real eigenvalues and real-valued eigenvectors, but without there being an associated symmetric matrix A, as in Eq. (5) (see [88,Sect. 2.3] for an example). A non-reversible transition matrix may also have only real eigenvalues but with geometrical multiplicities that not all agree with the algebraic multiplicities. Such a matrix is not diagonalizable. Finally, P may have real and then pairs of complex eigenvalues [88]. This generic transition matrix can be analyzed in terms of its eigenvalue spectrum, and expanded in terms of a basis of eigenvectors [83,29,17]. As nonreversible transition matrices may well have only real eigenvalues, the size of their imaginary parts does not by itself indicate the degree of non-reversibility [74]. Reversibilizations of non-reversible transition matrices have been much studied [24], but they modify the considered transition matrix.

Total variation distance, mixing times
In real-world applications, irreducibility and aperiodicity of a Markov chain can usually be established beyond doubt. The stationary distribution π of a transition matrix constructed to satisfy a global-balance condition is also known explicitly. The time scale for the exponential approach of π {t} to π-that always exists-is much more difficult to establish. In principle, the difference between the two distributions is quantified through the total variation distance: The distribution π {t} depends on the initial distribution, but for any choice of π {0} , for an irreducible and aperiodic transition matrix, the total variation distance is smaller than an exponential bound Cα t with α ∈ (0, 1) (see the convergence theorem for Markov chains [56,Theorem 4.9]). At the mixing time, the distance from the most unfavorable initial configuration, drops below a certain value ǫ t mix (ǫ) = min{t : d(t) ≤ ǫ}.
Usually, ǫ = 1 4 is taken, with t mix = t mix (1/4) [56]. Here, the value of ǫ is arbitrary, but smaller than 1/2 in order for convergence in a small part A of Ω (without exploration of Ω \ A) not to be counted as a success (an example is given in Subsection 3.2.1). Once such a value smaller than 1 2 is reached, exponential convergence in t/t mix sets in [56, eq. (4.36)]. The mixing time is thus needed to obtain a first sample of the stationary distribution from a most unfavorable initial configuration. It does not require the existence of a spectrum of the transition matrix P , and it can be much larger than the correlation time t corr given by the absolute inverse spectral gap of P , if it exists (the time it takes to decorrelate samples in the stationary regime). Related definitions of the mixing time take into account such thermalized initial configurations [60].

Diameter bounds, bottleneck, conductance
An elementary lower bound for the mixing time on a graph G = (Ω, E) (with the elements of Ω as vertices, and the non-zero elements of the transition matrix as edges) is given by the graph diameter L G , that is, the minimal number of steps it takes to connect any two vertices i, j ∈ Ω. The mixing time, for any ǫ < 1/2, trivially satisfies For the Metropolis or heatbath single-spin flip dynamics in the Ising model with N spins (in any dimension d), mixing times throughout the high-temperature phase are logarithmically close to the diameter bound with the graph diameter L G = N, the maximum number of spins that differ between i and j [65].
Mixing and correlation times of a MCMC algorithm can become very large if there is a bottleneck in the sample space Ω from which it is difficult to excape. Two remarkable insights are that in MCMC there is but a single such bottleneck (rather than a sequence of them), and that the mixing time is bracketed (rather than merely bounded from below) by functions of the conductance [13]. Also called "bottleneck ratio" [56], the conductance is defined as the flow across the bottleneck: where S = Ω \ S. Although it can usually not be computed in real-world applications, the conductance is of importance because the liftings which are at the heart of ECMC leave it unchanged.
For reversible Markov chains, the correlation time is bounded by the conductance as [13]: The lower bound follows from the fact that to cross from S into S, the Markov chain must pass through the boundary of S, but this cannot happen faster than through direct sampling within S, that is, with probability π i /π S for i on the boundary of S. The upper bound was proven in [84,Lemma 3.3]. For arbitrary MCMC, one has the relation where A is the "set time", the maximum over all sets S of the expected time to hit a set S from a configuration sampled from π, multiplied with the stationary distribution π(S). In addition, the mixing time, defined in [13] with the help of a stopping rule, satisfies: where the constants are different for reversible and for non-reversible Markov chains.
The above inequalities strictly apply only to finite Markov chains, where the smallest weight π 0 is well-defined. A continuous system may have to be discretized in order to allow a discussion in its terms. Also, in the application to ECMC, which is event-driven, mixing and correlation times may not reflect computational effort, which roughly corresponds to the number of events, rather than of time steps. Nevertheless, the conductance yields the considerable diffusive-to-ballistic speedup that may be reached by non-reversible liftings [13] if the collapsed Markov chain is itself close to the ∼ 1/Φ 2 upper bound of Eqs (11) and (13).

Factorization
In generic MCMC, each proposed move is accepted or rejected based on the change in total potential that it entails. For hard spheres, a move in the Metropolis algorithm [49,Chap. 2] can also be pictured as being accepted "by consensus" (no pair of spheres presenting an overlap) and rejected "by veto" (at least one overlap between two spheres). The "potential landscape" and the "consensus-veto" interpretations are here equivalent as the pair potential of two overlapping hard spheres is infinite, and therefore also the total potential.
The factorized Metropolis filter [70,22] generalizes the "consensus-veto" picture to an arbitrary system whose stationary distribution breaks up into a set M of factors M = (I M , T M ). Here, I M is an index set and T M a type (or label), such as "Coulomb", "Lennard-Jones", "Harmonic" (see [22,Sect. 2A]). The total potential U of a configuration c then writes as the sum over factor potentials U M that only depend on the factor configurations c M : and the stationary distribution appears as a product over exponentials of factor potentials: Energy-based MCMC considers the left-hand side of Eq. (14) and molecular dynamics its derivatives (the forces on particles). All factorizations M are then equivalent. In contrast, ECMC concentrates on the right-hand side of Eq. (14), and different factorizations now produce distinct algorithms.

Factorized Metropolis algorithm, continuous-time limit
The Metropolis filter [67] is a standard algorithm for accepting a proposed move from configuration c to configuration c ′ : where The factorized Metropolis filter [70] plays a crucial role in ECMC. It inverts the order of product and minimization, and it factorizes as its name indicates: The factorized Metropolis filter satisfies detailed balance. This is because, for a single factor, P Fact reduces to P Met (which obeys detailed balance [49]) and because the Boltzmann weight of Eq. (15) and the factorized filter of Eq. (17) factorize along the same lines. The detailed-balance property is required for proving the correctness of ECMC (see [22, eq. (27)]), although ECMC is never actually run in reversible mode.
The Metropolis filter P Met of Eq. (16) is implemented by sampling a Boolean random variable: where "True" means that the move c → c ′ is accepted. In contrast, the factorized Metropolis filter appears as a conjunction of Boolean random variables: The left-hand side of Eq. (19) is "True" if and only if all the independently sampled Booleans X M on the right-hand side are "True", each with probability min [1, exp (−β∆U M )]. The above conjunction thus generalizes the consensus-veto picture from hard spheres to general interactions, where it is inequivalent to the energy-based filters.
In the continuous-time limit, for continuously varying potentials, the acceptance probability of a factor M becomes: where is the unit ramp function. The acceptance probability of the factorized Metropolis filter then becomes and the total rejection probability for the move turns into a sum over factors: In ECMC, the infinitesimal limit is used to break possible ties, so that a veto can always be attributed to a unique factor M and then transformed into a lifting. The MCMC trajectory in between vetoes appears deterministic, and the entire trajectory as piecewise deterministic [11,10] (see also [69]). In the eventdriven formulation of ECMC, rather than to check consensus for each move c → c ′ , factors sample candidate events (candidate vetoes), looking ahead on the deterministic trajectory. The earliest candidate event is then realized as an event [77].

Stochastic potential switching
Whenever the factorized Metropolis filter signals a consensus, the system appears as non-interacting, while a veto makes it appear hard-sphere like. This impression is confirmed by a mapping of the factorized Metropolis filter onto a hamiltonian that stochastically switches factor potentials [63], in our case between zero and infinity.
In stochastic potential switching, a factor potential U M (see [63,Sect with c † = c. The switching affects both configurations, but is done with probability S M (c) [63]. The constant ∆U * M is chosen so that S M < 1 for both c and c ′ , for example as ∆U * M = max [U M (c), U M (c ′ )]+ ǫ with ǫ 0. If the potential is not switched, the move c → c ′ is subject to a pseudo-potential for both configurations c † ∈ {c, c ′ }. ECMC considers "zero-switching" towardsŨ (c) =Ũ (c ′ ) = 0. If U M (c ′ ) < U M (c), the zero-switching probability is ∼ 1 − βǫ → 1 for ǫ → 0.
If U M (c ′ ) > U M (c), the zero-switching probability is smaller than 1. If the potential is not switched so that the pseudo-potential U is used at c and c ′ , then U M (c) remains finite while U M (c ′ ) ∼ − log (βǫ) /β → ∞ for ǫ → 0. In that limit, the pseudo-potential barrier U (c ′ ) − U (c) diverges towards the hard-sphere limit.
with the unit ramp function of Eq. (21). In the infinitesimal limit U M (c ′ ) − U M (c) → dU M ∀M ∈ M, at most one factor fails to zero-switch. The factorized Metropolis filter of Eq. (22) follows naturally.

Potential invariants, factor fields
Invariants in physical systems originate in conservation laws and topological or geometrical characteristics, among others. Potentials V that are invariant under ECMC time evolution play a special role if they can be expressed as a sum over factor potentials V M . The total system potential then writes as

Event-chain Monte Carlõ
Although V is constant, the factor terms V M can vary between configurations. In energy-based MCMC and in molecular dynamics, such constant terms play no role. Subsection 4.3.3 reviews linear factor invariants that drastically reduce mixing times and dynamical critical exponents in one-dimensional particle systems.

Lifting
Generic MCMC proposes moves that at each time are randomly sampled from a time-independent set. Under certain conditions, the same moves can be proposed in a less-random order without changing the stationary distribution. Lifting allows one to formulate the resulting algorithm as a random process without memory, that is, a Markov chain. Sections 3 and 4 will later review examples of MCMC, including ECMC, where lifting substantially reduces mixing and autocorrelation times and their scaling with system size.
Mathematically, a lifting of a Markov chain (with its opposite being a "collapsing") consists in a mapping f from a "lifted" sample spaceΩ to the "collapsed" one, Ω, such that each configuration v ∈ Ω splits ("lifts") into copies i ∈ f −1 (v). There are two requirements [13,Sect. 3]. The total stationary probability of all lifted copies must equal the original stationary probability: Also, the sum of the flows between the lifted copies of two original configurations u and v must equal the flow between the collapsed configurations: ECMC usually considers a special lifted sample spaceΩ = Ω × L, with L a set of lifting variables σ. The liftings i of a collapsed configuration v then write as i = (v, σ), so that each collapsed configuration has the same number of copies. In addition, the lifting variables usually are without influence on the stationary distribution so that:π A lifted Markov chain cannot have a larger conductance than its collapsed counterpart, because the statistical weights and the flows in Eqs (28) and (29) remain unchanged, and therefore also the bottlenecks in Eq. (10). Also, from Eq. (29), the reversibility of the lifted Markov chain implies reversibility of the collapsed chains. In this case, lifting can only lead to marginal speedup [13]. Conversely, a nonreversible collapsed Markov chain corresponds to a non-reversible lifted Markov chain. Within the bounds of Eq. (12) and the corresponding expression for mixing times, the lifted Markov chain may be closer to the ∼ 1/Φ lower "ballistic" mixing-time limit, for the collapsed Markov chain at the ∼ 1/Φ 2 upper "diffusive" limit.

Sequential MCMC, particle lifting
The sequential version of Metropolis MCMC for N-particle systems consists in performing singleparticle moves with a consecutive rather than a random order of particle indices. It was proposed in 1953, in the founding publication of MCMC, which states ". . . we move each of the particles in succession . . . " [68, p.22]. Sequential MCMC is a non-reversible "particle" lifting of the Metropolis algorithm or any other reversible Markov chain with single-particle moves, with a lifted sample spaceΩ = Ω×{1, . . . , N}, where Ω = {x = (x 1 , . . . , x N )} is the collapsed sample space of particle positions. With the lifted configurations (x, k), with k the active particle, and with stationary distributionπ(x, k) = π(x)/N, the lifted transition matrix iŝ where the periodicity in the indices is understood. The factor N on the right-hand side of Eq. (31) compensates for the a priori probability of choosing the active particle k, which is 1 N in the collapsed MCMC and 1 in the lifted MCMC. Sequential MCMC satisfies the global-balance condition if the collapsed Markov chain is reversible (see Subsection 4.2.2 for an incorrect sequential version of a nonreversible MCMC).
In particle systems with central potentials, sequential MCMC is found to be slightly more efficient than its collapsed reversible counterpart [76,44,53]. Likewise, in the Ising model and related systems, the analogous spin sweeps (updating spin i + 1 after spin i, etc.) again appear marginally faster than the random sampling of spin indices [5]. Sequential MCMC is but a slightly changed variant of a reversible Markov chain, yet its particle lifting, the deliberative choice of the active particle, is one of the key features of ECMC.

Displacement Lifting
In the local MCMC algorithms which are the collapsed counterparts of ECMC, the proposed moves can (roughly) be characterized by the independent sampling of a particle (or spin) i ∈ {1, . . . , N} and of a displacement δ from a displacement set D. The "displacement lifting" refers to a less-random sampling of the elements of D in the lifted transition matrix. For a single particle, on a one-dimensional path graph P n displacement lifting, with D = {−1, +1} corresponding to positive and to negative hops so thatΩ = Ω × D), is reviewed in Subsection 3.1.1 (see [18]). The lifting scheme, the choice of the lifted transition matrix, strives to preserve δ from one move to the next as much as is compatible with global balance and aperiodicity. Displacement lifting applies in more than one dimensions, for example for systems of dipole particles [79], whose dynamics is profoundly modified with respect to that of its collapsed counterpart. The displacement lifting in ECMC repeatedly uses the same infinitesimal displacement.

Event-chain Monte Carlo
ECMC [8,70], in its use of factorized potentials and the factorized Metropolis algorithm, realizes the extreme limit of particle lifting and of displacement lifting. If compatible with global balance and irreducibility, it reutilises its moves, that is, keeps the active particles and their displacements from one time step to the next. ECMC moves are thus persistent, and their randomness is strongly reduced, rendering its trajectories piecewise deterministic (compare with Subsection 2.2.1). The persistence of moves requires the global-balance condition to be carefully checked at every infinitesimal time step.
In ECMC, factorization and the consensus principle of the factorized Metropolis algorithm ensure that a single factor triggers the veto that terminates the effectively non-interacting trajectory (as discussed in Subsection 2.2.2). This factor is sole responsible for the ensuing lifting. The triggering factor has at least one active particle in its in-state. The out-state can be chosen to have exactly the same number of active elements. In translationally invariant particle systems (and in similarly invariant spin systems), the lifting does not require a change of displacement vector δ, and concerns only the particle sector of the lifting set L. For a two-particle factor, the lifting scheme is deterministic, and the active and the passive particle interchange their roles. For a three-particle factor, the active particle may have to be sampled from a unique probability distribution, so that the lifting scheme is stochastic but uniquely determined [30]. For larger factors, the next active particle may be sampled following many different stochastic lifting schemes that can have inequivalent dynamics [30] (see also [22, Sect. II.B]).

Thinning
Thinning [57,16] denotes a strategy where an event in a complicated time-dependent random process (here the decision to veto a move) is first provisionally established on the basis of an approximate "fatter" but simpler process, before being thinned (confirmed or invalidated) in a second step. Closely related to the rejection method in sampling [49], thinning doubly enables all but the simplest ECMC applications. First, it is used when the factor event rates β [dU M ] + of Eq. (23) cannot be integrated easily along active-particle trajectories, while appropriate bounding potentials allow integration. Second, thinning underlies the cellveto algorithm [43], which ascertains consensus and identifies vetoes among ∼ N factors (the number of particles) with an effort that remains independent of N. The random process based on the untrimmed bounding potential is time-independent and pre-computed, and it can be sampled very efficiently.
Thinning thus uses approximate bounding potentials for the bulk of the decision-making in ECMC while, in the end, allowing the sampling to remain approximation-free. This contrasts with molecular dynamics, whose numerical imperfections-possible violations of energy, momentum, and phase-space conservation under time-discretization-cannot be fully repaired (see [93]). ECMC also eliminates the rounding errors for the total potential, that cannot be avoided in molecular dynamics for long-range Coulomb interactions.

Thinning and bounding potentials
Thinning is unnecessary for hard spheres and for potentials with simple analytic expressions that can be integrated analytically. In other cases, it may be useful to compare a factor event rate to a (piecewise) upper bound, as In the case where, for simplicity, the active particle k moves along the x direction, the bounding factor event rate allows one to write The factor M triggers the veto with the factor event rate on the left-hand side of this equation, which varies with the configuration. On the right-hand side, this is expressed as a product of independent probabilities corresponding to a conjunction of random variables, of which the first corresponds to a constant-rate Poisson process and the second to a simple rejection rate. The second factor must be checked only when the first one is "True". Thinning thus greatly simplifies ECMC when the mere evaluation of U M and its derivatives is time-consuming. The JeLLyFysh application of Subsection 6.2 implements thinning for potentials in a variety of settings (see [35]).

Thinning and the cell-veto algorithm, Walker's algorithm
The cell-veto algorithm [43] (see also [12]) permanently tiles two or three-dimensional sample space (with periodic boundary conditions) into cells with precomputed cell bounds for the factor event rates of any long-ranged factor that touches a pair of these cells. In this way, for example, a long-range Coulomb event rate for a pair of distant atoms (or molecules) can be bounded by the cell event rate of the pair of cells containing these atoms. For an active atom in a given cell, the total event rate due to other faraway atoms can also be bounded by the sum over all the cell event rates. Identifying the vetoing atom (factor) then requires to provisionally sample a cell among all those contributing to the total cell event rate. This problem can be solved in O(1) using Walker's algorithm [87,16]. After identification of the cell, the confirmation step involves at most a single factor. The cell-veto algorithm is heavily used the JeLLyFysh application (see Subsection 6.2). In essence, the cell-veto algorithm thus establishes consensus and identifies a veto among ∼ N factors in O(1) in a way that installs neither cutoffs nor discretizations.

SINGLE PARTICLES ON A PATH GRAPH
The present section reviews liftings for one-particle MCMC on a path graph P n = (Ω n , E n ), where the sample space is Ω n = {1, . . . , n}, and where the set of edges E n = {(1, 2), . . . , (n − 1, n)} indicates the non-zero elements of the transition matrix. The graph P n thus forms a one-dimensional n-site lattice without periodic boundary conditions. The stationary distribution is {π 1 , . . . , π n }. Two virtual vertices 0 and n + 1 and virtual edges (0, 1) and (n, n + 1), with π 0 = π n+1 = 0, avoid the separate discussion of the boundaries. The Metropolis algorithm on P n (with the virtual additions, but an unchanged sample space) proposes a move from vertex i ∈ Ω n to j = i ± 1 with probability 1 2 . This move is accepted with the Metropolis probability min(1, π j /π i ). If rejected, the particle remains at i. The equilibrium flows: satisfy detailed balance, and the total flow into each configuration i (including from i itself) satisfies the global-balance condition of Eq. (2), as already follows from detailed balance.
The displacement-lifted Markov chains in this section, withΩ n = Ω n × {−1, +1}, split the Metropolis moves i → i±1 into two copies (the virtual vertices and edges are also split). The lifted MCMC algorithm consists in the succession of a transport transition matrixP trans and a re-sampling transition matrixP res . P trans proposes a move from vertex (i, σ) to vertex (i + σ, σ) with probability 1 and accepts with the Metropolis filter min(1, π i+σ /π i ) given thatπ (k,σ) = 1 2 π k ∀k. In the "+1" copy, this corresponds to forward moves, and in the "−1" copy to backward moves. All rejections are replaced by liftings: If the move is not accepted, a "lifting" move (i, σ) → (i, −σ) takes place.P trans , which is not a Metropolis filter, satisfies the global-balance condition for each lifted configuration, as the flow into (i, σ) equals 1 2 π i =π (i,σ) : transport flows: In the re-sampling stage, the configuration (i, σ) is simply switched to (i, −σ) with a small probability ǫ = 1/n.P res satisfies detailed balance, and the re-sampling flows are: re-sampling flows: The lifted transition matrixP =P transP res satisfies global balance. In ECMC applications, the resampling stage after each transport stage is abandoned in favor of a carefully designed re-sampling after a large number of transport steps.
The speedups that can be achieved by lifting on the path graph P n depend on the choice of π. Bounded stationary distributions are reviewed in Subsection 3.1) and unbounded ones in Subsection 3.2. All Markov chains on P n naturally share a diameter bound n which is independent of π and that (up to a constant) is the same for the lifted chains.

Bounded one-dimensional distributions
In real-world applications of MCMC, many aspects of the stationary distribution π remain unknown. Even the estimation of the minimum or the maximum of π i is often a non-trivial computational problem treated by simulated annealing, a variant of MCMC [46,49]. Also, crucial characteristics as the conductance are usually unknown, and the mixing time can at best be inferred from the running simulation. Simplified models, as a particle on a path graph, offer more control. This subsection reports on results for stationary distributions for which (max i∈Ωn π i )/(min j∈Ωn π j ) remains finite for n → ∞. The conductance bound then scales as Φ ∼ n for n → ∞, and reversible local MCMC mixes in O(n 2 ) steps.

Constant stationary distribution, transport-limited MCMC
The constant stationary distribution π = { 1 n , . . . , 1 n } on P n (see [18]) illustrates the lifting concept. The collapsed Metropolis MCMC moves with probability 1 2 to the left and to the right. This Markov chain is aperiodic because of the rejections triggered by the aforementioned virtual vertices and edges (which are not counted in P n ). The collapsed Markov chain performs a diffusive random walk. Its O(n 2 ) mixing time is close to the upper mixing-time bound of Eq. (13). As t mix is on larger time scales than the lower conductance bound, the collapsed Markov chain appears "transport limited". Reorganizing its moves through lifting is thus promising.
In the lifted sample spaceΩ n = Ω n × {+, −}, the transition matrixP trans describes a deterministic rotation on the lifted graph which is equivalent to a ring with 2n vertices and periodic boundary conditions induced by the lifting moves (i = 1, σ = −1) → (1, 1) and (i = n, σ = +1) → (n, −1) at the boundaries (see Eq. (34)). Randomness is introduced by the re-sampling with rate ǫ = 1 n withP res , which switches the copies. The mixing time of the combined transition matrixP transP res is t mix = O(n) (see [18]).
The re-sampling ofP res is not the only strategy for preserving some randomness in the model. Restarts provide an alternative to small-probability re-sampling at each time step, where after a given number ℓ of moves with one type of lifting variables, the type is changed. Restarts can be formulated as Markov chains, with a lifting variable that performs a countdown. In ECMC, restarts are often required to ensure irreducibility.

Square-wave stationary distribution, inefficient liftings
On the path graph P n with the square-wave stationary distribution π 2k−1 = 2 3n , π 2k = 4 3n ∀k ∈ {1, . . . , n/2} (for even n), the conductance Φ = 2 3n for n → ∞ is of the same order of magnitude as for the flat distribution. Its bottleneck (the vertex where the conductance is realized), is again at i = n/2. The Metropolis algorithm proposes moves from vertex i with probabilities 1 2 to i ± 1 but rejects them with probability 1 2 if i is odd (see Eq. (16)). Its mixing time is O(n 2 ), on the same order as in Subsection 3.1.1. In the lifted sample spaceΩ = Ω × {+, −}, the transition matrixP trans generates lifting moves with probability 1 2 on odd-numbered vertices, and the persistence of the lifting variable is lost. In consequence, the lifting i → (i, ±1) is inefficient for the square-wave stationary distribution, and the mixing time scale unchanged as t mix = O(n 2 ). An additional set of height-type lifting variables, that decomposeπ into a constant "lane" and an another one that only lives on even sites, recovers O(n) mixing [33].
Lifting on the path graph is thus more difficult to put into practice for the square-wave than for the homogeneous distribution, illustrating lifted mixing-time exponents are sensitive to the details of π. ECMC in real-world applications is likewise simpler for homogeneous systems with a global translational invariance (or with a similar invariant for spin systems). Remarkably however, ECMC can apply the same rules to any continuous interparticle potential in homogeneous space. It achieves considerable speedup similar to what the lifting of Subsection 3.1.1 achieves for the flat stationary distribution on the path graph.

Unbounded one-dimensional stationary distributions
For unbounded ratios of the stationary distribution (max ij (π i /π j ) → ∞ for n → ∞) the conductance limits as well as the benefits that can be reaped through lifting vary strongly. Exact results for a V-shaped distribution and for the Curie-Weiss model are reviewed.

V-shaped stationary distribution, conductance-limited MCMC
The V-shaped stationary distribution on the path graph P n of even length n is given by π i = const| n+1 2 − i| ∀i ∈ {1, . . . , n}, where const = 4 n 2 (see [34]). The stationary distribution π thus decreases linearly from i = 1 to the bottleneck i = n 2 , with π( n 2 ) = π( n 2 + 1) = const 2 ∼ n −2 and then increases linearly from i = n 2 + 1 to i = n. The Metropolis algorithm again proposes a move from vertex i to i ± 1 with probabilities 1 2 , and then accepts them with the filter of Eq. (16) (with the virtual end vertices ensuring the correct treatment of boundary conditions). The conductance equals Φ = 2 n 2 , for the minimal subset S = {1, . . . , n/2} (see Eq. (10)). The Metropolis algorithm mixes in S on an O(n 2 ) diffusive timescale, but requires O(n 2 log n) time steps to mix in Ω [34]. However, even a direct sampling in S, that is, perfect equilibration, samples the boundary vertex n/2 only on a π i /π S ∼ n −2 inverse time scale. For the Vshaped distribution, the benefit of lifting is thus at best marginal (from O(n 2 log n) to O(n 2 )), as the collapsed Markov chain is already conductance-limited, up to a logarithmic factor.
The optimal speedup for the V-shaped distribution is indeed realized withΩ = Ω × {+, −}, and the transition matricesP trans andP res . The lifted Markov chain reaches a total variation distance of 1 2 and mixes in S on an O(n) timescale. A total variation distance of ǫ < 1 2 , and mixing inΩ is reached in O(n 2 ) steps only [34]. This illustrates that ǫ < 1 2 is required in the definition of Eq. (7).

Curie-Weiss model, mixing times vs correlation times
The Curie-Weiss model (or mean-field Ising model [4]) for N spins s i = ±1 has the particularity that its total potential only depends on the magnetization i s i = M. The mapping k = (M + N)/2 + 1, n = N + 1, places it on the path graph P n . The heat-bath algorithm (where a single spin is updated) can be interpreted as the updating of the global magnetization state. The model was studied in a number of papers [55,86,23,11].
The conductance bound of Eq. (10) for the Curie-Weiss model derives from a bottleneck at k = n/2. In the paramagnetic phase, one has Φ = O(n 1/2 ), and at the critical point βJ = 1, Φ = O(n 3/4 ) (see [55]). These bounds are smaller than n because the magnetization is strongly concentrated around M = 0 (that is, k = n/2), in the paramagnetic phase. In the ferromagnetic phase, βJ > 1, the magnetizations are concentrated around the positive and negative magnetizations (k 1 and k n) and in consequence Φ grows exponentially with n, and so do all mixing times.
The reversible Metropolis algorithm mixes in n log n steps for T > T c (see [55]). This saturates the upper bound of Eq. (13). At the critical point, the mixing time is as n 3/2 at T = T c . As is evident from the small conductance bound, the Metropolis algorithm is again limited by slow diffusive transport and not by a bottleneck in the potential landscape.
The lifted algorithm improves mixing times in both cases to the optimal. The correlation time can be shorter than n because of the strong concentration of the probability density on a narrow range of magnetizations [86,23,11], while mixing time on the path graph is always subject to the ∼ n diameter bound. By contrast, in the paramagnetic phase J > 1, the conductance bound is exponential in n, and all mixing times are limited by the central bottleneck and not by slow transport.

N PARTICLES IN ONE DIMENSION
The present section reviews MCMC for N hard-sphere particles on a one-dimensional graph with periodic boundary conditions (the path graph P n with an added edge (n, 1)), and on continuous intervals of length L with and without periodic boundary conditions. In all cases, moves respect the fixed order of particles (x 1 < x 2 , . . . , < x N −1 < x N , possibly with periodic boundary conditions in positions and particle indices). With periodic boundary conditions, uniform rotations of the configuration are ignored, as they mix very slowly [80]). The hard-sphere MCMC dynamics is essentially independent of the density for the discrete cases (see [50, eq. (2.18]) as well as in the continuum (see [44] and [53, Fig. 1]). Most of the hard-sphere results are numerically found to generalize to a class of continuous potentials (see [44,Supp. Item 5]). Subsection 4.1 reviews exact mixing and correlation-time results for reversible Markov chains and Subsection 4.2 those for non-reversible ones, including the connection with the totally asymmetric simple exclusion model (TASEP). Subsection 4.3 discusses particle-lifted Markov chains, like the lifted TASEP, the lifted-forward Metropolis algorithm as well as ECMC, for which rigorous mixing times were obtained. In many cases, non-reversible, and in particular non-reversible lifted Markov chains and ECMC can mix on faster time scales than their collapsed counterparts.

Reversible MCMC in one-dimensional N-particle systems
Although local hard-sphere MCMC was introduced many decades ago [67], its mixing times were obtained rigorously only in recent years, and this only in one spatial dimension. In the discrete case (the symmetric simple exclusion process), the especially precise description of the convergence process calibrates numerical approaches.

Reversible discrete one-dimensional MCMC
The symmetric simple exclusion process implements the local hard-sphere Metropolis algorithm in the discrete sample space Ω SSEP = {x 1 < x 2 <, . . . , < x N } with x i ∈ {1, . . . , n} (periodic boundary conditions understood for positions and indices). All legal hard-sphere configurations c have the same weight π c = π * . At each time step t, a random move x t i →x i = x t i ± 1 is proposed for a (single) randomly sampled particle i. If vertexx i is unoccupied, the move is accepted (x t+1 i =x i ). Otherwise the particle stays put (x t+1 i = x t i ). The process is time-reversible and it thus trivially satisfies detailed balance. The equilibrium flows into a legal configuration c = (x 1 , . . . , x N ) arrive from the configuration c itself and possibly from 2N neighboring configurations that may however not all be legal [44]: The configuration c + i contributes flow Total flow into c: The mixing time of the symmetric simple exclusion process from the most unfavorable initial distribution π {0} (the compact configuration) is known to scale as O (N 3 log N), whereas the mixing time from an equilibrium configuration scales only as O(N 3 ), that is, as the correlation time [50]. These behaviors are recovered by numerical simulation.

Continuous one-dimensional reversible MCMC
In the one-dimensional continuum, the scaling for the mixing times of reversible local hard-sphere MCMC has been obtained for the heat-bath algorithm. In this dynamics, at each time step, the position of a randomly sampled particle i is updated to a random value in between particles i−1 and i+1. The heatbath algorithm mixes in O(N 3 log N) moves [81] on an interval with fixed boundary conditions. The mixing time for the same model with periodic boundary conditions is between O(N 3 ) and O(N 3 log N) [80]. Numerical computations confirm O(N 3 log N) mixing on the continuous interval with periodic boundary conditions both for the heatbath and for the Metropolis algorithm [44].

Non-reversible MCMC in one-dimensional N-particle systems
Non-reversible MCMC in one-dimensional particle systems has a long history in physics and mathematics through the study of the totally symmetric simple exclusion process (TASEP) [14], which is one half of a directional lifting of the simple symmetric exclusion process. Even for that model, however, rigorous mixing-time results are very recent. The "particle-lifted" TASEP [44] is a variant of the TASEP and, at the same time, a finite-displacement variant of ECMC. With its chains of particles moving forward, it features intriguing mixing behavior. This introduces to the problem of restarts, a key subject for ECMC.
All systems considered in this subsection have periodic boundary conditions. In the setting where the collapsed algorithm features forward and backward moves, one may access the limit of vanishing switching rates, that is, retain only one of the equivalent sectors of the lifted transition matrix.

Discrete non-reversible MCMC in one dimension -TASEP, lifted TASEP
The TASEP is a displacement lifting of the simple symmetric exclusion process (see Subsection 4.1.1) with a lifted sample spaceΩ = Ω SSEP × {+, −}. Because of the periodic boundary conditions, neither the global-balance condition nor the irreducibility condition require liftings between the "+" copy, where particles attempt moves as x i → x i + 1 and the "−" copy, where particles attempt moves as x i → x i − 1. It suffices to consider only one half of the lifted sample space, say, the "+" copy.
The TASEP attempts at each time step to advance a random particle i from its current vertex x t i tõ x = x t i +1 (with periodic boundary conditions), in other words from c to c − i . If the vertexx i is unoccupied (if c − i is legal), the move is accepted (x t+1 i =x i ), and otherwise the particle stays put.
The equilibrium flow into a configuration c = (x 1 , . . . , x N ) now arrives from c itself and possibly from the N neighboring conditions c + i of Eq. (37). The configuration c + i contributes accept-flow A + i = 1 N π c + i which is 1 N π * if c + i is legal, and zero otherwise (compare with Subsection 4.1.1). Now, c itself contributes reject-flow of particle i − 1, namely R + i−1 , which is 1 N π * if c + i is illegal and zero otherwise. The total flow into configuration c (with periodic boundary conditions understood) is therefore [44]: Total flow into c: The TASEP is proven to mix in O(N 5/2 ) (see [3]). The equilibrium correlation time (the inverse gap of the transition matrix that, in this case, can be diagonalized) has the same O(N 5/2 ) scaling [17,29]. Numerical simulations (with an operational definition of mixing times [44, Supplementary Item S2]) recover this behavior.
In Eq. (39), the rejected moves of particle i − 1 compensate the accepted moves of particle i, because in the TASEP they are both sampled with the same 1 N probability. For this reason, a naive sequential TASEP, where particles are simply chosen one after the other, is incorrect. The flow into a given lifted configuration (c, i) is given by terms analogous to A + i and to R + i (see Eq. (39), but without the factors 1 N ). They do not add up to a constant π * . An improved lifting scheme however remedies the situation. The particlelifted TASEP [44] concerns a lifted sample spaceΩ = Ω SSEP × {1 + , . . . , N + ; 1 − , . . . , N − }. Again, the periodic boundary conditions allow the "+" and "−" sectors of lifting variables to separate, and only one of these sectors (and of the corresponding transition matrix) must be considered, say, the "+" sector. In the particle-lifted TASEP, the lifted configuration (c, i + ) attempts to move to (c − i , i + ) (by advancing i by one vertex). If c − i is illegal, the particle lifting (c, i + ), . . . , (c, (i + 1) + ) takes place. The flows into a configurationĉ = (c, i + ) is contributed to by two possible configurations, namely by (c, (i−1) + and by the configuration (c + i , i) of Eq. (37). Global balance can be established analogously to Subsection 4.2.1, with the uniform stationary distribution c * . The particle-lifted TASEP is deterministic, and is not irreducible (a compact state of N particles evolves at best into a compact state of N − 1 plus one detached particle). Randomness can be introduced through restarts [44].

Continuous non-reversible MCMC in one dimension: Forward, Lifted Forward
The discrete MCMC algorithms of Subsection 4.2.1 can all be generalized to the continuous case with Ω SSEP replaced by the N-dimensional hypercube with 0 < x i < x i+1 < x N < L (with periodic boundary conditions in L implied). The moves x i → x i ± 1 of the discrete Markov chains can now be generalized to x i → x i ± ǫ, where ǫ > 0 is taken from a certain probability distribution and the non-hopping condition x i < x i+1 (corrected for periodic boundary conditions in N) is enforced.
With the non-hopping condition, the MCMC algorithms for spheres of diameter σ on a ring of length L are equivalent to those for zero-diameter point particles in a ring of length L free = L − Nσ, for which the non-hopping condition becomes the only test of legality of a configuration. The global-balance condition can be checked for every value of ǫ separately, with configurations c +ǫ i and c −ǫ i generalizing the 2N neighboring configurations c + i and c − i of Eq. (37). The value of ǫ may serve as an additional lifting variable. This has however not been studied.
The forward Metropolis algorithm, where a move by ǫ > 0 is proposed to a random sphere i, is the continuous-space version of the TASEP, and the particle-lifted forward Metropolis algorithm generalizes the particle-lifted TASEP. As this algorithm now has a random element ǫ, it is generically irreducible and aperiodic, and its mixing time is numerically found to scale as t mix = O(N 5/2 ) (see [44]). With a suitable restart strategy, this mixing time can be decreased to O( 2 log N) scaling. The sequential version of the lifted forward Monte Carlo algorithm is again incorrect, although the local relabeling strategy outlined in Subsection 4.3.2 can be applied. Extensions of the hard-sphere case to nearest-neighbor, monotonous potentials with a no-hopping conditions were discussed [44, Supp. Item 5].

ECMC in one-dimensional N-particle systems
For one-dimensional N-particle hard-sphere systems, ECMC is the infinitesimal-displacement limit (ǫ → 0) of the particle-lifted forward Metropolis algorithm (see Subsection 4.2.2). Its validity follows, in addition to this limiting procedure and to an explicit check of global balance, through the equivalence of this specific MCMC algorithm with molecular dynamics with an initial condition where one sphere has velocity +1 and all other spheres are stationary (velocity v i = 0). In the ǫ → 0 limit, the details of the probability distribution of ǫ disappear, and only its mean value remains relevant. The randomness in ǫ that ensures the irreducibility of the forward Metropolis algorithm disappears. ECMC thus requires restarts after a given chain length ℓ = ǫ t in order to achieve irreducibility. These restarts replace the re-sampling transition matrixP res of Subsection 3.1. For specific chain-length distributions (probability distributions for ℓ), the ECMC relaxation can be described exactly. Non-trivial relabelings and factor-field variants again illustrate that ECMC is a family of algorithms with different mixing-time scalings rather than a fixed dynamics.

Standard ECMC, coupon-collector problem
Hard-sphere ECMC in one dimension (with periodic boundary conditions) consists in the limit of the forward Metropolis algorithm with restarts, where the step size ǫ t at time t approaches zero, and the number l of steps in the event chain (between restarts) approaches infinity such that their product the total chain length ℓ = l−1 t=0 ǫ t , remains finite. For hard spheres, ECMC is devoid of random elements except of the chain length ℓ (which may be a random variable) and the sphere which initiates each event chain.
In the reduced-circle representation of length L free , with spheres of zero diameter (see Subsection 4.2.2), one event chain of length ℓ effectively advances the sphere that initiates the event chain by a distance ℓ. For ℓ = ran[const, const + L free ], the initial sphere of the event chain is placed (up to a relabeling) at a random position on the ring. EMCMC obtains a perfect sample when each sphere has once initiated an event chain. The solution of the mathematical coupon-collector problem shows that it takes O(N log N) random samples of N elements to touch each of them at least once. In consequence, the number of event chains required for mixing scales as ∼ N log N, with a total event count of O(N 2 log N) (see [53,51]. This number of events represents the computational effort more faithfully than the number of time steps, which diverges (by construction) for each finite-length event chain.

ECMC with local relabeling
Pair-factor ECMC may be modified by having at each event the two particles of the vetoing factor not only switch their lifting status (active / target), but also their identities (particle indices). This relabeling is without incidence on the sequence of positions during a given event chain, throughout which the initial sphere remains active. The relabeling at each event is local. It suffices to render sequential ECMC correct [53]. (The same relabeling also allows the sequential forward MCMC to satisfy global balance.) With the initial sphere of each event-chain taken sequentially rather than randomly, the coupon-collector-related logarithmic slowdown is avoided so that the choice of random chain lengths ℓ = ran[const, const + L free ], the only remaining random element in the algorithm, produces a perfect sample after N chains and O(N 2 ) events rather than the O(N 2 log N) events of regular ECMC. ECMC with local relabeling again illustrates that reducing the randomness of moves can speed up mixing, in other words, bring about faster overall randomization. In higher-dimensional hard-sphere systems, the relabeling appears to speed up the approach to equilibrium by only a constant factor [51].

Factor-field ECMC in one dimension
The factorized Metropolis filter separates the total potential into independent terms (the factors). Large variations in one factor potential then remain uncompensated by those in other factors, as would be the case for the total potential (for energy-based filters). As a consequence, ECMC may possess event rates that are too high for efficient mixing. This potential shortcoming of standard ECMC was remarked at low temperature in one-dimensional N-particle systems [36] (see also [53, Fig. 2]). It is overcome in onedimensional particle systems through the addition of invariants and their subsequent breakup into factor fields. Such terms may decrease the event rate. Moreover, they can profoundly modify the characteristics of the event-chain dynamics and decrease mixing-time and correlation-time exponents.
For N particles on a one-dimensional interval of length L, the quantity i (x i+1 − x i ) (with periodic boundary conditions in L and in N) is an invariant. A linear function f (x) = ax can transform it into a Krauth Event-chain Monte Carlo sum of factor fields: factor field (40) The factor field of Eq. (40) may be added and its linear parameter a adjusted to any pair-factor potential. Attractive factor fields may thus be added to hard-sphere factors or to Lennard-Jones factors [54]. With the linear factor adjusted to compensate the virial pressure, O(N 3/2 ) autocorrelation times (rather than O(N 2 ) without factor fields) and O(N 2 ) mixing times (rather than O(N 2 log N) are found. One particular feature of event-chain dynamics at the optimal value of the factor field is that the chains have zero linear drift, which itself is a measure of the virial pressure [54,70].

STATISTICAL-MECHANICS MODELS IN MORE THAN ONE DIMENSION
The present section reviews statistical-physics models in more than one dimensions where, in ECMC, several factors are present at any time. The vetoing factor is chosen among them. Subsection 5.1 reviews the two-dimensional hard-disk model, the first application of ECMC [8,7,6]. The long-range physics of this model is intricate, with phonons coexisting with two types of topological excitations (dislocations and disclinations). Hard-disk ECMC is contrasted with related algorithms that have been applied to this system. Subsection 5.2 briefly reviews ECMC for the harmonic model of a solid (with a fixed list of neighbors for each particle) which also serves as the fundamental low-temperature theory for continuousspin systems, as the XY or the Heisenberg model. In these spin models, the behavior of ECMC has been thoroughly analyzed, and the interplay between spin waves and topological excitations for local MCMC algorithms including ECMC, as reviewed in Subsection 5.3, is better understood than for particle systems.

Two-dimensional hard disks
ECMC is conceptually simpler for hard spheres than for general interactions, as the potential switching between zero and infinity is deterministic rather than stochastic (see Subsection 2.2.2). Most studies have concentrated on the two-dimensional case, the hard-disk model, where fundamental results were obtained using ECMC [7,21]. (The original Fortran90 implementation used is publicly available [58].) Harddisk ECMC is found to be very fast, but elementary questions as for example the scaling of correlation times with system size, in the hexatic and solid phases, are still without even an empirical answer. The simplifications related to the existence of a constraint graph of possible liftings are at the heart of a recent implementation of multithreaded ECMC for hard disks [58].

Characterization of hard-sphere ECMC
Pair-factor ECMC for two-dimensional hard disks [8] consists in freely moving an active sphere in continuous time t along a fixed direction up to a lifting event. This event corresponds to contact with a target sphere, which then exchanges its lifting status with the active one. Two versions of hard-sphere ECMC are compatible with the conservation of phase-space volume: "Straight" ECMC, where the outstate displacement is along the same directions as the in-state displacement, reducing the effective spatial dimension to one. Between restarts, when all lifting variables are reset, straight ECMC is in higher dimensions no different from the one-dimensional case treated in Section 4.3. For straight ECMC, the lifted sample space is for exampleΩ = Ω × {1 ±x , . . . , N ±x ; 1 ±y , . . . , N ±y }, where 2N-dimensional collapsed sample space Ω describes the particle positions. Periodic boundary conditions are assumed, and the "+x" and "+y" copies separate from the motions in the negative directions, which can be neglected. The correctness of straight ECMC also follows from the continuum limit of a hard-disk lattice version as well as from a mapping onto molecular dynamics with two-dimensional sphere positions but with onedimensional sphere velocities (see for example [58,Lemma 1]). Alternatively, "reflected" ECMC uses as an out-state displacement direction the one reflected with respect to the line (hyperplane) of incidence (see [6, Fig. 2.8]. Straight ECMC is generally faster [8, Fig. 6]. A variant of straight ECMC with a larger direction space D containing many angles was found to bring no improvements for hard disks [90], while it has drastic effects for hard-sphere objects with internal degrees of freedom [79].
In hard-sphere systems, the computation of the pressure is notoriously complicated for generic MCMC, which does not allow one to compute the virial from a finite number of equiprobable configurations. Instead, the virial is extracted from a high-precision extrapolation of the pair-correlation function to contact [6,Sect. 3.3.4]. Continuous-time ECMC obtains an infinite number of highly correlated samples between any two lifting events, and it allows for an improved estimator [70,Eq. 20]) for the virial pressure from the expectations of basic geometrical properties of the ECMC trajectories. The estimator was tested in two-dimensional and three-dimensional hard spheres (see [70,21,38]. It generalizes from hard spheres to arbitrary potentials. Between restarts of hard-disk ECMC, any disk can lift with no more than three other disks, if they are monodisperse (see [42]). This generalizes from the one-dimensional case, where the order between spheres remains unchanged so that a sphere lifts with only a single other sphere that it cannot hop over. A constraint graph encodes these relations. It has at most three outgoing arrows for every disk, corresponding to the liftings it can provoke (see [42]). The optimal (sparsest) constraint graph is locally planar [58,Sect. 3.1]. The constraint-graph formulation of hard-sphere ECMC exposes its close connection with the harmonic-model ECMC (see Subsection 5.2), where the neighbor relations are likewise fixed.
In two-dimensional hard-disk MCMC algorithm, the total variation distance or spectral gaps cannot be estimated or evaluated, and theoretical bounds cannot be used [19]. Rigorous mixing-time scaling exponents are available for low density, but only for a non-local version of the Metropolis algorithm [41]. The analysis of all recent computations builds on the hypothesis that for two-dimensional disks the autocorrelation function of the global orientational order parameter is the slowest relaxation process in this system. The hypothesis has proven robust. Practical computations generally adopt a square box with periodic boundary conditions, where the expectation of the above autocorrelation function vanishes because of symmetry, simplifying the interpretation of time series [8,7].

Hard-disk ECMC and other algorithms
In ECMC, the continuous-time limit is usually required in order to break ties between candidate event times of different factors. For two particles with pair interactions in a box with periodic boundary conditions, and in particular for two hard spheres, ECMC satisfies global balance with finite displacements D = {δ x , δ y } (that would have to vary in order to ensure aperiodicity). Because of translational invariance, a given event chain, say, with a fixed displacement δ x then maps to the transition matrixP trans on a path graph (see Subsection 3.1.1), with the displacement of one hard disk corresponding to the "+1" sector, and the displacement of the other corresponding ot the "−1" sector.
For any N and in any spatial dimension d, hard-sphere ECMC is equivalent to modified molecular dynamics where, for one event chain, the positions are d-dimensional but the velocities are onedimensional. This corresponds to hard spheres on one-dimensional constraining "rails" that remain fixed in between restart. The ECMC events correspond to molecular-dynamics collisions, which conserve energy and momenta for an configuration with only a single non-zero velocity (the one of active disks).
A discrete-time precursor algorithm of ECMC [39,Sect. 5] chooses for each move a random initial disk, a random direction of displacement (such that a direction and its inverse are sampled with the same probability), and a total number n c of disks to be displaced. A chain move is then constructed (not unlike ECMC) by displacing n c − 1 disks (starting with the initial one) along the direction of displacement until they hit their successor disks. The n c th (final) disk is placed randomly between its initial position and its (hypothetical) successor disk. In order to satisfy detailed balance, this algorithm requires a Metropolis rejection step for the entire chain by the ratio of the intervals available for the first and for the final disks. One move of the precursor algorithm resembles the continuous-time ECMC evolution between restarts. However, it cannot be interpreted in terms of a continuous-time Markov chain, and it requires rejections and must allow for chains in a given direction and its inverse direction with equal probability. The convergence properties of this algorithm have not been analyzed.
The hard-disk model was also successfully simulated on graphics cards (GPUs) with a massively parallel implementation of the Metropolis algorithm using a four-color checkerboard scheme. In this scheme, a cell of any one color touches only differently colored cells [2,21]). At any time, disks in cells of only one color can thus be updated in parallel. Irreducibility is assured by frequently translating the fourcolor checkerboard by a random vector. This highly parallel algorithm overcomes its considerable speed handicap (roughly two orders of magnitude in single-processor mode), through massive parallelization.
Massively parallel MCMC has been tested to high precision against ECMC and against event-driven molecular dynamics for a number of physical quantities, as the pressure (that is, the equation of state) and the orientational and positional order parameters [21].

Parallel hard-disk ECMC
Hard-disk ECMC and event-driven molecular dynamics both identify the earliest one of a number of candidate events, the one that will be realized as an event and drive the dynamics of the system. The event then generates new candidate events, while some of the old ones continue to exist and yet others disappear. Modern event-driven molecular dynamics codes are optimized for the management and the update of a very large number of such candidate events, typically organized in a heap [82,37]. In event-driven molecular dynamics, an extensive number ∝ N of candidate events are present at any given moment. The possible scheduling conflicts among this large number of candidate events has long stymied attempts to parallelize the algorithm, that is, to handle several events independently from each other [62,61,28,48,64].
Domain-decomposition strategies for the candidate events in molecular dynamics present many problems of their own [72]. Within ECMC [42], domain decomposition leads to residual interactions that destroy the global translational symmetry, so that the very convenient separation into non-communicating sectors of lifting variables breaks down, and some of the ECMC efficiency is lost.
In ECMC, one may freely choose the number of active particles (those with non-zero velocities) and keep this number fixed throughout a simulation. If this number is √ N , the mathematical birthday problem shows that the candidate events are usually disjoint for any two active particles, cutting down on the degree of interference between different active particles. Using a framework of local times, this has allowed to prove for that scheduling conflict hard-sphere systems appear with finite probability in the N → ∞ limit for finite run times [58].

Harmonic model
The harmonic model [89] describes spin-wave excitations for N spins on a lattice, with a total potential U = 1 2 i,j (φ i −φ j ) 2 between neighbors i, j , with non-periodic angles φ i ∈ R, which approximate the small elongations |φ − φ j | ≪ 1 in the XY model with its total potential U = − i,j cos (φ i − φ j ). Spin waves are the dominant excitations in the XY model at low temperatures, notably in two dimensions [26]. The harmonic model also provides the quintessence of phonon excitations in particle systems, for small displacements from perfect lattice positions. In the harmonic particle model, each particle interacts harmonically with a fixed set of neighbors so that disclinations and stacking faults, and other excitations, cannot develop.
Besides its role of isolating phonon (resp. spin-wave) excitations of many systems, the harmonic model is of importance for hard-sphere ECMC whose sequence of constraint graphs [42,58], between resamplings, effectively defines a sequence of models with fixed neighborhood.

Physics of the harmonic model
The harmonic model is exactly solved [89]. It has a single phase for all finite temperatures, but the nature of this phase depends on dimension d. The differences in angle for two spins at positions distant by r for the harmonic spin model (and similarly the difference in elongations with respect to the lattice positions for the harmonic particle model) in equilibrium are given by: For a system of size L, the mismatch of two typical spins or particles is thus of the order of O( √ L) in one spatial dimension, it grows as the square root of the logarithm of L in two dimensions, and it remains constant in three dimensions and higher. The harmonic particle model, in two dimensions, features long-range orientational order but only power-law decay of positional order [66]. Only in more than two dimensions does it have long-range orientational and positional order.

ECMC algorithm for the harmonic model
The ECMC algorithm for the harmonic spin model performs in displacement and particle lifting. Each lifted configuration is described by {φ 1 , . . . , φ N } × {i, D}, where i describes the active spin and D ∈ {+, −} the direction of changing the angular variable φ (which is between −∞ and ∞). Because of the invariance of the pair potentials with respect to the absolute spin angles, no liftings between D = + and D = − are required, and the lifted sample spaceΩ again splits into two equivalent copies. It suffices to retain the "+" direction.o te The harmonic model essentially differs from the spin model in that the 2π periodicity is lost. ECMC monotonically increases each of the φ i , at a difference of Metropolis MCMC, which must allow for changes of the angles in both directions. In event-driven harmonic-model ECMC, for a finite number of neighbors per spin, the displacement per event of each active spin is finite. It can be expected that the total displacement of each spin or particle corresponding to the ∆φ in Eq. (41) decorrelates the system. It thus follows that each particle or spin moves between independent samples O(

ECMC for continuous spin models
ECMC readily applies to XY and Heisenberg-type spin models because their spin-spin pair potentials on neighboring sites i and j write as U(φ i − φ j ) with continuous spin angles φ i and φ j . The invariance of the pair potentials with respect to the absolute spin angles then simplifies ECMC in the same way that translational invariance does for particles. For spin models, the algorithm rotates the spin i in a given sense until the factorized Metropolis algorithm calls for a veto by a neighbor j of i. The spin i then stops and j starts to rotate in the same sense, with the cumulative rotation corresponding to the MCMC time.
In the XY model, the two-dimensional fixed-length continuous spins live on a d-dimensional spatial lattice. The pair potential U ij = − cos (φ i − φ j ) favors the alignment of neighboring spins i and j, and in ECMC, the activity (that is, the rotation) passes from one spin to one of its neighbors. ECMC needs no restarts, and its sequence of active spins realizes, in the two-dimensional lattice, an anomalous diffusion process [45]. ECMC on the Heisenberg model, where spins are three-dimensional, can be reduced to this case, with restarts, through rotations in two-dimensional sub-planes in spin-space. In both models, ECMC lowers the relaxation rates with respect to local reversible MCMC. Although more efficient algorithms are available [31,32,92] for these particular models, the comparison of ECMC with reversible MCMC may illustrate what can be achieved through non-reversible MCMC in real-world applications where, as discussed in Section 1, customized a priori move sets remain unavailable.
In the two-dimensional XY model, ECMC autocorrelation functions can be fully explained in terms of spin waves and topological excitations. The two-dimensional hard-disk system similarly features phonons and two types of topological excitations, but the precise relaxation dynamics of MCMC algorithms has yet to be clarified (for a synopsis, see [6, Chap. 1]).

Spin waves and topological excitations in the two-dimensional XY model
Although it possesses no long-range spin order at finite temperatures, the two-dimensional XY model famously undergoes a phase transition [47,31] between a low-temperature phase rigorously described by spin waves with bound vortex-antivortex pairs [89,26], and a high-temperature phase where these topological excitations are free.
In the low-temperature phase, vortices pair up with antivortices. The maximum d max of all vortexantivortex pair separations in a system of size L × L follows a Fréchet distribution [52], with a size-dependent scale s = L 2/α s 0 and a size-independent exponent α that depends on the inverse temperature β. Throughout the low-temperature phase, d max ∼ L 2/α , with α > 2 increases slower than the system size L for non-zero temperatures. Only in the zero-temperature limit does d max remain constant as the system size increases (α → ∞ for T → 0). It is extensive as the transition point is approached (d max ∼ L, that is, α → 2 for T → T c ).

Relaxation time scales in spin models
In spin models, local MCMC simulations must relax spin waves and topological excitations. In the two-dimensional XY model of length L, reversible energy-based MCMC can be expected to diffusively relax spin waves on a O(L 2 ) time scale at all temperatures [52] and topological excitations on a O(d 2 max ) time scale, so that the overall correlation time scales as O(L 2 ) from both contributions. ECMC relaxes spin waves on a much faster time scale that, from the analogy with the harmonic model, can be thought to be O(log L), leading to very fast (system-size-independent) initial decay of spin-spin autocorrelations. The diffusive relaxation of topological excitations by ECMC on the O(d max ) = O(L 4/α ) is dominant at large times throughout the low-temperature phase, with the two contributions leading to a two-timescale decay of spin autocorrelation functions [71,75]. The above arguments are borne out by numerical computations [52] with the α of Eq. (42) as the single nontrivial parameter.

ECMC AND MOLECULAR SIMULATION
In molecular simulation of classical atomic models with explicit solvents, sampling plays an all-important role [85], in addition to the study of explicit time dependence. However, the complexity of total potentials only allows for unbiased local moves, and intricate a priori choices or cluster MCMC algorithms [20,59] have not met with success. The Metropolis or the heatbath algorithms are energy-based: they require the evaluation of the total potential. In the Coulomb case, this comes with prohibitive time consumption for MCMC. In the field, molecular dynamics codes are therefore used exclusively. Great effort has been directed towards the computations of the forces (the derivatives of the total potential), which remains onerous and fraught with approximations. Molecular dynamics invariably implements the Newtonian relaxation dynamics. Several molecular-dynamics codes, developed over the last three decades (e.g. [78]), have found score of real-world applications from biology to medicine, physics, and chemistry. Molecular dynamics is thus highly successful, but it leaves little algorithmic freedom.
ECMC, while still in its infancy, may represent a new starting point for the use of MCMC for real-world applications in molecular simulation. First, to sample the Boltzmann distribution π ∝ exp (−βU) without approximations, the potential U (which may contain a Coulomb component) need not be evaluated. Second, in the field where, within MCMC, complex a priori choices for moves continue to remain out of reach, the tools of factorization, lifting, and thinning may well leverage the development of powerful algorithms. Subsection 6.1, reviews the theoretical aspects of ECMC in the context of molecular simulation while Subsection 6.2 provides an overview of the "JeLLyFysh" Python application, 1 an opensource ECMC implementation for molecular simulation.

Theoretical aspects of ECMC for all-atom models
Three main theoretical challenges stand out in the application of ECMC to molecular simulation. First is the choice of factors, in particular of the Coulomb factors. The different factors are independent and do not compensate each other. Different types of factorizations (especially for the Coulomb potential) therefore impact the total event rate and the MCMC trajectories (see Subsection 6.1.1). The second challenge for ECMC, reviewed in Subsection 6.1.2, is the handling of the Coulomb potential. The statistical independence of factors obviates the need for computing the total potential. However, the factor event rates (the derivatives with respect to the active-particle displacement of the factor potentials) must be provided. One such factor may consist in a single pair of charges together with all their periodic image charges ("merged-image" Coulomb factor) or else in one pair of specific periodic image charges ("separate-image" factor). Both choices are practical, with the latter leading to an infinite number of factors already for a two-charge system [43,22]. For both the merged-image and the separate-image formulation, a single Coulomb factor may collect potentials between a number of atoms, for example all those belonging to a pair of neutral molecules. Such a "dipole-factor" choice corresponds at large distances to the interaction of an atom with a dipole, or even to the interaction of two dipoles, leading to reduced event rates. Finally, the choice of lifting variables (in other words of the set L inΩ = Ω × L), as well as lifting schemes (how to select the elements of L) appear as main fixtures in a largely uncharted territory. As reviewed in Subsection 6.1.3, an ECMC event of the factor that triggers the veto is analogous to a collision with an in-state and an out-state. The out-state is not fixed by a physical dynamics and must only satisfy the global-balance condition, allowing for the setup of in-equivalent lifting schemes for larger factors. The influence of the choice of lifting variables (for example, the choice of directions, individual moves, etc.) on the mixing properties remains poorly understood (see Subsection 6.1.3).

Factors, Coulomb factors
The Coulomb problem illustrates the algorithmic inequivalence of different factorizations. With Coulomb pair factors between each individual pair of charges distant by r (both for merged-image or separate-image formulations) the pair-factor potential is ∼ 1/r, and the pair-factor event rate ∼ 1/r 2 . In a three-dimensional box of side L, at fixed density, so that N ∝ L 3 , a typical event rate between two charges is ∼ 1/L 2 and the total event rate for Coulomb pair factors is O(N/L 2 ) = O(N 1/3 ) [43] (see also [22, eq. (89)]). With this rate, one event corresponds to a single particle moving forward by O(1/N 1/3 ). To advance N charges by a constant displacement (for example 1Å), Coulomb pair factors require a computing effort of O(N 4/3 ). This is clearly borne out by numerical simulations (see [43] for details).
For a dipole factor made up of all the atoms in two charge-neutral molecules with finite dipole moments, the active particle in one of the molecules effectively interacts with a dipole in the other molecule. At a distance r between molecules, this gives a factor potential ∼ 1/r 2 and a factor event rate ∼ 1/r 3 . Integrated over three-dimensional space in a box of sides L, the total event rate scales as O(log N). Advancing N charges by a constant distance requires in this case O(N log N) events, as was observed in N-body simulations [22, Fig. 13].
In real-world applications such as the SPC/Fw water model [91], the interparticle potential contains, beyond the Coulomb term, bond-fluctuation, bending, and Lennard-Jones contributions. For large N, the O(N log N) Coulomb event rate dominates the local contributions, as was observed in simulations (see [22, Fig. 15]). In ECMC, it can be expected that the Coulomb dipole factors have mixing and correlation times that are a factor of N 1/3 / log N smaller than those for Coulomb atomic factors, although this has not been verified explicitly.

Computation of Coulomb event rates
The evaluation of the Coulomb potential and of its derivatives is a major challenge for molecular simulation. Molecular dynamics simultaneously requires the derivatives of the full Coulomb potential for all particles at each time step The dominant Ewald-sum-based solution to this problem consists in discretizing charges onto a fine spatial lattice and in solving Poisson's equation with fast Fourier transformation. The complexity of this algorithm is O (N log N) with a prefactor that depends on the numerical precision (see [25] and [22,Sect. 1A]) The difficulty in evaluating the difference of the Coulomb potential also frustrates energy-based local-move MCMC: Updating a single particle costs at least O(N 1/2 ) with known algorithms [9] so that a sweep of local moves costs a prohibitive O(N 3/2 ).
In ECMC with Coulomb interactions, a naive implementation of the factorized Metropolis algorithm would need to compute O(N) Coulomb factor event rates in order to decide on a possible veto of the consensus process. The thinning implemented in the cell-veto algorithm reduces this complexity to O(1), because in a first part, the position-dependent Coulomb factor event rates are "thickened" into timeindependent and precomputed rates. Specific computations are required only to confirm a provisional veto. Approximate information on the location (and possibly the orientation) of molecules is sufficient in order to compute the ("thickened") bounding potentials and event rates, while the confirmation step accesses the factor event rates to machine precision.
The separate-image and the merged-image formulations both lead to conditionally converging sums, that need to be regularized. In the separated-image formulation of the Coulomb problem, counter-linecharges [43] and their generalizations as, for example, the counter-volume-charges [22, Fig. 5], regularize each individual image charge. Likewise, in the merged-image formulation, the sum over all images with associated compensating line charges is equivalent to the standard tin-foil boundary condition [15] for a single pair of Coulomb charges in a finite periodic box (see [22,Sect. III]).

ECMC liftings and lifting schemes for molecular simulation
As emphasized throughout this review, the factorized Metropolis algorithm attributes each MCMC rejection to a single factor and transforms it, in ECMC, into a lifting. For a pair factor, this lifting is usually determined uniquely (the active and the target particles exchange roles). For a three-particle factor, the probability distribution for the new active particle is uniquely determined, whereas for larger factors, there is a choice among probability distributions, socalled "lifting schemes". For Coulomb dipole factors of two water molecules, each such factor contains six particles, and lifting schemes can differ in their inside-flow (active particles on the same molecule for the in-state as for the out-state) and in their outside-flow. This is exemplified in the "inside-first" and the "outside-first" lifting schemes [22].
For a factor in which the two molecules are separated by a distance r, the intra-molecular lifting rate scales as 1/r 3 for the inside-first lifting scheme whereas the lifting rate towards the other molecule scales as 1/r 4 . Integrated over the whole simulation box of length L, the lifting rate inside the original molecule scales as log L, whereas the outside lifting rate remains constant. For large L, the lifting remains inside the original molecule with probability 1. For the outside-first lifting scheme, in contrast, both the intramolecule and the inter-molecule lifting rates are ∼ 1/r 3 as a function of the distance between molecules, and the total intro-molecule and inter-molecule lifting rates both scale as ∼ log L for large L [22, Table 1]. In addition to the lifting schemes, the choice for the set L and for the choice of selections of its elements, remains poorly understood [79].

Jellyfysh application for ECMC
The open-source "JeLLyFysh" application [35] is the first general-purpose open-source implementation of ECMC. The configuration files of its Version 1.0 realize proof-of-concept of ECMC for interacting particle systems and for the handling of the Coulomb interactions. They also showcase different factorizations and implement various liftings and lifting schemes as well as thinning strategies through bounding potentials and through the cell-veto algorithm. The application provides a platform for future method development, for benchmarking and for production code. The application's architecture (using the mediator design pattern) mirrors the mathematical structure of ECMC. Its event-driven nature is far removed from the time-driven setup of present-day molecular dynamics codes.

Mediators, activators, and event handlers
The "JeLLyFysh" architecture is entirely based on the concept of events. A mediator [27] serves as a central hub for all other elements of the code. After each event, this mediator accesses the active particles (in fact, the "active global state") and then fetches from another element, the activator, socalled event handlers that each treat one factor (depending on the active particle but also on the previous event). For each of these event handlers (that is, roughly, for each factor), the mediator then accesses the in-state for each factor, which allows the event handler to compute a candidate event time. All candidate event times are compared in an element of the program called a scheduler. For the shortest time in this scheduler, the corresponding event handler is then asked to compute an out-state. If confirmed, the out-state is committed to the state of the system. Output may also be generated (see [35,Fig. 5]).
Besides those lifting events imposed by the global-balance condition, the JeLLyFysh application implements a number of pseudo-events. These pseudo-events correspond to sampling, for moving from one cell to another, for changing the direction of motion, etc. and even for the start and the end of the program. It also takes into account the fundamental modularity of ECMC. For example, each event handler, for a factor M, only receives the factor in-state c M in Eq. (14) in order to contribute a candidate event. It computes the out-state only if M is the factor that causes the veto, that is, the true event.

Configuration files, performance tests
The configuration files of JeLLyFysh Version 1.0 construct runs for two charged point masses, for interacting charged dipoles, and for two interacting SPC/Fw water molecules, each with a large choice of factors, and high-precision comparison with reversible MCMC. Four different configuration files treat two interacting SPC/Fw water molecules, with their harmonic oxygen-hydrogen bond interaction on one molecule, a three-body bending potential, a Lennard-Jones oxygen-oxygen potential as well as a Coulomb potential between atoms on different molecules [22, Fig. 12]. The configuration files illustrate the possibilities offered to choose between atomic and molecular factors, and to implement event-handlers that invert potentials to obtain direct event rates, to use bounding potentials, or even to invoke the cell-veto algorithm. Correlation functions are shown to agree in the 0.1% range.
Performance tests of ECMC for large systems of SPC/Fw water molecules will be the object of the next version of the JeLLyFysh application. It features substantially rewritten code, with many of the Python modules rewritten in a faster compiled language. Benchmarks will be provided against single-processor LAMMPS code.

PROSPECTS
This article has reviewed the mathematical and algorithmic foundations of ECMC and analyzed a number of illustrative examples. It has also discussed first applications of the method to key models in statistical physics as well as an ongoing initiative to apply it in the field of molecular simulation. It is now time to resume the principal characteristics of ECMC, to formulate the main working assumptions behind its development, and to define its major challenges.
ECMC reinterpretes three key aspects of reversible Markov-chain Monte Carlo methods. First, the detailed-balance condition is replaced by global balance, in other words the condition of vanishing probability flows is replaced by a steady-state condition, which may allow a swifter exploration of sample space. Second, the energy-driven Metropolis algorithm is replaced by the factorized Metropolis algorithm, which is consensus-driven. Systematically, in ECMC, the consensus must be broken by a single factor, and this usually brings about its formulation as a continuous-time Markov chain. Once the vetoing factor is identified, it is sole responsible for the time evolution into and out of the next event. Third, the trademark rejections of the Metropolis algorithm are reinterpreted as liftings, which enable persistent moves. ECMC has no rejections in the lifted sample space and, under certain symmetry conditions such as global translational invariance, allows one to separate certain lifting sectors from each other and to sample only one such sector.
The ongoing development of ECMC for molecular simulation and related fields builds on two assumptions. The first assumption is that Newtonian dynamics, while being physically realistic, is not the fastest possible relaxation algorithm towards the stationary distribution. The three-decade-long investment in molecular-dynamics codes will thus allow alternative approaches to be ultimately successful. The second assumption is that the elaborate a priori choices that lead to global MCMC moves cannot be adapted to real-world applications in molecular simulation and related fields, because of the sheer complexity of the total potentials. This leaves one with local algorithms but, as argued throughout this review, they must be non-reversible in order to be fast. The family of ECMC algorithms provides a framework for non-reversible local MCMC. It opens up numerous opportunities, from the choice of factors and liftings to the new approach to the Coulomb problem that avoids the computation of the potential and its derivatives.
ECMC enlarges the possibilities of MCMC algorithm and is therefore guaranteed to be useful for specific applications in physics and other disciplines. One promising field of application is polymer science [40,73]. As a general method for molecular simulation, ECMC however faces two critical short-term challenges. The first concerns the benchmarks against molecular dynamics for standard realworld systems. Such benchmarks, with single-processor codes for ECMC and molecular dynamics, will be the object of an upcoming version of the JeLLyFysh application. Their success will depend on the incorporation of appropriate lifting choices into this open-source code. The second challenge concerns the parallelization of ECMC. A road-map for multithreaded ECMC exists at present only for hard-sphere systems. The development of genuine parallel event-driven ECMC for generic potentials constitutes an open research subject. More generally, the mathematical understanding of the mixing and relaxation dynamics of reversible and non-reversible local MCMC remains rudimentary for next to all N-particle models in more than one dimension. A better mathematical grasp of such systems, beyond straightforward benchmarks, appears as a longer-term research goal as well as a prerequisite for the development of faster algorithms that are required for applications in physics and throughout the sciences.

FUNDING
The author acknowledges support from the Alexander von Humboldt Foundation.

ACKNOWLEDGMENTS
The author is indebted to S. C. Kapfer for discussions on Section 3 and to S. Todo for discussions on Section 2.2.2.