Acceleration of Convergence to Equilibrium in Markov Chains by Breaking Detailed Balance

We analyse and interpret the effects of breaking detailed balance on the convergence to equilibrium of conservative interacting particle systems and their hydrodynamic scaling limits. For finite systems of interacting particles, we review existing results showing that irreversible processes converge faster to their steady state than reversible ones. We show how this behaviour appears in the hydrodynamic limit of such processes, as described by macroscopic fluctuation theory, and we provide a quantitative expression for the acceleration of convergence in this setting. We give a geometrical interpretation of this acceleration, in terms of currents that are antisymmetric under time-reversal and orthogonal to the free energy gradient, which act to drive the system away from states where (reversible) gradient-descent dynamics result in slow convergence to equilibrium.


Introduction
In this paper we analyse the effects of breaking detailed balance for interacting particle systems (as described by Markov processes [32]), and their hydrodynamic scaling limits (as described by Macroscopic Fluctuation Theory [9]). The interacting particle systems represent microscopic descriptions of physical systems, in which the motion of each particle may be followed individually. The (fluctuating) hydrodynamic model of the same system describes its behaviour on large length and time scales, in which case the motion of the individual particles is no longer visible, and one works instead with a smooth density field, whose time evolution includes a deterministic element as well as a (weak) stochastic noise [28].
Among interacting particle systems, those with detailed balance are special-they correspond to Markov chains that are reversible with respect to an invariant measure π. Physically, these models are important because their steady states are time-reversal symmetric and lack any persistent currents, so they can be used to describe systems that relax to states of thermal equilibrium. They also have applications outside physics, because given a (possibly non-normalised) measure ν, it is straightforward to design a reversible Markov chain whose invariant measure π is proportional to ν. This construction is at the root of many Markov chain Monte Carlo (MCMC) methods [2,34], in which one typically aims to generate large numbers of uncorrelated samples from a prescribed distribution π. Such methods have widespread applications including Bayesian learning, protein folding and cryptography [14].
In both the physical systems and the MCMC methods, an important question is the rate of convergence to equilibrium of the relevant Markov chains. In MCMC, this rate controls the computational cost required to obtain independent samples from π, which is an important factor in the efficiency of the method. In the physical systems, the question of how fast a system converges to equilibrium controls many physical properties including fluid viscosities, and systems' abilities to respond to changes in external conditions, such as temperature.
Recently, several results have become available which show that for a given invariant measure π, reversible Markov chains have the slowest convergence [10,24,30,35,36,39]. Given that most common MCMC methods are based on such reversible models, and that faster convergence is linked to improved efficiency, this observation offers a route towards the development of new and more efficient methods, some of which are already becoming available [5]. Breaking reversibility can be achieved by an explicit modification of transition rates [36], or by an expansion of the state space (lifting) to incorporate persistence of motion or inertial effects [12,15]. The main physical feature of the resulting irreversible Markov chains is that they (generically) have non-equilibrium steady states characterised by finite entropy production and dissipation of energy. Compared to the equilibrium setting, the nature of fluctuations and convergence to steady state in non-equilibrium systems is much less understood, and is an area of important current activity [3,9,13].
To address these questions, this paper presents several new results. First, we revisit existing results for microscopic models, concentrating in particular on the spectral gap of the generator, and how it is affected when detailed balance is broken. Second, we investigate how breaking detailed balance affects the hydrodynamic limit of the model-in this latter case, convergence to equilibrium is most easily analysed via large deviations of the empirical measure [35,36]. Third, we illustrate our general results by numerical results of a simple interacting particle system-the zero-range process [38]. These numerical results are particularly relevant since the analytical results indicate that breaking detailed balance can never slow down convergence to equilibrium, but they provide rather little insight into how much this convergence can be accelerated, nor how this effect depends on the specific way in which detailed balance is broken. We provide some general remarks and comments in this direction.

Characterisation of Convergence to Steady State
A number of methods are available to analyse the time required for a system to reach its steady state. This section contains a brief review of some of them. For microscopic modelssuch as Markov processes on (finite) discrete spaces and SDEs-we mention some recent work showing how breaking detailed balance can accelerate convergence of systems to their steady states. These results serve as a foundation for our results here, which show how these effects manifest on the macroscopic scale.

Spectral Gap
The first-and most common-method for analysis of convergence to equilibrium is to estimate the spectral gap of the generator of the relevant stochastic process. In general, the eigenvalues {λ i } of the generator are complex numbers, there is a simple eigenvalue λ 0 = 0 and all other eigenvalues have negative real parts. The spectral gap α min is the minimal value of |λ r i | among the non-zero eigenvalues, where λ r i denotes the real part of λ i . Roughly speaking, the physical significance of the spectral gap is that the system converges exponentially fast to its steady state, with a characteristic time scale (1) For stochastic differential equations [24,30] and discrete-space Markov processes [36], it has been shown that irreversible processes generically have smaller time scales τ g , compared to reversible processes with the same invariant measure. We provide further results in this direction in Sect. 2.1.1 below, for the discrete space Markov processes that are relevant for interacting particle systems.

Asymptotic Variance
Another set of methods for the analysis of the convergence to steady state is based on empirical time averages. That is, let X t be the state of the system at time t and let f be an observable quantity (test function) whose value at time t is f (X t ). Then the empirical time average of f is The quantity f (T ) is a random variable which-under suitable conditions related to ergodicity-converges almost surely to the expectation value of f , which we denote by E π ( f ). Moreover the distribution of √ T (f (T )−E π ( f )) converges by the central limit theorem to a normal distribution with variance σ 2 f . The latter is referred to as asymptotic variance or time average variance constant (TAVC) which can be obtained as [35] and [40,Sect. 3.5]. Hence, the variance of f (T ) decays for large times as Var( f (T )) ∼ σ 2 f /T . It is then natural to identify a time scale τ for some rate function I f (which depends on the choice of test function f ). We use the notation in (3) throughout this work as an informal way to state large deviation principles: it means that the log probability that f (T ) takes a value in a small interval containingf can be bounded above and below by quantities related to the rate function I f [23,40]. The rate function achieves its minimal value of zero whenf is equal to E π ( f ), and the second derivative of I f at this minimum is related to σ 2 f . The function I f is a level-1 rate function [23]. A yet more detailed analysis is available by considering not just the large deviations of a single test function f but instead to consider large deviations of the empirical measure. That is, for a Markov chain on a discrete space , define the empirical measurē where δ x,y is a Kronecker delta. The empirical measure at time t is a vectorμ T = (μ T (x)) x∈ . For large enough T , ergodicity implies thatμ T (x) converges almost surely to π(x), and the fluctuations of the measure μ in this limit are described by a large deviation principle at level-2: where I 2 is the rate function [17], which now depends on a vector ν instead of a single real argumentf . Note that the (level-1) rate function I f for any observable f can be obtained by a contraction of this large deviation principle, so the function I 2 contains a great deal of information about the convergence of a system to its steady state. Moreover, as might be expected from the terminology "rate function", the quantity 1/I 2 (μ) has an interpretation as a μ-dependent time scale associated with the decay of an initial measure μ to the invariant measure π. Recent work by Rey-Bellet and Spiliopoulos [35,36] has motivated the analysis of I 2 as a measure of the rate of convergence of processes to their steady states. Their work, and that of Bierkens [10], show that breaking detailed balance accelerates this convergence. Note however that in contrast to the spectral gap-where a single number characterises the rate of convergence of the whole system-the rate function I 2 depends on the measure μ for which it is evaluated; similarly the asymptotic variance σ 2 f depends on the specific observable f . In this sense, the information available from the asymptotic variance and the large deviations is greater than that available from the spectral gap, but this extra information may also make these measures harder to interpret in terms of simple acceleration or slowing down of convergence to equilibrium. Of course, other useful measurements of convergence rates are available, such as mixing times [31], cutoff phenomena (see e.g. [29], where cutoff was recently established for the asymmetric simple exclusion process) and log-Sobolev constants (e.g. [16]), but these are not analysed in this work.

Outline
The remainder of this paper is organised as follows: Sect. 2 includes a theoretical analysis of the effects of breaking detailed balance on convergence to steady states, including both Markov chains (Sect. 2.1) and hydrodynamic limits (Sect. 2.2). Section 3 presents numerical results that illustrate this acceleration in the zero-range process: we provides examples in both one-dimensional and two-dimensional settings. Finally, Sect. 4 contains our conclusions.

Acceleration of the Microscopic Dynamics
In this section, we consider an irreducible Markov jump process on a finite state space which contains n states. In terms of interacting particle systems, this process describes the dynamics of a finite number of particles that move on some finite lattice. The process is defined by rates c(x → y) for states x, y ∈ . The condition of detailed balance (or reversibility) is that for some probability measure π and all x, y then In this case the (unique) invariant measure of the Markov process is π. Let the generator of the Markov process be L. The generator has a representation as an n × n matrix and the reversibility condition (6) corresponds to symmetry of L with respect to the L 2 (π) inner product f, g π = x f (x)g(x)π(x). If detailed balance is broken (nonreversible Markov chain), then L is not symmetric with respect to L 2 (π), but one may always write where L S is symmetric with respect to L 2 (π), while L A is antisymmetric. Moreover, L S is a generator for a reversible stochastic process, whose transition rates may be verified to be where π is the invariant measure of L, which is also the invariant measure of L S . (Recall that the original Markov process is finite and irreducible, which ensures that π(x) > 0 for all x). We also identify the off-diagonal elements of L A as Hence one has from (8) that Note that L and L S both are generators, whereas the operator L A is not a generator of a Markov chain. Alternatively one can think of the decomposition of L in L S and L A as follows: consider the Markov process η t (with t ∈ [−T, T ] for some T > 0) associated to L, with the initial condition distributed by the steady state π. The time reversed processη(t) := η(−t) is also associated to a generator, L * , say. The symmetric part of the generator can be recovered as Given these preliminaries, we can now be precise about the sense in which breaking detailed balance accelerates convergence: in all cases we compare the process L with the corresponding symmetrised process L S . (Equivalently, one may imagine starting from a reversible process L S and breaking detailed balance by adding an extra term L A to the generator.) The processes L and L S both converge to the same invariant measure π-one aims to prove that convergence times such as τ g or 1/I (μ) are smaller for L than for L S .

The Spectral Gap
To illustrate how breaking detailed balance accelerates convergence, we show in Proposition 1 below that breaking detailed balance can only increase the spectral gap, so that the convergence of the irreversible process is characterised by a smaller value of the time τ g . This result has been proven in greater generality in [24,36], but we provide a short proof here, for illustrative purposes.
To this end, consider an initial measure μ 0 , and represent it in terms of an eigendecomposition of L, so that where the α j ∈ C are μ 0 -dependent coefficients, while ν j are complex-valued measures which are left-eigenvectors corresponding to eigenvalues λ j of L. The overbar (e.g.ᾱ) denotes the complex conjugate. Decomposing the non-zero eigenvalues λ j into real and imaginary parts, as λ j = λ r j + iλ i j , the measure at time t is given by where · † denotes a matrix transpose. Note that for real-valued eigenvalues (with λ i j = 0) the term in brackets is equal to 2α j ν j , as in this case also the left (and right) eigenvectors are real-valued.
Moreover, λ r j < 0 for all j, since L is the generator of an irreducible finite Markov process. One sees immediately that this Markov process relaxes exponentially fast to its steady state. Moreover, the rate of this exponential decay is controlled by the non-zero eigenvalue of L whose real part is smallest in magnitude. Similar results to the following proposition have already been obtained in e.g. [26,37]: Proposition 1 Let L and L S be given as above. The non-zero eigenvalues of −L S are real and positive; let the smallest such eigenvalue be α min and the largest be α max . Then every non-zero eigenvalue λ of −L satisfies Proof Define the Dirichlet form for L as E( f, g) where π is the unique stationary distribution of L. Let λ be a non-zero eigenvalue of −L with corresponding right eigenvector f λ + ig λ .
This implies that both f λ and g λ are mean zero, so Var , the bilinearity of the Dirichlet form yields that the real part of λ is given by In addition, one has (for the min and max taken over the two cases h = f λ and h = g λ ) Define E S ( f, g) = f, −L S g π , and note that h,h π . Hence the left hand side of (14) is bounded below by α min . Applying a similar argument to the right hand side of (14) and combining with (13) finally yields (12).

Bounds on Level-2 Rate Functions for Discrete Markov Processes
From Proposition 1 and using (1), one clearly has That is, the irreversible process converges to its steady state at least as quickly as the reversible one. A similar argument [10] establishes that the level-2 rate functions for L and L S are related as again establishing a faster rate of convergence on breaking detailed balance. Recall that results of the form (16) yield information about the empirical measureμ T defined in (4), whereas the previous result (12) concerns the spectral gap and the convergence of μ t , the distribution of the process at time t as defined in (11). Note thatμ T is a random quantity, whereas μ t is the solution to a deterministic differential equation. We now show (Proposition 2) that the rate of convergence of the irreversible model has an upper bound, as well as the lower bound given by I S 2 (μ). That is, I 2 (μ) is bounded both above and below, just as the spectral gap is bounded in (12). This limits the acceleration that is available by breaking detailed balance for (finite) discrete Markov processes, in contrast to the situation for diffusions [35]. The proof for the following proposition is based on the variational formula for the level-2 LDP [23]. Whilst the lower bound, which is known in the literature, see e.g. [10,36], follows from the variational representation of the rate function, the upper bound is (to our knowledge) a novel result. Proposition 2 Consider a finite-state continuous-time Markov chain with generator L = L S + L A and transition rates c(x → y) = c s (x → y) + c a (x → y), as defined in Sect. 2.1. The level-2 rate functional I 2 (μ) is bounded as follows: (17) where the rate functional I S 2 (μ) for the reversible process with generator L S is given by Proof The rate functional is given by a variational formula [23]: In the symmetric case (L = L S ) the maximum is I S 2 (μ), which is attained when f = √ μ/π .
In general we write f = √ μ/π e V for some potential V . A direct computation yields and with If V is a constant function, then I A (μ, V ) = 0 so clearly sup V I A (μ, V ) ≥ 0. Hence, (19) yields the lower bound in (17), as in [10].
For the upper bound, it is convenient to define m(x, y) : π(x)π(y) and q(x, y) := π(x)c(x → y). This yields where we have symmetrised the summand with respect to x, y. For positive constants a, b, one may easily establish the general inequality ae V + be −V ≥ 2 √ ab. Applying this inequality to the summand in (21) yields substituting these results into (22) yields and the combination with (19) establishes the upper bound in (17).

Discussion
Our intuition for the (bounded) acceleration by breaking detailed balance is as follows: for reversible processes we can think of μ t (the distribution of the process at time t) undergoing a steepest descent process (gradient flow) for the free energy , within a particular geometric setting [33]. The precise nature of this geometry is immaterial for this discussion: the key point is that relaxation to equilibrium is fast when the free energy gradient is steep, and tends to be slow when it is shallow. On breaking detailed balance, the free energy still decreases monotonically, but its motion is no longer restricted to the direction of steepest descent. This can have several possible effects and the rate of change of F(t) may either increase or decrease on breaking detailed balance. However, we argue that an important contribution to the acceleration of convergence arises because the irreversible component of the dynamics drives the system away from regions where the free energy gradient is shallow and into regions where it is steeper. We will demonstrate this effect explicitly at the hydrodynamic level, in Sect. 2.2.3.
Notice however, that while slow processes associated with L S are accelerated by breaking detailed balance, the inequality involving α max in Proposition 1 implies that fast aspects of the relaxation tend to be slowed down. Indeed, tr(L A ) = 0 so tr(L) = tr(L S ): since the trace is equal to the sum of the eigenvalues, one sees that if some (slow) processes are accelerated by breaking detailed balance another set of (faster) processes must be slowed down by a similar amount. Within the intuitive picture, our interpretation is that the irreversible component of the dynamics acts to push the system away from regions where the free energy gradient is very steep, so the differences between very fast and very slow processes tend to be smoothed out by the irreversibility.

Accelerating Macroscopic Processes
In this section we consider hydrodynamic limits of interacting particle systems, as described by the macroscopic fluctuation theory (MFT) [9]. We will demonstrate that the large deviation result (16) has a counterpart at the hydrodynamic level. We also explore the geometrical interpretation of this result, and we connect our result to earlier work related to SDEs that describe the motion of single particles [35].

Macroscopic Fluctuation Theory
We first recall the core parts of the macroscopic fluctuation theory (MFT). For a detailed review we refer to [9]. Let ⊆ R d be a connected domain with boundary ∂ . For simplicity, we choose here the domain = [0, 1] d . If we consider a microscopic particle process (indexed by L), its description within MFT involves two random fields, the empirical particle densityρ L t and the empirical current j L t . Roughly speaking, for x ∈ then ρ L t is the local particle density and j L t is a vector that indicates the rate of particle flow. The idea of the hydrodynamic limit is that if we observe an interacting particle system on suitably large scales of length and time, then the system can be described in terms of sufficiently smooth fields ρ and j, instead of requiring a microscopic description in which all particle positions are taken into account. The deterministic quantities ρ and j are then related by a continuity equation given by The domain is fixed in the hydrodynamic limit. The relevance to large length and time scales in the microscopic model is that one considers a large number of particles N within a domain L of linear size L. One takes N , L to infinity together for a fixed densitỹ ρ 0 = N /L d . The domain is obtained by rescaling the (increasingly large) domain L , so that remains fixed as L → ∞.
Within this hydrodynamic limit, the behaviour of the system on suitably large scales of space and time becomes increasingly deterministic. For example, given a time interval [0, T ] and initial and final densities ρ 0 and ρ T , the probability measure for paths connecting these initial and final states concentrates (in the hydrodynamic limit) on a single most likely path. This result can be expressed as a large deviation principle for paths, which can, following [9], be written as with whenever ∂ t ρ t = −∇ · j t is satisfied, and I(ρ, j) = ∞ otherwise. We refer the reader to the review [9] for details on the validity of (24) for a large class of particle systems including the symmetric exclusion process and zero-range processes [28,38].
Note that in contrast to the large deviation principle in Sect. 1.1.3 which is concerned with large times, this principle involves a limit of large L, with a fixed time interval [0, T ].
Physically, we interpret J (ρ t ) in (25) as the most likely current field j t , given that the system has density ρ t . Within MFT, the current is assumed [9, Eq. (2.6)] to have the form where χ(ρ) and D(ρ) are symmetric positive definite d × d matrices that depend on the local density ρ, and E is a fixed (x-dependent) vector field. Physically, D and χ correspond to a density-dependent diffusivity and mobility, while E corresponds to an external force. For a given interacting particle system, the parameters D, χ and E can (in principle) be derived from the microscopic rules of the model. These parameters (along with appropriate boundary conditions associated with ∂ ) fully specify the rate function (25) and they fully describe the hydrodynamic limit of the interacting particle system. To fix the ideas precisely, it may be useful to note that J (ρ) in (26) is itself a field, whose value at position Since J (ρ) is the most likely current for a given density ρ, it follows that for a given initial condition, the path measure is dominated by paths These paths have I = 0 and are said to satisfy the hydrodynamics.
As well as the large-deviation principle for paths (24), the MFT also provides a largedeviation principle for the fluctuations of the instantaneous density, in the steady state of the system. That is, if the time T is large enough that the system has converged to its steady state, one has where V is called the quasipotential: it determines the probability of fluctuations in the density. Eq. (27) is derived under the assumption that the adjoint dynamics satisfy a further Large Deviation principle for a rate functional I * . We refer to chapter II in [9] for a detailed discussion.
We assume throughout that our system has a unique steady state, for which the most likely (x-dependent) density is ρ. In this case V(ρ) = 0 and V(ρ) > 0 for all ρ = ρ.

Reversible and Irreversible Systems
For the microscopic dynamics, we already observed that the detailed balance condition (6) describes an important special case. By starting from this case, the generator was decomposed into two components (7), corresponding to a reversible process and a correction term that captures the irreversibility. At the hydrodynamic level, there is a corresponding decomposition which takes place at the level of the current: one writes The symmetric part of the current is defined [9, Equ. (2.19)] as where δV δρ denotes the functional derivative of the quasipotential introduced in Eq. (27). The antisymmetric part of the current is orthogonal to J S , in the sense that which is sometimes referred to as a Hamilton-Jacobi equation. Note that this is an orthogonality in the space of fields: the presence of the integral implies that the currents J S and J A do not have to be orthogonal at any specific point x. We note that on combining (29) and (30), one has J A (ρ) · ∇ δV δρ dx = 0; integrating by parts and using (28) one sees that which is independent of J A . Hence the quasipotential is non-increasing for paths satisfying the hydrodynamics, and (for any given ρ t ) its time derivative is independent of J A . The special case in which the microscopic model is reversible has two implications for the hydrodynamic limit as described by MFT. First, reversible models lead to J A = 0, so J = J S . Second, assuming that correlations in the particle model occur only on the microscopic scale, the quasipotential within the MFT takes the simple (local) form [9, Eq. (2.25)] where f (ρ) is the free energy per unit volume. (The dependence of f on ρ is fixed by the microscopic model of interest; note also that both ρ and ρ depend in general on the position Hence for reversible microscopic models, the hydrodynamic current obeys In this case consistency with (26) requires The second of these conditions is required within MFT. It is known as the local Einstein relation since it relates the mobility χ to the diffusion constant D. Note that the equations (34) are consistent with the hydrodynamic limit for a large class of particle systems of 'gradient type', see [9, Chap. VIII, Sect. G]. We end this section with a brief comment on the boundary conditions within MFT. If the boundary is associated with coupling of the system to a reservoir at chemical potential λ, the density at the boundary is fixed such that f (ρ) = λ. If particles cannot penetrate the boundaries, one requires D∇ρ = χ E (and j = 0) on ∂ . Paths (or configurations) that do not respect these boundary conditions have I = ∞.

Breaking Detailed Balance Accelerates Convergence
We now state the sense in which breaking detailed balance accelerates convergence of interacting particle systems at the hydrodynamic scale. For the microscopic models, we compared two Markov chains, with the same invariant measure and generators L and L S . At the hydrodynamic scale, we will compare two systems with the same quasipotential (this corresponds to comparing two microscopic models with the same invariant measure). One system is irreversible and has a general J given by (26); the second system is reversible and so J A = 0. In order to ensure a fair comparison, we also assume that the two models have the same mobility χ(ρ): for Markov processes the equivalent condition was that we always compared models with the same L S . Since V and χ are the same for both models, they both have the same symmetric current J S which is given by (29).
For each of these systems, we consider the large deviations of the time-averaged density, following Sect. 1.1.3. Large deviation principles of the form apply in both reversible and irreversible models. This large deviation principle applies on taking the large-T limit after the hydrodynamic limit: one should take L → ∞ before T → ∞. To obtain bounds on I 2 , we introduce the so-called level-2.5 large-deviation principle for the joint fluctuations of the empirical current and empirical measure [4,8]. That is, If we assume that the paths that dominate the level-2.5 LDP are constant in time, the relevant rate function can be obtained from (24) as if ∇ · j = 0, and I 2.5 = ∞ otherwise. The assumption of time-independent paths is equivalent to assuming that no dynamical phase transition takes place [6,11]. Using this assumption, we now calculate a bound (Proposition 3) for the level-2 rate functionals, which is analogous to (16) in the microscopic case.

Proposition 3
Let the level-2.5 rate functional be given by (37) and let I 2 be the level-2 large deviation rate functional obtained from I 2.5 by contraction. We write I rev 2 for this rate functional if the current is symmetric, J = J S , and we write I irrev 2 for the rate functional for the general case J = J S + J A as in (48). Then Remark Note that this result will be strengthened later. We will obtain in equation (51) an exact identity for I irrev 2 as the sum of I rev 2 and a non-negative quantity.
Proof We write I 2 for I irrev 2 . The rate functional at level-2 can be obtained by a contraction of the level-2.5 rate functional, Note that I 2.5 (ρ, j) as given in Eq. (37) is [using (30)] equal to the sum of the following three summands: The summand in the first line coincides with the symmetric rate functional I rev 2.5 (ρ, j) and the second line is the part that corresponds to the anti-symmetric dynamics. Dropping the first summand in the second line (which is non-negative), we obtain An expansion of the square shows that the right hand side is equal to and the last summand vanishes under the assumption that ∇ · j = 0, as by Eq. (29) We obtain with (39) that To establish (39) we now show that the right hand side of this expression coincides with I rev 2 (ρ). Note that again for j such that ∇ · j = 0, by the same argument as in (42), the reversible level-2.5 rate functional is equal to As one would expect for the reversible case, the infimum in (39) is clearly attained for a vanishing current ( j = 0), so that which completes the proof.
Of course, given the acceleration at the microscopic scale, the result (38) that this acceleration is preserved at the hydrodynamic limit may not be surprising. However, we show below that the geometric structure underlying the MFT allows some stronger results for this acceleration to be established.

Splitting the Current
To understand the geometrical origin of (38) in more detail, we now show that as well as the decomposition (28), the antisymmetric current J A has a further decomposition into two parts which are orthogonal to each other, and are both orthogonal to J S . [Here, orthogonality should be understood in the sense of (30). ] We consider the problem with the boundary condition ψ = 0 on ∂ . For any fixed ρ (such that χ(ρ) and J A (ρ) are sufficiently regular) Eq. (45) has a unique strong solution ψ (see for example Theorem 6.24 in [20]). This solution ψ is therefore a functional of ρ we will denote with ψ(ρ). Equation (45) motivates us to decompose J A (ρ) as where J F (ρ) is a new vector field, which is again a functional of ρ. From (45) we see that for all ρ.
We arrive at the following structure for the hydrodynamic current: Of the three terms on the right hand side, the first is familiar as the symmetric current, while the third is divergence free and so does not transport any density. The remaining term (involving ψ) specifies how the density is transported by the antisymmetric current, and also determines the large deviations at level-2. The latter will be established below as a consequence of the following proposition.

Proposition 4
The three terms on the right hand side of Eq. (48) are all orthogonal in the sense of Eq. (30). Moreover, J S (ρ) and −χ(ρ)∇ψ(ρ) are orthogonal to all divergence free vector fields that vanish on the boundary.
Combining Eq. (48) and Eq. (47), the dynamics of the density is given by The first term on the right hand side describes steepest descent (gradient flow) of the quasipotential, within a (modified) Wasserstein metric [1,27]. The second term describes a current that is orthogonal to the gradient flow (within the same metric), and leads to an evolution of ρ within the level sets of the quasipotential: this is the geometric result anticipated in Sect. 2.1.3, but in this hydrodynamic setting the geometrical objects are more explicit. We now derive exact formulas for the level-2.5 and level-2 rate functionals based on the splitting in Proposition 4.

Proposition 5
Let the level-2.5 large deviation rate functional be given by (37). Further let ρ be such that Eq. (45) has a unique classic solution (up to a constant) and j such that ∇ · j = 0. Then, Moreover, the level-2 rate functional is given by Proof The proof of equation (50) follows from Proposition 4 and the representation of the rate functional (40). The second result (51) follows readily as j = J F (ρ) is the minimiser of (50).
Note that these results are consistent with (43) and (44), where the minimising current was given by j = 0. In the general case, the minimising current is given by j = J F (ρ). We moreover can recognise the first term on the right hand side of (51) as I rev 2 (ρ), so the second term on the right hand side is an exact formula for the difference in rate for reversible and irreversible processes. This shows that the convergence rate for the irreversible process is strictly faster, unless the force (−∇ψ) vanishes. We recognise this as a condition that the antisymmetric part of the current contributes to the time derivative of the density (otherwise the convergence to equilibrium of the density can not be accelerated).
Note that the objects ∇ δV δρ and ∇ψ should be interpreted as forces acting in the space of densities. In order to sustain a large deviation of the density, the stochastic forces within the system must act to resist these (deterministic) forces. One sees from (51) that the probability of this rare event (or large deviation) is given by the norms of the two forces, within a metric that depends on the mobility χ.

An Example
We have discussed the status of the MFT as a theory for the hydrodynamic limit of interacting particle systems. For a concrete example of this approach, we consider an interacting particle model known as the zero-range process (ZRP) [38]. A microscopic description of the ZRP is given in Sect. 3.1. For the purposes of this section, the important features of the ZRP are that its hydrodynamic limit is described by the MFT and that irreversible ZRPs have local quasipotentials of the form (32). This latter fact allows straightforward comparison between reversible and irreversible models with the same quasipotential.
The hydrodynamic limit of the ZRP is a non-linear drift-diffusion where φ is a function that depends on the local density [that is, φ(ρ)(x) = φ(ρ(x))], and E is a drift term. The specific function φ that appears in the MFT depends on how the particles interact within the ZRP. A formal derivation of this hydrodynamic limit can e.g. be found in [9]. If φ(ρ) = ρ, then the model corresponds to drift-diffusion of non-interacting particles. One sees immediately from (52) that the hydrodynamic current is given by (26) with χ(ρ) = φ(ρ)I and D(ρ) = φ (ρ)I , where I is the identity matrix. Moreover, the quasipotential for the ZRP is given by (32) with f (ρ) = log φ(ρ), consistent with (34). The ZRP may be either reversible or irreversible: one sees that reversible ZRPs lead to E = −∇V for some potential V . In this case (34) shows that V (x) = log(φ(ρ(x))) + λ, where ρ is the steady state density profile and λ is a constant (independent of x). Hence one identifies the irreversible current as Examining the rate function (51) for the specific case of the ZRP, one can interpret the result as a generalisation of a result in [35]. One has δV/δρ = log φ(ρ) − log φ(ρ). Hence where ψ is the solution of ∇ · φ(ρ)∇ψ = −∇ · [φ(ρ)(E + ∇ log φ(ρ))]. If we now consider the special case φ(ρ) = ρ then we recover the same rate function as in Theorem 2.2 of Ref. [35]: the non-gradient force C in that work is here replaced by E + ∇ log ρ (note that this is independent of ρ). The condition that ∇ · (ρC) = 0-which ensures that the invariant measure is unchanged by breaking detailed balance-is satisfied within the MFT because ∇ · J A (ρ) = 0 and setting φ(ρ) = ρ yields J A (ρ) = ρ(E + ∇ log ρ). Note however the setting discussed in this work is different to that in [35]: here we consider the hydrodynamic limit of many particles on a lattice while that work considers a single particle in a compact manifold without boundary. For non-interacting particles, the result is the same: the reason that for the many-particle system, the rate function I N associated with all the particles undergoing the same rare fluctuation is equal to N I 1 . So the only difference between the one-particle and many-particle systems arises in the prefactors (speeds) of the large deviation principles (35), (36).

The Zero-Range Process
The ZRP [38] is a system in which interacting particles move on a finite lattice L = {0, . . . , L − 1} d ⊆ Z d where L ∈ N is the linear system size. The particles are assumed to be indistinguishable and each particle is located at one of the sites x ∈ L . The number of particles on site x is η(x) and the configurations of the system are η = (η(x)) x∈ L . We will assume that the total number of particles is conserved such that no particles are added or removed over time.
The interaction of the particles is encoded in a function g(k), with g(0) = 0. The rate of particle transfer from site x to site y is g(η(x))c(x → y), where the function c determines the connectivity of the sites. The case g(k) = k corresponds to non-interacting particles. The model is referred to as zero-range because particles interact only when they are on the same site. For example, if g(k) = k α for k > 0, then α < 1 means particles on the same site attract each other (suppressing jumps away from that site) while α > 1 means that particles on the same site tend to repel each other.

Reversible and Irreversible ZRP
The behaviour of the ZRP depends strongly on the choice of the connectivity function c as well as the interaction function g. We assume that particles hop only to nearest neighbour sites, so c(x → y) > 0 only if x and y are nearest neighbours. At the boundaries of the lattice, the system has either reflecting boundaries (particles cannot leave the lattice) or periodic boundaries.
It is easily verified that the model obeys the detailed balance condition (6) if one takes (for nearest neighbour sites) for some potential function V . In this case the model is reversible.
To arrive at a class of irreversible models, we take with k x,y = −k y,x . In this case positivity of transition rates requires |k x,y | < e − 1

Generator and Invariant Measure
We denote the configuration of the ZRP at time t with η t . The generator acts on the test function f as Here η x,y denotes the configuration obtained from η by removing one particle from position x and adding it at position y. If η(x) = 0 we simply set η x,y = η and hence leave the configuration unchanged. Note that the ZRP as defined so far is reducible, since the number of particles is a conserved quantity under the dynamics. This setting is useful because it is easily verified (directly from the definition (56) and using that the invariant measure π satisfies η π(η)L f (η) = 0 for all f ) that the reversible model with rates defined in (54) has a family of invariant measures, the so called grand-canonical measures, which are parameterised by the chemical potential λ and given by with the fugacity ϕ(x) = e −V (x)−λ for some λ ∈ R; the notation g!(k) indicates the generalised factorial g!(k) is a normalisation constant [19,28]. We here assume that V , λ and g(·) are such that z(ϕ(x)) < ∞ for all x ∈ L . This is in particular the case for any V and λ, when g(·) satisfies g(k) ≥ ck for some constant c > 0 [28]. On restricting the model to a fixed number of particles N , the invariant measure π (which is called the canonical measure) can be obtained by a conditioning of (57). Note that (57) has the structure of a product measure. Also if g(k) = k then one recovers the case of non-interacting particles and the local marginals of (57) are Poisson distributions.
To make the comparison between reversible and irreversible models described in Sect. 2.1, we require an irreversible model whose invariant measure is (57). Again using that η π(η)L f (η) = 0 for all f , we take f = η(x) to be the number of particles on site x, from which we see that the irreversible rates (55) are also consistent with the invariant measure (57) if we take y:y∼x where the notation y ∼ x indicates that sites x and y are nearest neighbours. (If we imagine a system with just one particle, this constraint states that the rate of hopping onto site x is balanced by the rate of hopping away from that site. For the ZRP, this same balance condition ensures that the invariant measure (57) is still valid even for many interacting particles). Finally then, the conditions on the perturbations k x,y required for a meaningful comparison between reversible and irreversible models can be summarised as: y:y∼x k x,y = 0, k x,y = −k y,x , and |k x,y | < e −(V (x)+V (y))/2 .
The rates k x,y can be interpreted as elements of a matrix, which coincides (up to the factor 1/2) with the vorticity matrix introduced in [10].
In terms of the splitting (7) the symmetric part of the dynamics is given by c s (x → y) = e 1 2 [V (x)−V (y)] and the anti-symmetric part by c a (x → y) = k x,y e V (x) , such that the symmetric part (corresponding to L S ) is independent of k x,y .

Hydrodynamic Limit
The hydrodynamic limit of the ZRP is defined as follows. For a ZRP on a lattice L with L d sites, one takes N = ρ 0 L d particles, where ρ 0 is the average density. The lattice L is mapped into the domain [0, 1] d by identifying each site x ∈ L with a positioñ x ∈ with = [0, 1] d . Hence the site x with integer co- ordinates (i, j, . . . ) has a positioñ x = (i/L , j/L , . . . ). Roughly speaking, the density ρ t (x) in the MFT is equal to the typical number of particles on site x, and the normalisation of the density is ρ t (x) dx = ρ 0 . The hydrodynamic limit corresponds to a sequence of models in which L → ∞ at fixed ρ 0 , so N → ∞.
The hydrodynamic limit corresponds to observing a system on increasingly large length and time scales. Note that since the number of sites in L is diverging (proportional to L d ) in the hydrodynamic limit, the diffusion constant for a single particle (in ) vanishes as L −2 . For this reason, when the lattice L is mapped into the fixed domain , it is also convenient to scale the hop rates for all particles, by taking c(x → y) → L 2 c(x → y). This ensures that the diffusive behaviour characteristic of the hydrodynamic limit is observed, and the hydrodynamic limit is consistent with MFT.
To fix the hop rates between sites in the ZRP, one fixes a smooth potential functioñ V : → R on the hydrodynamic scale, and one considers a sequence of ZRPs of increasing sizes L with potential functions V (x) =Ṽ (x), wherex is the image in of the discrete site x ∈ L . Similarly one fixes a vector fieldk : → R d and takes k x,y =k(x) · (ỹ −x) where the dot indicates a scalar product in R d .
The relation between the ZRP and the MFT is discussed in e.g. [7], [22] and in the review paper [9]. In particular, for both reversible and irreversible ZRPs one arrives at the situation described in Sect. 2.2.5. The hydrodynamic limit (52) depends on the drift function E : → R d which is given by The MFT description of the ZRP also depends on a function φ which can be obtained as the solution of We identify the right hand side of this equation as the mean local density associated with the measure (57), at fugacity ϕ = φ(ρ). The quasipotential V for the ZRP is given by [9],

Simulation Results
We present numerical results for one-dimensional and two-dimensional systems, showing how breaking detailed balance [that is, taking k x,y = 0 in (55)] accelerates convergence to equilibrium. The simulations are performed using the Gillespie algorithm [21]. The results illustrate several aspects of the theoretical analysis in Sect. 2. First, the results of that section do not rely on how detailed balance is broken: we show that there are several possible choices and discuss their consequences. Second, our numerical results show in what contexts we expect to see significant acceleration of the dynamics on breaking detailed balance, and in what contexts we expect the acceleration to be mild. In all cases, we show results that are scaled to be consistent with the hydrodynamic limit. That is, we map the lattice L into [0, 1] d and we rescale the microscopic hop rates by a factor of L 2 so as to recover diffusive behaviour in the hydrodynamic limit.
In practical situations where the rate of convergence to equilibrium is important, a common situation is that the potential function V is not convex, but includes several (or many) minima, separated by high barriers. From a physical perspective, the temperature of our systems is a parameter that has been absorbed into the function V . In general, high barriers are linked with long (Arrhenius) time scales that are proportional to e V . In order to understand whether breaking detailed balance can accelerate convergence in such non-convex problems, we consider cases where the function V has two minima, with longest time scale in the system corresponding to motion between these minima.

Characterisation of Convergence
We perform numerical simulations starting from a fixed (deterministic) initial condition η 0 . To analyse convergence to equilibrium, we perform numerical simulations of the ZRP, and we track the time-dependence of several different quantities. For any configuration η, the mean potential energy is We generate several trajectories (sample paths) η t of the ZRP and we estimate the mean potential energyV by taking the mean value of η t , V over these trajectories. For systems of non-interacting particles (where φ(ρ) = ρ), we also estimate the macroscopic relative entropy as which can be seen as an approximation to the quasipotential, which is for an independent random walk given by where we used the fact that z(ϕ) = e −ϕ in (61) and the last identity follows from the fact that the density is conserved: For numerical purposes, we estimate E μ 0 (η t (x)) as the average occupancy of site x over the sample paths that we generate, and we calculate E π (η(x)) by direct construction of the invariant measure (whenever possible). Finally, we estimate the Gibbs entropy which is large if particles are delocalised throughout the system, and small if they are concentrated on a small number of sites. Again, we estimate E μ 0 (η t (x)) as the average occupancy of site x over the sample paths that we generate, which provides an estimator of S. These three quantitiesV , D, S all converge as a function of time to stationary values, providing differing information as to the rates of convergence. Note that for non-interacting particles, ρ(x) = E π (η(x)) = e −V (x) /z for some constant z, so D(t) = −S(t)+V (t)+log z.

One-Dimensional Case: Results
We consider periodic boundaries for a model on a one-dimensional strip, this is equivalent to motion on the perimeter of a circle (flat torus in one dimension). In this case condition (59) requires k x,x+1 = k x−1,x , so we set k x,x+1 = c with some constant c that is independent of x. The choice c > 0 corresponds to a fixed force c e V that is forcing the particles to travel around the circle. For a hydrodynamic limit consistent with macroscopic fluctuation theory, we require c to vary with the system size L as c = E/L with E a fixed constant [9].
We note in passing that the use of periodic boundaries is essential for breaking balance in these closed systems: on a finite strip with reflecting boundary conditions, (59) has no solutions except k x,y = 0 so there is no way to break detailed balance.
Thus, returning to the case with the periodic boundaries, the generator is where the addition is periodically extended on L = {0, . . . , L − 1}, i.e., (L − 1) + 1 = 0 and 0 − 1 = L − 1. We take g(k) = k so that the particles do not interact. The potential is with A = 3/2 and B = 3/4 so that the global minimum of the potential is atx ≈ 0.888 with V ≈ −2.052. The height of the barrier is approx 2.609. The initial condition has all particles on a single site, x 0 = L/4, in the vicinity of the secondary minimum. The stationary state has ρ(x) = E π (η(x)) ∝ e −V (x) with a proportionality constant determined by the total density (which in this case is z ≈ 2 377). The parameter E in Eq. (66) is set to E = 36. For the lattice size L = 300, the maximal value allowed for E to ensure that c s + c a ≥ 0 is slightly above 38.4. In principle one can choose larger values for E by increasing the lattice size L.
The results in Fig. 1 are for a domain of size L = 300; we also compared this to simulations for L = 150, L = 300 and L = 450 for the value E = 18 (to ensure positiveness of the transition rates for L = 150). We found the results to be qualitatively very similar, see the bottom right panel in Fig. 1. Figure 1 shows the convergence to equilibrium of the mean potential energy and the entropy. One sees that convergence of both the energy and the entropy is significantly faster when detailed balance is broken. To illustrate the mechanism for this effect, Fig. 2 shows how the mean density E μ 0 (η t (x)) varies with time.
In the irreversible case, the non-gradient part of the drift force E acts to the right and is equal to c e V , so it is large near the maxima of the potential. This prevents the system from becoming localised in the secondary (local) minimum and aids convergence to the steady state. By contrast, in the reversible system, the particles need to diffuse over the maxima of the potential, which is a slower process. This difference explains the much faster convergence  to the steady state observed in Fig. 1. The overshoot of the entropy for the reversible case in Fig. 1 occurs because the state where the particles are distributed evenly between the two minima has a higher entropy S than the steady state (where they are localised primarily in the global minimum). The state where the particles are distributed evenly between the minima is an example of a situation where the gradient of the free energy is small (within the relevant metric), so that steepest descent of the free energy leads to slow changes in the density. Note also that (64) implies that D(t) → 0 at long times, as the system converges to its steady state. However, in Fig. 1 one sees that our estimate of D(t) converges instead to a small positive constant. This offset arises because our estimator of D(t) is biased: it is based on m independent numerical simulations (each with N particles) and the expectation value of our estimator converges to D(t) only as m → ∞. Specifically, we estimate E μ 0 (η t (x)) as is the number of particles on site x at time t in the k-th simulation. Inserting this estimate into the (nonlinear) expression (64), it is easily shown that the resulting estimator of D(t) has in general a finite bias. However, as m → ∞, ϑ obeys a law of large numbers and converges almost surely to E μ 0 (η t (x))-hence our estimator converges to D(t) as m → ∞.

One-Dimensional Case: Discussion
This one-dimensional model is useful for illustrative purposes and establishes the general principles derived in Sect. 2. However, the restriction to one dimension means that detailed balance can only be broken by applying a driving force c e V (otherwise the invariant measure would be changed). If barriers are large, one sees that the driving force near the top of the barrier must be very large indeed: it is hard to see how this can be realised in practical applications. Physically, the idea is to drive a constant current around the periodic system, and this requires the drift velocities (and hence forces) to be largest at the top of any barriers, where the density is least. In this sense, it is perhaps not surprising that by applying large forces to quickly drive particles over all barriers in the system, one can significantly speed up mixing of the particles between the two minima of the potential.
For these reasons, we turn to a two-dimensional system, where there are many more ways of breaking detailed balance while preserving the same invariant measure.

Two Dimensional Case: Model and Results
In two dimensions, there is considerably more freedom in the choice of the rates k x,y . If one again assumes periodic boundaries, it is always possible to have all non-gradient forces acting in a single direction: for example k x,x+e 1 = c where e 1 is a lattice vector, as in the previous one-dimensional example. However, this requires driving forces that depend exponentially on the value of the potential, as in one dimension. We therefore pursue a different strategy.
Denoting the Euclidean basis for L with e 1 , e 2 , Eq. (59) implies that both k x,x±e j = −k x±e j ,x and k x,x+e 1 + k x,x−e 2 + k x,x−e 1 + k x,x+e 2 = 0 have to be satisfied. One way to choose appropriate k x,y is to consider the plaquettes of the square lattice as in Fig. 3 and to define a vorticity W at the centre of each plaquette. The value of W on the plaquette centred at x + 1 2 (e 1 + e 2 ) is W (x). One then can choose the rates k x,y as the following differences:  This choice satisfies both conditions k x,y = −k y,x and y k x,y = 0. The quantity W can be identified as a vorticity, in the sense that taking W (x) = W 0 δ x 0 ,x with W 0 > 0 causes particles to circulate clockwise around plaquette x 0 . Any choice of the function W is possible, and should lead to acceleration of the dynamics, following the theoretical analysis of Sect. 2. Here we concentrate on a case where W is related to the potential V , so that the rates c(x → y) depend only on the gradients of the potential in the vicinity of site x. (The physical idea is that particle motion is naturally sensitive to local potential gradients since these correspond to forces acting on the particles. On the other hand, the motion of a particular particle should not be sensitive to the total energy V , since this depends on the state of the system far away from that particle). To arrive at forces that depend only on potential gradients, we take , where a is a parameter that sets the scale of the vorticity.
On taking the hydrodynamic limit, this gives rise to the driving force where a > 0 (recall from Sect. 3.1.3 thatx is the image in of the discrete site x ∈ L ). We recognise the second term on the right hand side as a force that is obtained by rotating ∇V clockwise by π/2 radians, so that it acts to drive the system around the level sets of V . The following simulations are on a two dimensional closed domain with L = 140 and zero flux at the boundary, i.e., the domain has 140 × 140 = 19 600 sites and the particles cannot leave the domain.
We consider three different ZRPs, corresponding to different choices for g(k). Firstly, we consider the linear case (independent particles), where g(k) = k. We further consider the superlinear case with g(k) = k 3/2 , such that the particles repel each other (the hop rates away from site x is increased when that site contains more particles). Finally we investigate the sublinear case with g(k) = k 5/6 in which the particles prefer to cluster together. For each setting, we simulated the process with both reversible and irreversible dynamics, with L 2 /2 = 9 800 particles averaged over 16 simulations. The potential, which is also depicted in Fig. 4, is for shifted coordinates x = (x 1 , x 2 ) ∈ [−1/2, 1/2] 2 given by with a cut-off at a given height V * . For the simulations we chose the parameters A = 500, B = 0.085, C = 30, D = 2.5 and V * = 5 (that is, the potential used is max (V (x 1 , x 2 ), V * )).
The parameter in (69), which sets the strength of the non-gradient term of the driving force, rev., t=0.1 x 1 x 2 irrev., t=0.1 x 1 x 2 rev., t=0.9 x 1 x 2 irrev., t=0.9 x 1 x 2 rev., t=1.7 x 1 x 2 irrev., t=1.7 x 1 x 2 rev., t=2.5 x 1 x 2 irrev., t=2.5 was set to a = 0.4. This value is again close to the maximal allowed value (which is slightly above 0.405). For all simulations, the particles start at position (0.5, 0.75) ∈ [0, 1] 2 close to the local minimum of the double well potential. The particles then try to leave this well and move to the global minimum (on the left) as can be seen in the plots in Fig. 5 for the linear case. The test observables for the linear/superlinear/sublinear case can be found in Figs. 6, 7 and 8, respectively. Depending on the chosen configuration, the simulation time on a HPC node with 16 cores using Matlab took between 10 and 13.5 hours.
As in the one-dimensional case, the particles are under the irreversible dynamics able to leave the minimum faster than it is the case for reversible dynamics (compare the bottom row with the top row in Fig. 5).

Two Dimensional Case: Discussion
We close this section with Table 1, which quantifies the acceleration in the models where particles attract, repel, or have no interactions. For this, we consider the average energyV and the average x 1 position of the particles. Assuming that the final values of these observables in the irreversible simulations are close to their steady-state values, we consider the distance  V (resp. x 1 ) of both the reversible and irreversible process and keep track of the first time where the distance is below a given threshold. Denoting this time for the reversible process with t s and for the irreversible process with t a , we can use the ratio t s /t a as an estimator for the acceleration. From the data in the table, on sees that the processes are typically accelerated by factors about 1.75 independent of the choice of g(k). We checked different thresholds (here we displayed V = 0.3 and x 1 = 0.2) which all lead to the same conclusions.
These are significant accelerations, although considerably less than the dramatic speedup of order 10 observed in one dimension. However, the physical mechanisms for the acceleration are different in the two cases. In one dimension, the drift forces which act to push particles up and over the barrier, so the forces are very large at the top of the barrier. In two dimensions, the effect is more subtle: returning to Fig. 4 and recalling that the drift force in (69) is obtained by a rotation of the potential gradient, one sees that in the vicinity of the saddle point of the potential, there is a net drift to the left in the top part of Fig. 4b, and a drift to the right in the bottom part. A natural analogy is a gentle stirring motion that happens in the vicinity of the saddle point, and tends to accelerate mixing. This seems a much more plausible mechanism for accelerating convergence to equilibrium in practical situations, compared with the large forces required in one dimension.
Finally, we note that transport between the minima of a non-convex potential energy always involves a slow time scale proportional to e V , since a particle must still reach the barrier in order to cross it, and the probability that a particle visits the barrier is proportional to e − V . However, the results here show that mixing of particles between energy minima can be accelerated by enhancing the probability that if a particle reaches a region with high V , it takes advantage of this excursion in order to cross the barrier. The mechanisms for this enhanced probability differ between the models considered here-it would be interesting to investigate this effect further, so as to understand how general these mechanisms are and how they can be exploited in practical applications.

Conclusion and Outlook
We have considered interacting particle systems described by Markov chains, and their hydrodynamic limits, as described by macroscopic fluctuation theory. We compare reversible and irreversible processes: for an irreversible system with generator L, the corresponding reversible process is the one identified in (7), whose generator is L S . At the microscopic level, it is known that the irreversible process then converges to its steady state at least as fast as the reversible one-this can be demonstrated by considering either the spectral gap or the (level-2) large deviations of the empirical measure. In the hydrodynamic limit, Eq. (38) shows that this property is preserved, by considering the large deviations of the empirical density. Moreover, Eq. (51) gives a quantitative expression for the acceleration of convergence, which may be seen as a generalisation of previous results for single-particle diffusions [35].
Our numerical results for the ZRP reinforce the observation that for a given reversible system, there is a large family of irreversible systems for which convergence to equilibrium is faster (or, at least, equally fast). We considered two cases: either a drift force in a single direction, which acts to drive a system around a circle (Sect. 3.2.2) or the introduction of a force that drives the system around the level sets of the potential (Sect. 3.2.4). In both cases, we observe acceleration of convergence, as expected.
The results within MFT provide a geometrical interpretation of the acceleration, in terms of forces that act in directions perpendicular to the free energy gradient, as shown by orthogonality relations for currents such as Eq. (46). We have argued that such forces can act to accelerate convergence by driving the system away from regions where the free energy gradient is shallow, in which cases reversible processes exhibit slow convergence.
We offer two perspectives on future application of these ideas. First, we have shown that breaking detailed balance generically accelerates convergence, but of course there are very many ways to write down irreversible models, and it is not clear what choices are most practical in applications, nor which ones lead to the fastest convergence. In particular, the choice considered for ZRP examples shown here are rather specific to systems in one or two dimensions. (We emphasise however that the configuration spaces of the ZRP are very high-dimensional since we consider N interacting particles, so the methods are not restricted to systems with low-dimensional configuration spaces.) Second, we gave a geometrical interpretation in which the symmetric dynamics correspond to the gradient flow (steepest descent) of the free energy and the antisymmetric dynamics are in some sense orthogonal to this gradient flow. This offers a potentially new perspective on hydrodynamic limits in irreversible systems, which it would be interesting to investigate further, for example with a view towards obtaining analytic estimates for the rate of convergence.
Supporting data for this manuscript and the code used for the simulations will be made available short after publication on the University of Bath data archive (DOI:10.15125/ BATH-00365).
(HPC) Service at the University of Bath. The authors thank the anonymous referees for very helpful comments and suggestions.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.