Large deviations and optimal control forces for hard particles in one dimension

TL;DR

This paper investigates large deviations in one-dimensional diffusing hard particle systems, analyzing effective forces for biased trajectory sampling and identifying distinct regimes with different interaction characteristics.

Contribution

It introduces an optimal-control framework to study effective interactions and sampling improvements in large deviation regimes of one-dimensional hard particle systems.

Findings

High activity leads to hyperuniform, long-ranged repulsive interactions.

Low activity exhibits near-equilibrium behavior with long-ranged attractive forces.

Far-from-equilibrium regime features macroscopic gaps and strongly size-dependent interactions.

Abstract

We analyse large deviations of the dynamical activity in one-dimensional systems of diffusing hard particles. Using an optimal-control representation of the large-deviation problem, we analyse effective interaction forces which can be added to the system, to aid sampling of biased ensembles of trajectories. We find several distinct regimes, as a function of the activity and the system size: we present approximate analytical calculations that characterise the effective interactions in several of these regimes. For high activity the system is hyperuniform and the interactions are long-ranged and repulsive. For low activity, there is a near-equilibrium regime described by macroscopic fluctuation theory, characterised by long-ranged attractive forces. There is also a far-from-equilibrium regime in which one of the interparticle gaps becomes macroscopic and the interactions depend strongly on the size of this gap. We discuss the extent to which transition path sampling of these ensembles is improved by adding suitable control forces.

Equations

$$p(X)=\frac{1}{Z}\exp [-\beta U(X)]$$

$$p_{\mathrm{acc}}=\frac{1}{1+\exp (\beta \Delta U)}$$

$$v(x)=\cases{0,&$x>l_{0}$\\ \infty,&$x<l_{0}$}$$

$$r_j^M=\frac{1}{2M}\left[\min (M,|x_j-x_p|)+\min (M,|x_j-x_q|)\right]$$

$${\dot{x}}_i=-\beta D_0{\nabla }_{i}U+\sqrt{2D_0}{\eta }_i$$

$$L_{\mathrm{r}}=L-N{l}_0.$$

$${\Phi }_a=\frac{Na}{L_{\mathrm{r}}}$$

$$K[\mathbf{x}]=\sum _{i=0}^N{\int }_0^{t_{\mathrm{obs}}}{r}_i^a(t)\mathrm{d}t$$

$$k[\mathbf{x}]=\frac{K[\mathbf{x}]}{N{t}_{\mathrm{obs}}},$$

$$⟨r_i^a{⟩}_0={\int }_0^a(1/a){\int }_z^{\infty }(N/L_{\mathrm{r}}){\mathrm{e}}^{-yN/L_{\mathrm{r}}}\mathrm{d}y\mathrm{d}z$$

$$⟨k[\text{x}]{⟩}_0=\frac{1}{{\Phi }_a}(1-e^{-{\Phi }_a})$$

$$p(k|t_{\mathrm{obs}},N)\sim {\mathrm{e}}^{-t_{\mathrm{obs}}I(k)}$$

$$⟨O{⟩}_0=\int O[\mathbf{x}]\mathrm{d}{P}_0[\mathbf{x}]$$

$$\mathrm{d}{P}_s[\mathbf{x}]=\frac{e^{-sK[\mathbf{x}]}}{Z_s}\mathrm{d}{P}_0[\mathbf{x}],$$

$$⟨O{⟩}_s=\int O[\mathbf{x}]\mathrm{d}{P}_s[\mathbf{x}].$$

$$\psi (s)=\underset{t_{\mathrm{obs}}\to \infty }{\lim }\frac{1}{t_{\mathrm{obs}}}\log ⟨{\mathrm{e}}^{-sK[\mathbf{x}]}{⟩}_0.$$

$$k(s)=⟨k[\mathbf{x}]{⟩}_s$$

$$I\left(k(s)\right)=-sk(s)-\psi (s).$$

$$P_{s,\mathrm{end}}(z)={⟨\delta [z-\hat{z}(t_{\mathrm{obs}})]⟩}_s.$$

$$P_{s,\mathrm{ave}}(z)=\frac{1}{t_{\mathrm{obs}}}{\int }_0^{t_{\mathrm{obs}}}{⟨\delta [z-\hat{z}(t)]⟩}_s\mathrm{d}t.$$

$${\dot{x}}_i=-D_0{\nabla }_i(\beta U+V_{\mathrm{opt}})+\sqrt{2D_0}{\eta }_i$$

$$\mathrm{d}\tilde{P}^{V}_{s}[\mathbf{x}]\propto\exp\left(\frac{1}{2}\Big{[}V(X(0))-V(X(t_{\rm obs}))\Big{]}\right)\,\mathrm{d}{P}_{s}[\mathbf{x}]\;.$$

$$\mathrm{d}{\tilde{P}}_s^V[\mathbf{x}]\propto \exp \left({\mathcal{A}}^{\mathrm{sym}}[\mathbf{x}]-sK[\mathbf{x}]\right)\mathrm{d}{\tilde{P}}^V[\mathbf{x}],$$

$${\mathcal{A}}^{\mathrm{sym}}[\text{x}]=\frac{1}{4}\sum _i{\int }_0^t{}_{\mathrm{obs}}{\nabla }_{i}V\cdot D_0({\nabla }_{i}V+2\beta {\nabla }_{i}U)-2D_0{\nabla }_i^{2}V\mathrm{d}t.$$

$${\tau }_{\mathrm{B}}=l_0^2/(2D_0).$$

$$L_r=N/\overline{\rho },t_{\mathrm{obs}}={\gamma }_{\mathrm{obs}}\frac{L_{\mathrm{r}}^2}{2D_0}.$$

$$\lambda =s{L}_{\mathrm{r}}^2/D_0$$

$$\dot{\rho }$$

$$j$$

$$(\partial /\partial \hat{t})\rho$$

$$\hat{ȷ}$$

Full text

1 Introduction

Large deviations of time-averaged quantities are becoming increasingly useful for understanding dynamical fluctuations in physical systems [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. For example, consider an ergodic system in which time-averaged quantities converge almost surely to ensemble-averaged values. Given some large time scale, the probability of a significant deviation between the time-average and the ensemble average is small but finite – these rare events are described by large-deviation theory [14, 15]. Despite their scarcity, analysis of these events has led to new insight into the behaviour of physical systems, and their dominant fluctuation mechanisms[16, 10, 17, 9, 11, 18].

Early studies of these events focused on the entropy production in non-equilibrium systems, which is intrinsically linked to fluctuation theorems [1, 19, 2, 20]. Another direction has been the analysis of time-averaged currents, aiming towards a general theory of transport in non-equilibrium systems [3, 5, 21, 22]. Yet another line of enquiry has focused on glassy systems [8, 23, 9, 24], which have long-lived metastable states that hinder equilibration.

In studies of large deviations, there are numerous examples of dynamical phase transitions [17, 8, 23, 9, 24, 25, 13, 26, 27, 28, 29, 10]. In simple terms, these occur when deviations from ergodic behaviour occur by mechanisms that differ qualitatively from the typical behaviour of the model. For example, the rare events may involve spontaneous symmetry breaking, as in [13, 10, 29]. In other cases, one encounters the phenomenology of first-order phase transitions, including phase coexistence[8, 23].

Here we focus on a simple Brownian hard particle model (BHPM), which has rich fluctuation behaviour, including dynamical phase transitions [30, 31]. It consists of many hard particles diffusing in one dimension, and therefore has some similarities with the simple symmetric exclusion process (SSEP), generalised to continuous space. Dynamical phase transitions in the SSEP have been analysed in detail [5, 32, 29, 13, 33]. In particular, its behaviour on very large (hydrodynamic) length scales is described by macroscopic fluctuation theory (MFT) [21]. The general applicability of this theory means that its predictions apply also in the BHPM [30, 31]. A key prediction is that for large deviations with low values of the dynamical activity, the system becomes macroscopically inhomogeneous, which was identified in [30] as a form of phase separation.

This article extends previous work [30, 31] on the BHPM in two main directions. First, we analyse the macroscopically inhomogeneous regime, both numerically and analytically. We generalise existing MFT predictions [17, 25] to the BHPM and we show quantitative agreement with numerical results. This theory is accurate for fluctuations that involve smooth large-scale modulations of the density. We also analyse fluctuations that are too large to be captured by MFT, where one observes the formation of an extensive region without any particles at all: a macroscopic gap. Based on previous results for kinetically constrained spin models [34, 35, 36], we propose a simple interfacial model that captures qualitative features of this regime.

The second direction of this work is to show how the addition of control forces [37, 38] to the equations of motion of the system can be used to improve numerical convergence. It is known [38, 39] that for any given biased ensemble there is an optimal set of control forces for which numerical sampling of the rare events becomes trivial. While these optimal forces cannot usually be computed in complex physical systems, it is expected that adding non-optimal control forces can also improve the convergence of numerical calculations, via a form of importance sampling [40, 36, 41, 42]. We use theoretical arguments to derive approximations to the optimal control force, in two regimes that we have identified. We show that these control forces do indeed improve numerical performance, and this improvement is increasingly strong when we consider large systems. (This is because our approximations to the optimal control forces are increasingly accurate for large systems.)

The form of this paper is as follows. Sec. 2 introduces the models that we consider, and some of the quantities that we will measure. Sec. 3 collects properties of biased ensembles of trajectories. Sec. 4 gives an overview of the main theoretical results, before Sec. 5 and Sec. 6 describe the detailed calculations for the two macroscopically inhomogeneous regimes that we identify. We summarise our conclusions in Sec. 7.

2 Models

Consider $N$ hard particles moving in one dimension with periodic boundaries. Each particle has size $l_{0}$ and the position of particle $i$ at time $t$ is $x_{i}(t)$. We write $X=(x_{1},x_{2},\dots,x_{N})$ for a configuration of the system and x for a trajectory of the the system, over a time interval $[0,t_{\rm obs}]$. The particle motion is stochastic and obeys detailed balance with respect to an equilibrium distribution

[TABLE]

where $\beta$ is the inverse temperature, $Z$ is a normalisation constant, and $U$ is a pairwise additive potential energy $U=\frac{1}{2}\sum_{i\neq j}v(x_{i}-x_{j})$. As in [30, 31], we consider two variants of the system, which have very similar behaviour.

2.1 Monte Carlo dynamics

The MC variant of the model is a discrete-time Markov process. On each step, a particle (say $i$) is chosen at random and we propose to move it to a new position $x_{i}+\Delta x$ where $\Delta x$ is uniformly distributed in $[-M,M]$. Here, $M$ is a parameter of the model. We use Glauber dynamics so the move is accepted with probability

[TABLE]

where $\Delta U$ is the difference in energy between the current configuration and the proposed configuration. If the move is rejected then the configuration remains the same. After each attempted move, the time is incremented by $\delta t=M^{2}/(12ND_{0})$ where $D_{0}$ is the single-particle diffusion constant, which is also a parameter of the model. This ensures that in the dilute limit, particles diffuse independently between collisions, with diffusion constant $D_{0}$. We use Glauber dynamics because this facilitates later analysis, when we add control forces to the system, see Sec. 3.3.

This variant of the model considers hard particles, so the interaction potential is

[TABLE]

Suppose that the $j$th particle is selected to be moved in a given MC step, and suppose that the neighbouring particles have indices $p,q$. The probability that the proposed move does not result in two particles overlapping is

[TABLE]

The superscript $M$ indicates that this quantity depends on the MC step size, it is a label (and not any kind of exponent).

2.2 Langevin dynamics

All numerical results in this work use MC dynamics. However, it is convenient for theoretical analysis to consider a Langevin equation

[TABLE]

where $\eta_{i}$ is a standard Brownian noise. In this case the pair potential should be differentiable: we assume a regularised version of (3) such that $v(x)=\infty$ for $|x|<l_{0}$. Also there is some $l_{1}$ such that $v(x)=0$ for $|x|>l_{0}+l_{1}$, with $v(x)$ a continuous function for $l_{0}<x\leq l_{1}$, diverging as $x\to l_{0}$. This choice ensures that the separation between any pair of particles is always larger than $l_{0}$.

The similarity between the MC and Langevin models can be justified in the following way. In the Langevin model, take $l_{1}=l_{0}+M$ where $M$ is the step size in the MC model. The two models behave equivalently in the limit $M\to 0$: this can be verified by constructing the Fokker-Planck equation for the Langevin model and the corresponding master equation for the MC model, then taking the relevant limits. See [43].

2.3 Rescaled representation

Since we consider hard particles in $d=1$, the ordering of particles in the system is preserved.111 The particles exclude a volume $l_{0}$ which we assume throughout is larger than the MC step $M$. One may always map such a model to a system of point particles that move in a spatial domain of size

[TABLE]

We insist that the particle positions are ordered with $x_{1}<x_{2}<\dots$ (modulo periodic boundaries) in which case the position of the $j$th point particle is $\tilde{x}_{j}=x_{j}-jl_{0}.$ For the MC variant of the model, the equilibrium distribution (1) reduces to an ideal-gas distribution for the positions $\tilde{x}$. In some cases, this means that the rescaled system is simpler to analyse. However, we emphasise that the rescaled system and the original system contain exactly the same information.

In this rescaled representation, it is easy to see that varying $l_{0}$ in the original model simply shifts particles’ positions by constants that are independent of time. It follows that many properties of the system (including trajectories of individual particles) are independent of $l_{0}$.

2.4 Dynamical activity

In the following, we focus on ensembles of trajectories that are biased to low (or high) values of time-averaged measurements of dynamical activity. The definition of activity used in this work differs from [30, 31] – the choice used here does not change the qualitative behaviour but it makes it easier to analyse, both numerically and computationally. Large deviations for a different kind of dynamical activity have recently been analysed in a similar model [44].

The activity measures motion on a characteristic length scale $a$. We introduce a dimensionless parameter

[TABLE]

which is the ratio between $a$ and the average interparticle gap. For a trajectory x, we define

[TABLE]

where $r_{i}^{a}$ is the acceptance probability for an MC move of size $a$, as defined in (4). We allow the parameter $a$ that appears in the definition of $K$ to be different from the parameter $M$ that determines the size of MC moves, although our numerical results take $a=M$. Note also that while $K$ is defined in terms of the MC acceptance rate, it can be evaluated directly from particle trajectories, using (8). Thus, $K$ is a well-defined quantity in the Langevin variant of the model, as well as in the MC variant. Also, the value of $K$ only depends on gaps between adjacent particles and is therefore the same in the rescaled representation, or the original representation.

It is useful to define an intensive (and dimensionless) version of $K$ by dividing by the number of particles and by $t_{\rm obs}$:

[TABLE]

In large systems, the gaps between adjacent particles are exponentially distributed with mean $L_{\rm r}/N$. Hence the mean of $r_{i}^{a}$ is the probability that a randomly chosen gap ($y$) is larger than the proposed step ($z$):

[TABLE]

(Here and throughout, $\langle\cdot\rangle_{0}$ indicates an average in the equilibrium state of the system.) The integral gives $\Phi_{a}^{-1}(1-{\rm e}^{-\Phi_{a}})$, so one has (for large systems, $N\to\infty$)

[TABLE]

At low concentrations (small $\Phi_{a}$), particles diffuse almost independently and the activity $k$ is equal to unity. For high concentrations the mean activity (per particle) is reduced; it approaches zero as $\Phi_{a}\to\infty$ (in which case particles do not move at all).

3 Biased Ensembles of Trajectories

This work focusses on the distribution of the intensive activity $k[\mathbf{x}]$ as $t_{\rm obs}\to\infty$. In a system with $N$ particles, large deviation theory for this time-averaged quantity means that its probability density scales as

[TABLE]

where $I$ is the rate function. This is a large deviation principle, which holds for $t_{\rm obs}\to\infty$ at fixed $N$.

Evaluation of $I(k)$ gives the probability of rare events where the time-averaged activity takes a non-typical value. This section outlines several results from large deviation theory as it applies to ensembles of trajectories [5, 6, 39, 45], including connections to optimal-control theory [38, 46], and its application for numerical sampling [47, 40]. Readers familiar with this material may prefer to skip directly to the summary of main results in Sec. 4.

3.1 Biased ensembles

We define biased ensembles of trajectories according to standard methods [5, 6, 45], by modifying the probabilities of trajectories of the system. We use $\mathrm{d}P_{0}[\mathbf{x}]$ to indicate the (infinitesimal) probability that the system follows trajectory x. The meaning of this notation is that the expectation value of some observable $O$ can be expressed as

[TABLE]

where the integral runs over all possible trajectories, weighted by their probabilities. Now consider an ensemble in which the probability of trajectory x is biased according to its activity:

[TABLE]

with $Z_{s}=\langle{\rm e}^{-sK[\mathbf{x}]}\rangle_{0}$ for normalisation. By analogy with (13), averages in the biased ensemble are given by

[TABLE]

Since $K$ is extensive in time, it is useful to invoke an analogy between these biased ensembles and canonical ensembles in statistical mechanics, see [48, 15] for a discussion. This motivates us to define the dynamical free energy,

[TABLE]

The average of the intensive activity in the biased ensemble is denoted by

[TABLE]

(There should be no confusion between the mean activity $k(s)$ and the activity of an individual trajectory $k[\mathbf{x}]$.) The average is over trajectories of fixed length $t_{\rm obs}$ so $k(s)$ depends implicitly on $t_{\rm obs}$, as well as the parameters of the model. Note also that $\lim_{t_{\rm obs}\to\infty}k(s)=-\psi^{\prime}(s)/N$, where the prime indicates a derivative. One reason that these biased ensembles are useful is that for large $t_{\rm obs}$, typical trajectories taken from (14) are representative of the rare events associated with (12), evaluated at $k=k(s)$ [6, 23, 48, 45].

In analogy with thermodynamics, this derivative $\psi^{\prime}$ corresponds to the average value of an order parameter. The free energy is related to the rate function $I$ by Legendre transform $I(k)=\sup[-sk-\psi(s)]$. If results for ${k}(s)$ and $\psi(s)$ are available from numerical data then the rate function may be estimated parametrically as

[TABLE]

3.2 Dependence of averages on time

From the definition of the biased ensemble in (14), it follows that this ensemble has transient regimes when the time $t$ is close to $t=0$ or $t=t_{\rm obs}$. These transients can be characterised theoretically following [23, 39, 45]. To this end, consider a general observable quantity $\hat{z}$ that can be measured at some single time $t$ (for example, $\hat{z}$ might be the distance between two particles). The probability density for this quantity when evaluated at time $t=t_{\rm obs}$ is

[TABLE]

The probability density for $\hat{z}$ can also be averaged along the whole trajectory, which gives

[TABLE]

This is the probability that $\hat{z}$ has value $z$, if we measure at a time $t$ chosen uniformly at random from $[0,t_{\rm obs}]$. These distributions depend implicitly on $t_{\rm obs}$; their limits are well-defined as $t_{\rm obs}\to\infty$. The two distributions $P_{s,\mathrm{ave}}$, $P_{s,\mathrm{end}}$ are different in general; in particular, they have different limits as $t_{\rm obs}\to\infty$ because $P_{s,\mathrm{end}}(z)$ characterises the transient regime while $P_{s,\mathrm{ave}}(z)$ characterises typical times, away from the transient regimes.

3.3 Conditioning of Doob, guiding forces, and optimal control theory

It has been shown in recent years [39, 45] that properties of biased ensembles of the form (14) can be reproduced by considering the typical (unbiased) dynamics of an “auxiliary process” that has been modified to include additional “control forces”. For the Langevin process (5), the auxiliary process is

[TABLE]

where $V_{\rm opt}$ is an optimal control potential whose determination is discussed below.

In the following, we make extensive use of (non-optimal) control forces, to improve convergence of our numerical algorithms, following [47, 40, 36, 41]. We now present the associated theory. Some details of derivations are given in A.

The object of primary interest in this study is the biased probability distribution $P_{s}$ of (14). It is convenient to define a new biased distribution that is very close to $P_{s}$, but differs in the transient regimes close to $t=0$ and $t=t_{\rm obs}$. Let $V$ be a control potential, similar to $V_{\rm opt}$ in (21), but not necessarily optimal. Then define

[TABLE]

The constant of proportionality in this equation is fixed by normalisation, we do not write it explicitly in order to have a compact notation. Since they differ only in transient regime, the distributions $P_{s}$ and $\tilde{P}^{V}_{s}$ are equivalent for long trajectories, in the sense that they yield the same results for $k(s)$ and $\psi(s)$. Marginal distributions $P_{\rm ave}$ are also identical for $P_{s}$ and $\tilde{P}^{V}_{s}$ (as $t_{\rm obs}\to\infty$), but $P_{\rm end}$ is different in general.222One way to characterise the optimal controlled model is that it leads to $P_{\rm ave}=P_{\rm end}$ [40], so there is no transient regime in that case.

The equivalence of $P_{s}$ and $\tilde{P}^{V}_{s}$ means that we are free to choose $V$ in such a way that $\tilde{P}^{V}_{s}$ is easy to analyse, either numerically or theoretically. To this end, let $\tilde{P}^{V}[\mathbf{x}]$ be the probability of trajectory x under the Langevin dynamics (21), with $V_{\rm opt}$ replaced by $V$. Then, we show in A that $\tilde{P}^{V}_{s}$ can be interpreted as a biased ensemble for this controlled process, that is

[TABLE]

with

[TABLE]

Similar results have been derived in [38, 40, 36, 49, 41]. Note that (23) applies for the Langevin model (5), the analogous result for MC dynamics is given in A.

Comparing (23) with (14), one sees that ${\rm d}P_{0}$ has been replaced by ${\rm d}\tilde{P}^{V}$, which indicates that the model has been modified by including the control potential $V$ in the equation of motion. Also the exponential biasing factor in (14) has been modified to include the action ${\cal A}^{\rm sym}$. In numerical work, these modifications are simple to implement, so algorithms for analysing $P_{s}$ can also be used to analyse $\tilde{P}^{V}_{s}$. This holds for any $V$ which allows an enormous flexibility [40, 36, 41]. In particular, the optimal control ($V=V_{\rm opt}$) is the potential $V$ for which the factor ${\cal A}^{\rm{sym}}[\mathbf{x}]-sK[\mathbf{x}]$ in (23) evaluates to a constant value $\psi(s)t_{\rm obs}$, independent of x. This allows $V_{\rm opt}$ to be obtained by solving an eigenvalue problem, see A. If this optimal control potential is known then sampling from $\tilde{P}^{V}_{s}$ is trivial [47, 50]. More commonly, the optimal potential is not available, but even non-optimal controls can greatly improve the performance of numerical schemes [40, 36, 41]. The nature of optimal control forces in some physical model systems is discussed in [46].

3.4 Sampling of biased path ensembles

We use transition path sampling (TPS) [51] to generate representative trajectories from $P_{s}$ and ${\tilde{P}}^{V}_{s}$. Our TPS methodology is the same as [31, 9]. To summarise, TPS is an MC method for sampling trajectories of a fixed length $t_{\rm obs}$. In each step, one starts with a trajectory and proposes to change it in some way. This is called a TPS move. The proposal is generated by direct simulation of the model of interest (see below). The proposed trajectory is accepted or rejected according to a Metropolis criterion based on the relevant weighting factor [for example, ${\rm e}^{-sK}$ in the case of (14)]. This MC method is designed to obey detailed balance, which ensures that it samples (14) or (23), as required. This is a key strength of the method; another advantage is that standard MC tests for numerical convergence can be applied, see for example Sec. 5.3.

There is some flexibility as to the specific choice of TPS moves. In this work, we use shifting moves [51, 31, 9] and the size of each shift is chosen uniformly from the range $\tau_{\mathrm{B}}\pm 0.5\tau_{\mathrm{B}}$, except where stated otherwise. As usual in TPS, proposing larger shifts is desirable for rapid exploration of trajectory space, but tends to lead to more TPS moves being rejected. The best choice of shift size is a compromise between these two effects.

We note that population dynamics (cloning) methods [52, 53, 54] have also been widely used for numerical studies of large deviations, and guiding (control) forces have also been used in that case [40, 36, 41]. We comment on the strengths and weaknesses of the two approaches at the end of this work, in Sec. 7.

4 Overview of main results

We consider fluctuations of the dynamical activity $K$ in the BHPM, as defined in Sec. 2.4. All numerical results are obtained using the MC variant of the model. The behaviour of the mean activity is illustrated in Fig. 1. We represent the data in two different ways. Fig. 1(a) shows $\langle K/(Lt_{\rm obs})\rangle_{0}=\langle k\phi\rangle_{0}$, which is the average activity per unit length, as a function of the volume fraction $\phi=Nl_{0}/L$. Fig. 1(b) shows the activity per particle $\langle k\rangle$, and its dependence on $\phi$. Note that the activity of a typical particle $\langle k\rangle$ decreases with volume fraction, but the activity per unit length is non-monotonic. (At small volume fractions, the activity is proportional to the number of particles and hence to $\phi$; on the other hand, it decreases for large volume fractions, because particles start to obstruct each other.)

We now consider large deviations of $K$. Fig. 2 shows the behaviour of $k(s)$ and the corresponding estimate of the rate function [obtained by (18)], for a representative state point $\phi=0.7$ in systems of $N=28$ and $N=42$ particles. These results were obtained by TPS, we note that they depend on the trajectory length $t_{\rm obs}$ which is quoted in units of the Brownian time,

[TABLE]

This is a natural unit of time in the model, and is comparable to the time required for a particle to diffuse its own size, $l_{0}$.

As noted in Section 3, $k(s)$ is analogous to an order parameter in thermodynamics. This quantity decreases sharply for positive $s$. As explained in [30, 31], this is a signature of a dynamical phase transition in the BHPM, which occurs in the limit $N,t_{\rm obs}\to\infty$. Before embarking on a detailed analysis, we give a brief summary of the associated phenomena. The qualitative behaviour of $k(s)$ in a system with finite $N$ is shown in Fig. 3, which also shows typical trajectories of the system, as one passes through the phase transition. At the phase transition, the system becomes inhomogeneous [30, 31]. In this work, we emphasise that (for this model) there are two distinct classes of inhomogeneous state. There are states where the density is modulated in space, but particle spacings remain of order unity as $N\to\infty$. However, for larger $s$ (smaller $k[\textbf{x}]$), there are states where a significant fraction of the available space in system is taken up by a single interparticle gap. The two classes of inhomogeneous state are discussed in Secs. 5 and 6.

In Sec. 5 we review and extend some previous work [55, 56, 17, 32, 25, 30], which shows that inhomogeneous states with spacings of order unity appear on taking $N\to\infty$ with $s={O}(N^{-2})$. This is the regime described by macroscopic fluctuation theory (which can also describe the behaviour for small negative values of the bias). In this regime, the optimal control forces are long-ranged; they are attractive for $s>0$ and repulsive for $s<0$. It is the attractive forces that drive the phase separation transition. We show that using control forces in numerical sampling significantly improves their efficiency. For $s<0$ the system always remains homogeneous; as explained in [30] it is hyperuniform when $s$ is negative and of order unity, see also Sec. 5.2.3.

In Sec. 6, we discuss the behaviour on taking $N\to\infty$ with $s={O}(N^{-1}$). We explain that this is the regime in which we expect a macroscopic gap to take up a finite fraction of the system. By applying such control forces in numerics, we show that computational efficiency is significantly improved. In fact, this improvement is much larger than for the MFT regime. We discuss how parameters of the control force can be optimised for efficient sampling.

We note that all these results apply in limits where $s\to 0$ as $N\to\infty$. The tractability of these limits arises because the biases that are applied to these ensembles of trajectories are weak. For example, a central assumption of MFT [21] is that the system is in “local equilibrium”, which means that any finite region of the system can be characterised through its local density and current, and that the bias has a negligible effect on the short-ranged correlations between microscopic particles. In the regime with a macroscopic gap, there are deviations from local equilibrium but our central assumption is that these are localised near the interfaces, at the edges of the gap.

5 Diffusion governed (MFT) regime

This section discusses the regime where MFT applies [21]. We work in the rescaled representation of Sec. 2.3. MFT is valid on large (hydrodynamic) length and time scales, which are related by a diffusive scaling. That is, we introduce the (average) density $\overline{\rho}$ and define

[TABLE]

The hydrodynamic limit is $N\to\infty$ at fixed $\overline{\rho},\gamma_{\rm obs}$. One then takes a second limit, $\gamma_{\rm obs}\to\infty$, in order to access the relevant large deviations. To arrive at a consistent theory, we also rescale the biasing parameter [32] as

[TABLE]

which is held constant as $N\to\infty$.

5.1 Theoretical analysis of density fluctuations and optimal control potential, using MFT

We consider the statistics of the local density and current, which are the relevant hydrodynamic fields within MFT [21]. We analyse large deviations of the activity using a physical argument based on fluctuating hydrodynamics – the same conclusions can also be reached other methods [32], the specific case of the BHPM is discussed in [30]. The statistics of the density and current may be characterised by writing Langevin equations:

[TABLE]

where $D(\rho)$ and $\sigma(\rho)$ are a diffusivity and a mobility, and $\eta$ is a space-time white noise with mean zero and $\langle\eta(x,t)\eta(x^{\prime},t^{\prime})\rangle=\delta(x-x^{\prime})\delta(t-t^{\prime})$. The specific forms of $D$ and $\sigma$ for the BHPM are discussed below. We use Ito calculus.

The usual approach in MFT is to rescale the spatial domain $[0,L_{\rm r}]$ into the unit interval $[0,1]$ and also to rescale time. That is, define dimensionless coordinates on the hydrodynamic scale as $\hat{x}=x/L_{\rm r}$ and $\hat{t}=t/L_{\rm r}^{2}$, the corresponding current is $\hat{\jmath}=jL_{\rm r}$ (there is no rescaling of the density). Then (28) becomes

[TABLE]

where $\hat{\nabla}$ is a gradient with respect to $\hat{x}$, and $\hat{\eta}$ is a noise with zero mean and $\langle\hat{\eta}(\hat{x},\hat{t})\hat{\eta}(\hat{x}^{\prime},\hat{t}^{\prime})\rangle=\delta(\hat{x}-\hat{x}^{\prime})\delta(\hat{t}-\hat{t}^{\prime})$. (Hence one sees that $\hat{\eta}=L_{\rm r}^{3/2}\eta$.) It is apparent from (29) that MFT is a weak-noise theory that is valid on large length scales. For later comparison with numerics, it is convenient to quote all results in the original co-ordinates (without hats), but we emphasise that they are valid only on the hydrodynamic scale (which in this case will mean $s={O}(N^{-2})$, as usual in diffusive systems [32]). Within MFT it is consistent (by the local equilibrium assumption [21]) to approximate the activity $K$ from (8) as

[TABLE]

where $\kappa(\rho)$ is the average activity (per unit volume) of an equilibrium system with density $\rho$. That is, $\kappa(\rho)=\frac{1}{L_{\rm r}t_{\rm obs}}\langle K\rangle_{0,\rho}$; for the BHPM then (11) and $\rho=N/L_{\rm r}$ imply that

[TABLE]

For activity fluctuations in the SSEP, one takes instead $\kappa=2\rho(1-\rho)$ in which case the following theoretical analysis is similar to [32, 25]. For the purposes of this discussion, the important feature is that $\kappa^{\prime\prime}(\rho)<0$. For positive $s$, this means that the system undergoes a continuous phase transition accompanied by spontaneous symmetry breaking [17, 32, 30], see also [13, 57, 58]. For $s=0$ the system is homogeneous but for positive $s$, it becomes inhomogeneous, see Fig. 3. We now adapt previous MFT results to this setting before comparing with numerical results in Sec. 5.2.

5.1.1 Homogeneous phase

As in [17, 32, 59, 30], we analyse the homogeneous phase by writing $\rho(x,t)=\overline{\rho}+\delta\rho(x,t)$, and assuming that $\delta\rho$ is small. From (28) we have $(\partial/\partial t){\delta\rho}=-\mathrm{div}\,j$ with (at leading order in $\delta\rho$):

[TABLE]

From (30) then

[TABLE]

Fourier transforming as

[TABLE]

one has (for $q>0$)

[TABLE]

where $\tilde{\eta}_{q}$ is a complex-valued Brownian noise. [There is one noise for each positive wavevector, each noise is independent of all the others, and $\langle\tilde{\eta}_{q}(t)\tilde{\eta}_{q}^{*}(t)\rangle=\delta(t-t^{\prime})$.] Also,

[TABLE]

In (35,36), the different wavevectors are completely decoupled from each other. The result is that each Fourier component of the density evolves independently by a (complex-valued) Ornstein-Uhlenbeck (OU) process with a bias proportional to $|\rho_{q}|^{2}$. Biased ensembles for these OU processes can be analysed exactly by standard methods, see B for details. Using (97) with $\alpha=\kappa^{\prime\prime}(\overline{\rho})$ and $\omega=D(\overline{\rho})q^{2}$ and $\gamma=\sigma(\overline{\rho})q^{2}$, the result is that

[TABLE]

which is equivalent to the results obtained in [32, 30]. The corresponding result for the activity is

[TABLE]

which will be compared in Sec 5.2 with numerical data.

Note however that (37) is valid only if the argument of the square root is positive which requires $2s\kappa^{\prime\prime}(\overline{\rho})\sigma(\overline{\rho})q^{2}>-(Dq^{2})^{2}$; otherwise the OU process predicts a divergence in the density fluctuations which signals a breakdown of the quadratic expansion in $\delta\rho$. Recalling that $\kappa^{\prime\prime}<0$, this criterion is most stringent for the smallest wavevector $q=q_{1}=2\pi/L_{\rm r}$ so one sees that the range of validity is $s<s_{c}$ with

[TABLE]

as in [30, 32]. For $s<s_{c}$, the sum in (38) converges to a finite value so $k(s)\to\kappa(\overline{\rho})/\overline{\rho}$ as $N\to\infty$, while the sum gives a finite-size correction that has a universal form within MFT [32]. In this regime, the optimal-control potential required to generate typical trajectories of the biased ensemble is obtained from Equ. (99) as

[TABLE]

with

[TABLE]

This corresponds to a pairwise-additive interaction whose pair potential $v(x)$ is given by the inverse Fourier transform of $v_{q}$. If $s\neq 0$ then $v_{q}$ diverges as $q\to 0$ indicating that this interaction is long-ranged. As discussed in [60, 61], the pair potential decays as $v(x)\sim 1/(\log x)$ for separations $x$ that are large compared to the particle spacing (but small compared to $L_{\rm r}$). The potential is attractive if $s\kappa^{\prime\prime}<0$ and repulsive if $s\kappa^{\prime\prime}>0$. We again emphasise that this analysis requires a weak bias $s<s_{c}$.

5.1.2 Inhomogeneous phase

For $s>s_{c}$, a slightly different approach is required, which is related to the Landau-like theory of [13, 62] as well as earlier work [17, 25]. For simplicity, we assume in this calculation that $D=D_{0}$ is a constant (independent of $\rho$): this situation holds for both the BHPM and the SSEP. The generalisation to density-dependent $D$ is straightforward.

We write the probability distribution (14) for the biased ensemble as

[TABLE]

where the action $S$ can be obtained from (28) or directly from MFT [21] as

[TABLE]

where the integrand depends on $(x,t)$ through $\rho=\rho(x,t)$. We consider the system close to the transition and we derive a result analogous to Sec 3.2 of [25]. Our method is slightly different from that work; note also that [25] considers specifically the case of the SSEP, where $\kappa=\sigma$ is a quadratic function of $\rho$. The following calculation can be interpreted as the derivation of a Landau theory for a suitable order parameter: similar calculations for systems with open boundaries are considered in [13, 62, 63].

We find the path that minimises the action $S$. The minimum has $j=0$ and the associated density $\rho$ is independent of time, as one might expect since the biased ensemble is time-reversal symmetric. Hence

[TABLE]

where now $\rho=\rho(x)$ and we used (26). This action functional can be minimised numerically over the profile $\rho$, in this work we make an expansion [17, 25] that is valid close to the critical point and allows an analytic treatment. In this regime one can capture the behaviour of the inhomogeneous phase by considering a density profile

[TABLE]

where $A,B$ are variational parameters. As noted above, the instability of the system originates in the smallest wavevector $q_{1}$ but it is necessary [25] to consider also the second-smallest wavector $2q_{1}$ to obtain accurate results for the inhomogeneous phase, even at leading order. We take $A>0$ without loss of generality; the system has translational symmetry so one may equivalently replace $x\to x+\delta$ to obtain an equivalent profile, shifted in space.

Close to the transition, we anticipate that $A$ is small, and $B=O(A^{2})$ is even smaller [25]. Inserting (45) in (44) we expand up to terms of order $A^{4},A^{2}B,B^{2}$. We introduce the short-hand $\sigma_{0}=\sigma(\overline{\rho})$ and similarly for derivatives such as $\sigma^{\prime}_{0}=\sigma^{\prime}(\overline{\rho})$, and also for $\kappa$. The result is

[TABLE]

with $\lambda_{c}=\frac{D_{0}(q_{1}L)^{2}}{-2\kappa^{\prime\prime}_{0}\sigma_{0}}$, consistent with (39,27) [recall $\kappa^{\prime\prime}<0$ and $q_{1}L=2\pi$] and

[TABLE]

where $\kappa^{(n)}_{0}$ is the $n$th derivative of $\kappa(\rho)$, evaluated at $\rho=\overline{\rho}$. The action is straightforwardly minimised over $B$, followed by minimisation over $A$. We assume (as usual) that $\kappa_{0}^{\prime\prime}<0$ and that $(\lambda-\lambda_{c})\ll\lambda_{c}$ is small so that $\mu_{2}>0$; we also require that $\mu_{4}>\mu_{3}^{2}/\mu_{2}$, which is true for the BHPM (see below). Then for $\lambda<\lambda_{c}$ the action is minimised at $A=B=0$ and the system is homogeneous; while for $\lambda>\lambda_{c}$ the minimum occurs for

[TABLE]

This predicts the degree of inhomogeneity for $\lambda>\lambda_{c}$. One sees that $A=O(\lambda-\lambda_{c})^{1/2}$ as one should expect, since $A$ is the order parameter for a $\phi^{4}$-like theory and we are making a mean-field analysis of the critical point, similar to [13, 62, 63]. Recall, we assumed in this derivation that $D$ is a constant, independent of $\rho$.333To make contact with the analogous calculation for the SSEP in [25, Sec 3.2] we take $\kappa=2\rho(1-\rho)$, $\sigma=\rho(1-\rho)$, $D=1$, noting that the definition of $\sigma$ in that work differs from ours by a factor of $2$. Hence $\lambda_{c}=\pi^{2}/[2\rho(1-\rho)]$. After some algebra one finds $\mu_{4}-(\mu_{3}^{2}/\mu_{2})=\lambda_{c}/[8\rho^{2}(1-\rho)^{2}]$ and hence $A^{2}=8\rho^{2}(1-\rho)^{2}(\lambda-\lambda_{c})/\lambda_{c}+O(\lambda-\lambda_{c})^{2}$, consistent with that work

These results closely resemble those of [29, 13, 62, 63] which apply in systems with open boundaries (not periodic). Those works concentrated on transitions where a $Z_{2}$ particle-hole symmetry is spontaneously broken. Here we do not assume particle-hole symmetry, instead the (periodic) system spontaneously breaks translational symmetry, which corresponds to a $U(1)$ symmetry (see below).

The average activity may then be estimated from (30) by plugging in the most likely density profile, which yields (for $\lambda>\lambda_{c}$):

[TABLE]

The two corrections appear because the computation of $A,B$ is only valid when $\lambda-\lambda_{c}$ is small, and the restriction to the most likely profile neglects fluctuations which enter as corrections at $O(1/L)$.

5.2 Application to the BHPM

5.2.1 Connection to MFT:

The results of Section 5.1 are general within MFT, in the sense that we did not specify the functional dependence of $\sigma,\kappa$ on $\rho$ (we did assume $\kappa^{\prime\prime}<0$ but the case $\kappa^{\prime\prime}>0$ is a simple generalisation.) We now consider the specific case of the BHPM, in the rescaled representation of Sec 2.3. For the Langevin variant of the model in the limit of a hard-core potential, this means that the statistics of the density field are identical to an ideal gas. (Particles are indistinguishable so collisions between hard particles are equivalent to events where the particles pass through each other.) In this case

[TABLE]

As discussed in Sec. 2, the MC variant of the model is equivalent to the Langevin one in the limit of small steps $M\to 0$. Since $M$ is non-zero for our numerical work, we expect corrections to (50), but we neglect these in the following. (We expect them to affect the quantitative predictions of the theory, but they are unimportant at the level of accuracy that we consider.) The expression for $\kappa$ is given in (31) which is independent of $M$; this yields $\kappa^{\prime\prime}(\rho)=-a{\rm e}^{-\rho a}$.

An interesting feature of the BHPM is that while the functions $D$ and $\sigma$ have ideal-gas behaviour, the nonlinear behaviour of $\kappa$ is still sufficient to drive the transition to an inhomogeneous state. However, we emphasise that none of our theoretical analysis relies on the fact that $\sigma$ is linear. (Recall from above that the finite step size $M$ in our numerical work should result in corrections to $\sigma$, but this is not expected to affect the qualitative behaviour.)

The MFT analysis requires that $\lambda$ is held constant as $L_{\rm r}\to\infty$. Noting from (7) that $\rho a=\Phi_{a}$ and using (39) one sees that the homogeneous state is stable if

[TABLE]

Hence, the calculation of Sec. 5.1.1 is valid in the range $0<\lambda<\lambda_{\rm c}$, for which the control potential is attractive (so density fluctuations are enhanced). It is also valid for negative $\lambda$, where the control potential is repulsive and density fluctuations are suppressed.

For the inhomogeneous phase $0<(\lambda-\lambda_{c})\ll L_{\rm r}^{2}$ the results of Sec. 5.1.2 are relevant. A suitable (complex-valued) order parameter for the phase transition is

[TABLE]

Recalling (34,45), we identify $|{\cal M}|^{2}$ with $(A^{2}/4)$ in Sec. 5.1.2. The normalisation of (34) means that typical values of $\tilde{\rho}_{q}$ are $\mathcal{O}(1)$ in the homogeneous phase, so ${\cal M}\to 0$ as $L_{\rm r}\to\infty$. For the inhomogeneous phase then ${\cal M}$ is of order unity – it is a complex number and its phase indicates the location of high and low-density regions in the system. The system is invariant under translation so there is a $U(1)$ symmetry for the phase of ${\cal M}$, which is spontaneously broken when the system becomes inhomogeneous.

5.2.2 Numerical results

To characterise the dependence of the mean activity on the bias $\lambda$ we define

[TABLE]

where the function $k(s)$ was introduced in (17). We use TPS calculations to sample ensembles $P_{s}$ and $P_{s}^{V}$ as defined in Sec. 3.3. These computations are performed at fixed (finite) values of $N,t_{\rm obs}$. Numerical errors are discussed in Sec. 5.3, below. The methodology provides accurate results for the given values of $N,t_{\rm obs}$; we compare these with the predictions of MFT that are valid in the hydrodynamic limit $N,\gamma_{\rm obs}\to\infty$

Numerical results for the BHPM are shown in Fig. 4. In this case, results for $\lambda<250$ were obtained using TPS without any control forces, effects of control forces in this regime are discussed in Sec. 5.3. For $\lambda>250$ we used the control force defined in (80) below, see section 6 for a discussion.

The results of Fig. 4 are consistent with the asymptotic predictions of MFT. The state point is $\Phi_{a}=0.233$ so Equ. (51) predicts $\lambda_{\rm c}\approx 107$, consistent with the data. On general grounds one would expect $\lambda_{\rm c}$ to be of order unity; its large numerical value in this case arises partly from the factor of $2\pi^{2}$ in (51) and partly from the fact that $\kappa^{\prime\prime}(\rho)$ in (31) is numerically small, for these parameters.

Considering the results in more detail, Fig. 4(a) shows that for fixed $\gamma_{\rm obs}$, the function ${\cal K}(\lambda)$ shows a scaling collapse as $N$ is varied, consistent with the expected diffusive scaling. (This scaling was less clear in [31]. We suspect that this difference arises because the values of $t_{\rm obs}$ used in [31] were not scaled with system size.)

Fig. 4(b) shows data for a single system size, and increasing $t_{\rm obs}$. Taking $t_{\rm obs}\to\infty$ at fixed $N$, one expects convergence of ${\cal K}$ to a limiting function: that is ${\cal K}(\lambda)\to\tilde{\cal K}_{N,\infty}(\lambda)$, where the subscripts indicate that $N$ is finite but $t_{\rm obs}\to\infty$. The system is finite so $\tilde{\cal K}_{N,\infty}(\lambda)$ is smooth (analytic) [45]. The data in Fig. 4(b) are consistent with this theoretical prediction, as $t_{\rm obs}\to\infty$. In principle, convergence of this limit is expected for $\gamma_{\rm obs}\gg 1$; in practice, it is notable that convergence appears to be already achieved for $\gamma_{\rm obs}\approx 0.7$. We expect that this small numerical value occurs for similar reasons to the large numerical value of $\lambda_{c}$, particularly fact that the largest diffusional time scale is $\tau_{L}=1/(D_{0}q_{1}^{2})=L_{\rm r}^{2}/(4\pi^{2}D_{0})$ so that $D_{0}\tau_{L}/L_{\rm r}^{2}=1/(4\pi^{2})$ is numerically small.

MFT makes predictions about the behaviour of $\cal K$ at large $N$. In particular, taking $N\to\infty$ leads to $\tilde{\cal K}_{N,\infty}(\lambda)\to\tilde{\cal K}_{\infty,\infty}(\lambda)$, where the function $\tilde{\cal K}_{\infty,\infty}(\lambda)$ is predicted to be singular at $\lambda=\lambda_{c}$. For $\lambda<\lambda_{c}$, Equ. (38) predicts that $\tilde{\cal K}_{\infty,\infty}$ is independent of $\lambda$; the same equation also predicts the first correction to the large-$N$ limit, as in [32]. For $\lambda>\lambda_{c}$, a simple MFT prediction for $\tilde{\cal K}$ can be read from (49), see also [25]. The validity of this result is restricted to small $(\lambda-\lambda_{c})$, because of the simple ansatz (45). When comparing with numerics, we obtain a similar prediction by substituting (45) into (30) and using (48) to fix $A$. The integral in (30) is performed numerically and yields a prediction for $\tilde{\cal K}_{\infty,\infty}$. This is the prediction based on (48) that is shown in Fig. 4(c); it matches (49) when $(\lambda-\lambda_{c})$ is small.

Fig. 4(c) compares a numerical estimate of $\tilde{\cal K}_{N,\infty}(\lambda)$ with these MFT predictions. The finite-size correction term in (38) is negative and diverges at $\lambda_{c}$, where the homogeneous theory is breaking down. One sees that (38) gives the correct qualitative behaviour for small $\lambda$, but a quantitative agreement with numerical data would require consideration of higher-order corrections, see also [58]. The theory behind (48) is valid as $N\to\infty$ and does not include any finite-size corrections; it captures the steep decrease in ${\cal K}(\lambda)$ but is not quantitative. Following [13, 62, 63], one expects a critical region $(\lambda-\lambda_{c})=O(N^{-2/3})$ where neither of (38,48) is applicable; this is consistent with the data but a more detailed finite-size scaling analysis would be required to confirm it.

Fig. 5(a) shows the behaviour of the order parameter $\langle|{\cal M}|^{2}\rangle_{s}$, which increases sharply at the transition, and takes a value of order unity in the inhomogeneous phase, consistent with the theory. Fig. 5(b) shows a smoothed representation of the density for the trajectory in Fig. 3(c), defined as

[TABLE]

where the normalisation constant $z$ is chosen such that $\int\rho_{\rm smooth}(x)\mathrm{d}x=1$. This shows that the density is macroscopically inhomogeneous, but we emphasise that the density is positive everywhere, which means that there is no macroscopic gap (see Section 6). We note that this system is quite far from the critical point, so the analysis of Sec. 5.1.2 is not sufficient to predict the density profile, consistent with the fact that it does not show a sinusoidal dependence on $x$. Instead, the behaviour (in this rescaled representation) is that the density shows a rather narrow peak. This might be analysed by minimisation of the action in (44) but we postpone such a calculation to future work.

D shows similar results to those presented here, using the definition of the dynamical activity that was used in [30, 31]. The qualitative behaviour is the same.

5.2.3 Negative $s$ and hyperuniformity:

Within MFT, biasing this system towards higher activity leads to hyperuniformity, as discussed in [30]. This means that density fluctuations on large length scales are strongly suppressed [64]. To measure this, define the structure factor

[TABLE]

Consider the limit $L\to\infty$ so that $q$ can take arbitrarily small values, and write $S_{\infty}(q)=\lim_{L\to\infty}S(q)$. A hyperuniform state is one where $\lim_{q\to 0}S_{\infty}(q)=0$ [64]. Such states are not expected in finite-temperature equilibrium systems with short-ranged interactions (they require that the system should have a vanishing compressibility), but there are many interesting examples that occur in systems away from equilibrium [65, 66, 67, 68, 69]. Within the MFT analysis of Sec. 5.1.1, hyperuniformity arises because the optimal control potential for $\lambda<0$ includes long-ranged repulsive forces, as may be deduced from (41). Hyperuniformity is a well-known property of systems with such long-ranged forces [70, 71], where it is sometimes referred to as super-homogeneity [68]. Within the framework of Sec. 5.1.1 and using (100), one sees that

[TABLE]

see also [59, 30]. This indicates that the system is hyperuniform for $s<0$. Fig. 6(a) compares this prediction with the results from simulations, the suppression of $S(q)$ at small $q$ is clearly apparent. The agreement is good – we attribute the differences between theory and simulation to the fact that MFT requires $N,\gamma_{\rm obs}\to\infty$ but these quantities are both finite in the numerical results. Fig.6(b) shows the pair correlation function

[TABLE]

which is proportional to the probability that two particles have separation $x$ (in the rescaled representation of Sec. 2.3). For the unbiased case ($s=0$) then $g(x)=1$ for all $x$. On biasing to high activity $s>0$ one sees a reduction in $g(x)$ for small $x$, since particles feel an effective repulsion, which enhances the activity via (8). Similarly, for a bias to low activity then $g(x)$ is enhanced for small $x$, consistent with an effective attraction.

5.3 Improved TPS by adding control forces

5.3.1 Convergence of TPS and error analysis

As discussed in Sec. 3.3, we expect the addition of control forces to improve the efficiency of TPS sampling. Since TPS is an MC method (which in mathematics would be called a Markov chain Monte Carlo (MCMC) method), analysis of convergence and numerical errors is straightforward [72]. To characterise the efficiency of the method, it is useful to compute how many TPS moves are required for trajectories to decorrelate from each other. Let $\hat{K}_{n}$ be the value of the activity for the $n$th trajectory generated by TPS. We define a block-averaged activity

[TABLE]

As $m\to\infty$, this block average converges to $\langle K[\mathbf{x]}\rangle_{s}$. Its variance behaves as

[TABLE]

where $\sigma^{2}_{\rm TPS}$ is the asymptotic variance [72]. Smaller values of $\sigma_{\rm TPS}^{2}$ correspond to more efficient TPS sampling: in particular $\sigma_{\rm TPS}^{2}/\mathrm{Var}(\hat{K})$ can be used as a rough estimate of the number of TPS moves required to generate an independent sample.

For small $m$ then all trajectories in the block are similar and one expects $\mathrm{Var}\left(\overline{K}_{n,m}\right)$ to be close to $\mathrm{Var}(\hat{K})$, independent of $m$. In our numerical analysis, we often plot

[TABLE]

as a function of $m$, for which the expected behaviour is of the qualitative form

[TABLE]

This quantity approaches $\sigma_{\rm TPS}^{2}$ as $m\to\infty$, as it should.

A suitable error bar for a numerical estimate of $K$ is then $\Delta K=\sigma_{\rm TPS}/N_{\rm TPS}^{1/2}$ where $N_{\rm TPS}$ is the total number of TPS moves over which the data is averaged. This error estimate accounts for the correlations between TPS moves. All error bars for TPS measurements in this work are computed in this way, estimating $\sigma_{\rm TPS}$ by (59).

5.3.2 Numerical results for accelerated convergence

As explained above, the results for $\lambda<250$ in Fig. 4 were obtained using TPS without any control forces. We now show how the control forces that can be derived from MFT can lead to a more efficient estimate of the same result, in the homogeneous regime $\lambda<\lambda_{c}$. We have computed the asymptotic variance for the BHPM in the regime where the system is homogeneous. As a simple control potential we take the first term in (40), so the control potential only depends on the first Fourier component of the density:

[TABLE]

where $q_{1}=2\pi/L_{\rm r}$ is the smallest allowed wavevector. This choice for the control potential has no free parameters. For these homogeneous states, it successfully captures the essential physical effect of the long-ranged control potential.

We emphasise once again that the TPS method is valid as a method for sampling from $P_{s}$, independent of whether control forces are used [9, 24]. The question that we address here is the rate of convergence, as characterised by the asymptotic variance. Fig. 7 shows the improvement in TPS sampling obtained using the control potential (62), which is significant for positive $\lambda$. All these results are in the homogeneous regime $\lambda<\lambda_{c}$, where the simple control potential (62) is applicable. They provide clear evidence of a significant speedup, and demonstrate proof-of-principle for the method. However, the small values of $s$ in this regime means that this is a regime where numerical sampling is relatively easy, and the speedup by the control forces is relatively modest. Alternative methods for analysing convergence of TPS are discussed in C.

We note again that the results of Fig. 7 are restricted to $\lambda<\lambda_{c}$. For the inhomogeneous regime ($\lambda>\lambda_{c}$) the next section considers control forces that apply for $s=O(1/N)$, which corresponds to $\lambda=O(N)$. We have not attempted to derive control forces for $\lambda=O(1)>\lambda_{c}$, this would be an interesting question within MFT.

6 The regime with a single macroscopic gap

The results of Sec. 5 are based on MFT which is valid for $N\to\infty$ at fixed $\lambda$, as discussed above. We note that minimisation of the MFT action $S$ in (44) always predicts that the density $\rho$ is finite everywhere, which means in turn that the gaps between particles almost surely have sizes of order $\overline{\rho}^{-1}$, as $L_{\rm r}\to\infty$.

In contrast to this, Fig. 3(c) shows a trajectory in which a single gap takes up a finite fraction of the system. This section focusses on that regime. As before, we work in the rescaled representation of Sec. 2.3. At time $t$, suppose that the largest gap in the system has size ${\cal Y}(t)$. We define $Y(t)={\cal Y}(t)/L_{\rm r}$, which is the fraction of the system occupied by this gap. If $Y(t)$ is order unity then the gap is macroscopic, in the sense that ${\cal Y}(t)=O(L_{\rm r})$.

To investigate this regime, define a new rescaled biasing parameter

[TABLE]

This rescaled bias $h$ is analogous to $\lambda$ of Section 5. We consider the behaviour on taking $L_{\rm r}\to\infty$ at fixed $h=O(1)$.

6.1 Numerical results

Fig. 8(a) shows that for small $h$, the gap size $Y$ remains close to zero (in particular, for small fixed $h$, the average $\langle Y\rangle_{s}$ decreases with $L_{\rm r}$). However, for larger $h(\gtrsim 2)$ there is a sharp increase in $\langle Y\rangle_{s}$, which we interpret as opening of a single macroscopic gap. [Recall again Fig. 3(c).]

To understand the behaviour for small $h$, we use extreme value theory to estimate the expected size of the largest gap. The distribution of interparticle gaps is exponential with mean $L_{\rm r}/N=\overline{\rho}^{-1}$. Hence for large $N$ the largest gap $\cal Y$ has a Gumbel distribution with mean $(\log N+\gamma_{\rm E})/\overline{\rho}$ where $\gamma_{\rm E}\approx 0.577$ is the Euler-Mascheroni constant [73]. Hence $\langle Y\rangle_{0}=(\log N+\gamma_{\rm E})/N$ which for $N=90$ is $\approx 0.06$, consistent with Fig. 8.

Fig. 8(b) shows the behaviour of the activity. As $h$ increases from zero, there is an initial sharp decrease in activity which corresponds to the MFT transition to an inhomogeneous state. As $L_{\rm r}\to\infty$, this transition would move towards $h=0$, because the critical point $\lambda=\lambda_{\rm c}$ discussed in Sec. 5 corresponds to $h=O(1/L_{\rm r})$. However, the systems considered here are only moderately large, and the numerical value of $\lambda_{\rm c}$ is also quite large – the result is that the MFT transition happens at $h\approx 1$ for the system sizes considered here. In contrast, the largest gap opens at $h\approx 2.5$, where an additional feature in $k(s)$ is also observed (in the larger systems). We now present a theoretical analysis of this regime, and we compare the resulting theory with these numerical results.

6.2 Theory – interfacial model

We define a simple model that captures the qualitative behaviour of the system in the regime with a single large gap, building on recent work on kinetically constrained models [34, 35, 36]. We separate the system into a dense region and a large gap, and we focus on the behaviour at the edge of the gap, which is the interface between the two regions. Hence we refer to this as an interfacial model.

6.2.1 Derivation of interfacial model:

To motivate the model, assume that configurations containing a large gap have all the particles distributed in some (dense) region of size $L_{\rm r}[1-Y(t)]$, and that they are distributed at random throughout this region. The mean distance between particles within the dense region is

[TABLE]

with $\overline{\rho}=N/L_{\rm r}$ as above. We model the dynamics of $Y$ by a Langevin equation where both the bias and the diffusion constant depend on $Y$:

[TABLE]

Here $\eta$ is a standard Brownian noise. To fix the functions $b$ and $D_{y}$ we use the MC variant of the BHPM to estimate the first and second moments of the change in the gap size $Y$, in a single MC move.

The gap size changes only when one of the particles on the edge of the gap has an accepted move. Proposed MC moves that reduce $Y$ involve particles moving into the largest gap: these are accepted with probability $(1/2)$, by (2). Proposed MC moves that increase $Y$ involve particles moving towards the dense region of the system: some of these moves will be rejected due to collisions between particles. Since we assumed that particles are distributed at random in the dense region, the distance between neighbouring particles in this region is exponentially distributed with mean $\ell_{Y}$. Hence, for MC moves that act to increase $Y$, the fraction that is accepted is

[TABLE]

where the factor of $2$ again comes from (2). Hence, for MC moves in which the proposed particle is on the edge of the macroscopic gap, the mean change in the gap size is

[TABLE]

where we consider separately the situations where the gap size decreases (first term) or increases (second term). The integrals can be computed exactly but we focus on the limit where $M/\ell_{Y}$ is small (small MC moves). This limit is sufficient to explain the main features of the model. It yields

[TABLE]

Similarly the mean square displacement is

[TABLE]

The relevant MC moves happen with rate $w_{y}=2/\tau_{0}$ where the factor of $2$ arises because particles on either side of the macroscopic gap can both affect its size, and $\tau_{0}=M^{2}/(12D_{0})$ is the time increment associated with one attempted MC move per particle (see also Sec. 2.1).

Using that the macroscopic gap is of size ${\cal Y}=YL_{\rm r}$ and taking $L_{\rm r}\to\infty$ one arrives at the Langevin equation (65) with $b(Y)=\overline{\Delta x}w_{y}$ and $D(Y)=\overline{(\Delta x)^{2}}w_{y}/2$. Hence (assuming as above that $M/\ell_{Y}\ll 1$):

[TABLE]

Since we assume that the particles are distributed at random in the dense region, the activity of a trajectory is [by analogy with (30)]

[TABLE]

with $\kappa_{y}(Y)=(1-Y)\kappa(\overline{\rho}/(1-Y)).$ [To derive this, recall that $\kappa(\rho)$ is the activity per unit length for a system with density $\rho$. Here, the dense region of the system has size $L_{\rm r}(1-Y)$ and density $\overline{\rho}/(1-Y)$.] Hence from (31)

[TABLE]

The Fokker-Planck equation corresponding to (65) is

[TABLE]

where $\mathcal{P}=\mathcal{P}(Y)$ is the probability density for $Y$.

In kinetically constrained models [34, 35, 36], a similar interfacial model was derived, which gives semi-quantitative predictions for the system behaviour if $b$ and $D_{y}$ are taken as constants. The system considered here is different in that $b(Y)$ has a diverging negative value as $Y\to 1$ – this reflects the fact that as the largest gap approaches the size of the system, all the particles end up confined in a very small region.

6.2.2 Biased ensembles for the interfacial problem:

We now analyse the effects of biasing to low dynamical activity in the interfacial model. The dynamical free energy $\psi(s)$ of the interfacial model is obtained by finding the largest $\psi$ that solves the following eigenproblem

[TABLE]

The diffusive term is suppressed by a factor of $1/L_{\rm r}$ so we identify this as a small-noise problem that may be solved by saddle point methods. It is convenient to transform to a self-adjoint (Hermitian) form by ${\cal U}(Y)=\int_{Y_{0}}^{Y}b(Y^{\prime})L_{\rm r}/(2D_{0}){\rm d}Y^{\prime}$ and defining $\mathcal{Q}(Y)=\mathcal{P}(Y){\rm e}^{-{\cal U}(Y)}$. (The reference point $Y_{0}$ can be chosen arbitrarily, so $\cal U$ is fixed only up to an additive constant.) The eigenproblem (74) becomes

[TABLE]

with a dimensionless potential

[TABLE]

The final step of the derivation used (63).

For large $L_{\rm r}$, this eigenproblem can be solved by saddle-point methods. The last term in (75) is negligible when $L_{\rm r}$ is large. Also, the dominant eigenfunction $\mathcal{Q}$ is sharply-peaked at the minimum of ${\cal V}$, which we denote by

[TABLE]

Also $\psi=-{\cal V}(Y^{*})\overline{\rho}^{2}D_{0}$. Using (70,72) we obtain

[TABLE]

Minimising this potential we find that $Y^{*}=0$ for small $h$, but there is a threshold $h_{\rm c}$ above which $Y^{*}$ becomes non-zero. At the threshold, $Y^{*}$ increases continuously which means that ${\cal V}^{\prime}(0)=0$ for $h=h_{\rm c}$. The existence of a threshold is consistent with Fig. 8, the accuracy of the detailed predictions will be discussed below. Before that, we derive the effective potential that describes the state with $Y^{*}>0$.

6.2.3 Optimal control potential:

We present two possible methods for estimating the optimal control potential introduced in Sec. 3.3. The first is based on a physical argument: observe that the dense region of the system contains particles that are distributed as an ideal gas, so their pressure is

[TABLE]

Maintaining a gap of size $Y^{*}$ requires a control force that balances the pressure. Since the eigenvector $\mathcal{Q}$ is sharply-peaked at $Y^{*}$, the fluctuations of $Y$ are small in the biased ensemble, so the behaviour is relatively insensitive to the form of the control potential, as long as it produces the correct force in the typical states (which have $Y=Y^{*}$). Hence, a control potential that reproduces the correct statistics for $Y$ is

[TABLE]

where $c>0$ is a constant with units of inverse length – its interpretation is that there is a constant force $c/\beta$ that acts to increase the gap size. To determine $c$ we equate the force to the pressure required to stabilise a gap of size $Y^{*}$:

[TABLE]

In order to use (80) with (21), the potential $V$ must be expressed as a function of the particle positions: this is straightforward because $Y$ is the size of the largest interparticle gap, which is a simple function of the the particle positions.

The second method for deriving a suitable control potential is the standard mathematical approach: consider the adjoint (Hermitian conjugate) of the eigenproblem (74) which is

[TABLE]

Since the noise is weak, the expected solution is of the form ${\cal F}(Y)={\rm e}^{-L_{\rm r}g(Y)}$ and the optimal control potential may be identified from (88) as $V_{\rm opt}(Y)=2L_{\rm r}g(Y)$. Inserting the expected form for $\cal F$, retaining terms at leading order in $L_{\rm r}^{-1}$, and using $\psi=-{\cal V}(Y^{*})\overline{\rho}^{2}D_{0}$ from above, one recovers $V_{\rm opt}^{\prime}(Y^{*})=-\overline{\rho}L_{\rm r}/(1-Y^{*})$. This is consistent with (80,81) which together imply $V_{\rm opt}(Y)=-\overline{\rho}L_{\rm r}Y/(1-Y^{*})$ for $Y\approx Y^{*}$.

6.2.4 Comparison with numerical results

For the parameters shown in Fig. 8, Equ. (78) predicts $h_{c}=\Phi_{a}/(1-e^{-\Phi_{a}}-\Phi_{a}e^{-\Phi_{a}})\approx 4.9$. This overestimates the value of the bias at which a macroscopic gap appears. The reason is clear if one considers the behaviour close to the threshold. In the interfacial model, the state with $Y^{*}=0$ has the particles distributed homogeneously but the MFT analysis of Sec. 5 has already established that the system is not homogeneous for these values of the bias.

If the state with $Y^{*}=0$ is already inhomogeneous, one sees that the probability of opening up a macroscopic gap will be enhanced, because the gap will likely appear at a location where the density is already low. Our conclusion is that the interfacial model predicts the existence of a threshold $h_{\rm c}$ at which a macroscopic gap appears, which is consistent with the numerical data. However, the assumption within the model that the dense region of the system is homogeneous is not accurate enough for the model to deliver quantitative predictions. In the following subsection, we show that despite these shortcomings, the optimal control potential predicted by the interfacial model is sufficiently accurate to significantly improve numerical sampling. In this sense, the interfacial model does capture the essential physical features of the regime with a macroscopic gap.

6.3 Improvement in sampling by control forces

We have performed TPS sampling using the control potential (80). The relation (81) is confirmed by our numerical results, in that a control potential with this value of $c$ leads to a typical largest gap of size $Y^{*}$. Fig. 9 shows the improvement in TPS sampling that is obtained with this control potential, which is more than an order of magnitude, even for small systems. The parameter $c$ in (80) is varied, in order to obtain the maximal speedup. For larger systems, the improvement increases rapidly – we are not able to quantify the speedup because (for example) the results shown in Fig. 8 would require a prohibitively large computational effort, if control forces were not used. The reason is that the macroscopic gaps that appear in those systems are extremely rare under the natural dynamics, so that TPS moves tend to be rejected if one uses a system without a control potential. We also note from Fig. 9 that significant speedup is possible for control forces that are not optimal, as emphasised in [40].

This improvement that is available from control forces also enables us to investigate what value of $c$ is most effective for improved sampling. As noted in Sec. 3.3, if one uses the optimal guiding force, the distributions $P_{\rm ave}$ and ${\tilde{P}}_{\rm end}$ of (19,20) coincide with each other, for all observable quantities. Recall that $P_{\rm ave}$ is independent of the guiding force but ${\tilde{P}}_{\rm end}$ is evaluated in a system with control forces, which does depend on the choice of these forces. It was suggested in [40, 36] that a suitable method for choosing approximate (non-optimal) control forces is to adjust their parameters to make the distributions $P_{\rm ave}$ and ${\tilde{P}}_{\rm end}$ as similar as possible.

This hypothesis is tested in Fig. 10. We first consider the distribution of $Y$, the largest interparticle gap. In this case one sees that the control force that gives the best overlap of $P_{\rm ave}$ and ${\tilde{P}}_{\rm end}$ is $c=5$, which is larger than the force which gives the most efficient sampling (this is $c=3$, from Fig. 9). We also consider the distribution of $Y_{2}$, which is the second largest interparticle gap, measured relative to the system size $L_{\rm r}$. For this quantity, the distributions overlap best at $c=3$, where the sampling is most efficient.

The conclusion of this analysis is that maximising the overlap of $P_{\rm ave}$ and ${\tilde{P}}_{\rm end}$ for any single observable does not guarantee that the distributions for other observables should overlap. This cautions against placing too much faith in the overlap of any single distribution, as an indicator of where sampling is most efficient.

As an alternative method for estimating which control force is optimal, one may consider the statistics of the action, as suggested in [10]. Let $\langle B\rangle_{s,V}=\int B[\textbf{x}]d\tilde{P}^{V}_{s}[\textbf{x}]$ be the average of the observable $B$ with respect to the distribution $\tilde{P}^{V}_{s}$ of (22). If $V$ is the optimal control then

[TABLE]

The suggestion of [10] is that optimising $V$ to achieve equality in (83) can be used to obtain good sampling. Note that there are many control forces that can achieve equality in (83). This situation is to be contrasted with the general inequality [38, 42]

[TABLE]

where the average is with respect to the controlled dynamics, without any biasing: $\langle B\rangle_{V}=\int B[\textbf{x}]d\tilde{P}^{V}[\textbf{x}]$. In (84), equality can only be achieved if $V$ is the optimal control potential, this can be checked by noting that $\psi(s)=\lim_{t_{\rm obs}\rightarrow\infty}t_{\rm obs}^{-1}\log\langle\mathrm{e}^{-sK+\mathcal{A}^{\mathrm{sym}}}\rangle_{V}$ and using Jensen’s inequality. On this basis one might expect that maximising the right hand side of (84) would give the best sampling.

Results for the averages in (83,84) are shown in Fig. 11. Contrary to the situation in [10], there is no value of $c$ for which equality is achieved in (83). However, we note that the most efficient sampling takes place for $c=3$, which is the value where the average on the right hand side of (83) is closest to $\psi(s)$, consistent with the proposal [10] that equality in (83) is a desirable feature. One also sees that the right hand side of (84) is decreasing in $\psi$ for all $c>1$. Thus, $c=1$ gives the best bound on $\psi$ but it does not achieve the best sampling, contrary to the intuitive expectation stated above.

Based on Figs. 9,10,11, our conclusion in this Section is that no single prescription seems satisfactory for determining the best choice of control force $V$ in practical situations such as this one, and some trial-and-error is still necessary in this process.

7 Conclusions

We have given a detailed analysis of the behaviour that was summarised in Fig. 3, including discussion of hyperuniform states that appear when states are biased to high activity, and inhomogeneous states with low activity. We have discussed the existence of two inhomogeneous regimes, with $s=O(N^{-2})$ (MFT regime) and $s=O(N^{-1})$ (macroscopic interparticle gap).

We have shown that control forces can be used to improve numerical sampling of these ensembles [40, 36, 41, 44]. In the MFT regime where the system is homogeneous, these optimal control forces are long-ranged. On biasing to low activity, these forces are attractive and drive the formation of inhomogeneous states. Using these guiding forces (in the homogeneous state) leads to an improvement in sampling efficiency. However, for these small values of $s$, the effect of the bias is weak, so sampling is already possible without these forces.

In the regime with a macroscopic interparticle gap, we have argued that a form of interfacial model can capture some features of the system, similar to [36]. Using this model to infer a suitable control force leads to an improvement in sampling efficiency that is more than a factor of 10 in small systems. For large systems, the computations that we present would be prohibitively expensive without these control forces. We have discussed how the parameters of the control force might be optimised. In particular, we find that the simple criterion of [36, 40], to match the distributions $P_{\rm ave}$ and $P_{\rm end}$ is not optimal for the cases considered here: since the control force is very simple we have instead optimised its free parameter by hand. Further work would be valuable, to understand how to infer control forces that improve sampling efficiency. Some generic aspects of optimally-controlled models are discussed in [46], where they are referred to as auxiliary models.

Comparing the TPS method with cloning methods [52, 53, 54, 40, 36, 41], we note that while TPS gives results for finite $t_{\rm obs}$, cloning provides direct access to the limit $t_{\rm obs}\to\infty$, which is the limit where large-deviation theory applies. On the other hand, the detailed-balance property of TPS [51] means that it samples directly from $P_{s}$ or $\tilde{P}_{s}^{V}$ as defined in (14,22). By contrast, cloning methods do not sample directly from a target distribution; they do allow estimation of averages with respect to these distributions but the associated statistical estimators have systematic errors (bias) which only disappear as the population size tends to infinity [74, 33]. Estimation of statistical uncertainties is also simpler for TPS, see Sec. 5.3. There are also other methods for sampling large deviations, some of which require accurate representations of an optimally-controlled dynamics in order to achieve accurate results [47, 75, 42], see also [44]. Such methods are attractively simple, but accurate representations of optimally-controlled dynamics may be challenging in complex systems. We emphasise once more that the role of control forces in TPS and cloning is to improve convergence, but accurate results are still available without obtaining the optimal control force.

On physical grounds, it is notable that the rich physics of inhomogeneous and hyperuniform states in the BHPM occurs for very small values of the bias parameter $s$, which are at either $O(N^{-2})$ or $O(N^{-1})$. The strong response of the system to these biasing fields has its origin in hydrodynamic modes. Many theories of biased ensembles assume the existence of a gap in the spectrum of the generator of the relevant stochastic process. Here the gap size is vanishing as $N\to\infty$, because of slow (diffusive) hydrodynamic modes. The MFT approach is to rescale (speed up) time so that one is restricted to hydrodynamic time scales, but the gap for the generator is restored.

The fact that the behaviour originates on the hydrodynamic scale also explains why MFT predictions are universal, in that they depend on diffusive scalings but not on microscopic details of the model. The predictions for the behaviour for $s=\mathcal{O}(N^{-1})$ are not universal in the same sense, but the simplicity of the interfacial model indicates that they may arise generically in systems with sharp interfaces between coexisting phases (see also [36]).

We thank Takahiro Nemoto and Vivien Lecomte for discussions about the use of control forces to aid numerical sampling. Jakub Dolezal was supported by a studentship from the EPRSC, reference EP/N509620/1. The figures in this article were made using Matplotlib [76].

Appendix A Biased path ensembles

This appendix includes a derivation of (23) and its analogue for MC dynamics.

A.1 Langevin dynamics

As in the main text, let $\mathrm{d}\tilde{P}^{V}[\mathbf{x}]$ be the probability for trajectory $\mathbf{x}$ under the controlled process (21), but with $V_{\rm opt}$ replaced by some general (possibly non-optimal) potential $V$. Then standard path-integral arguments (see eg [40]) show that

[TABLE]

where the normalisation constant has been omitted for ease of writing and

[TABLE]

Combining (14) with (85) yields

[TABLE]

which means that the ensemble (14) which was obtained from $P_{0}$ by biasing $K$ can also be obtained (exactly) by a suitable biasing of $\tilde{P}^{V}$. Recalling that the product of $\dot{x}$ and $\nabla V$ in (86) is to be interpreted in the Ito sense and using Ito’s formula for $\mathrm{d}V/\mathrm{d}t$, we obtain ${\cal A}[\mathbf{x}]=\frac{1}{2}[V(t_{\rm obs})-V(0)]+{\cal A}^{\rm{sym}}[\textbf{x}]$ where ${\cal A}^{\rm{sym}}$ is given by (24). Using this with (87) and (22) yields (23).

Finally, using (24) with the observation that ${\cal A}^{\rm{sym}}[\mathbf{x}]-sK[\mathbf{x}]$ is constant and equal to $\psi(s)t_{\rm obs}$ for $V=V_{\rm opt}$, one sees that

[TABLE]

where $u$ is the solution with largest eigenvalue $\psi$ of the eigenproblem

[TABLE]

This is a tilted Fokker-Planck equation in its adjoint form, see for example [45].

A.2 MC dynamics

Analogous formulae hold for the (discrete-time) MC variant of the model. Let $p(X_{k}|X_{k-1},i_{k})$ be the probability that the system is in state $X_{k}$ at step $k$, given that it was in state $X_{k-1}$ at step $k-1$, and that the particle proposed to be moved on step $k$ was $i_{k}$. Note that this $p$ is a normalised probability for $X_{k}$ and $p(X_{k}|X_{k},i_{k})$ is generically finite. Also let $p^{V}(X_{k}|X_{k-1},i_{k})$ be the analogous quantity for the controlled model. Then the analogue of ${\cal A}^{\rm{sym}}$ is

[TABLE]

For the logarithm to be finite, it is important that $p(X_{k+1}|X_{k},i_{k})$ should not be zero (except if $p^{V}(X_{k+1}|X_{k},i_{k})=0$ also). This is the reason to use the Glauber criterion in (2) instead of the Metropolis condition (because using Metropolis may result in $p(X|X,i)=0$ for some choices of $X,i$ but $p^{V}(X|X,i)\neq 0$).

Note also that ${\cal A}^{\rm{sym}}_{\rm MC}$ depends on which moves were proposed, as well as the actual sequence of states in the trajectory. If the control force $V$ has a very simple form then it is possible to write an equation similar to (90), in which the action depends only on the actual sequence of states. This gives some improvement in numerical sampling but is restricted to simple cases, for example where $V$ is a linear potential.

For the numerical results for the MC variant of the BHPM, we use control forces that are introduced by replacing $\beta\Delta U$ in (2) by $\beta\Delta U+\Delta V$ where $\Delta V$ is the change in the control potential, for the proposed move. We note that the optimal auxiliary model for such a system would require that we take instead

[TABLE]

where $r^{a}(t_{\pm})$ are the values of $r^{a}$ just before and after the proposed move. For small $\Delta V$ and small $\delta t$, this is equivalent to replacing $\beta U\to(\beta U+V)$ in (2) and it is also equivalent to the Langevin case.[43].

We emphasise that for any $V$, the TPS method targets the distribution $\tilde{P}^{V}_{s}$ and provides accurate results (as long as sufficiently many TPS moves are performed). However, it is possible that we might have observed faster convergence if we had used (91) instead of simply including the control potential in (2). To check this, we tested an algorithm based on (91) for several representative cases; the differences in performance were within the statistical uncertainty of our estimates of the asymptotic variance.

Appendix B Complex Ornstein-Uhlenbeck processes

We collect some results for biased ensembles constructed from Ornstein-Uhlenbeck (OU) processes, see for example [45, Sec. 6.2]. Suppose that $z$ is a complex number which evolves by the complex OU process

[TABLE]

where $\omega,\gamma$ are real positive constants and $\eta$ is a complex-valued white noise. That is, $\eta=\eta_{r}+{\rm i}\eta_{i}$ with real-valued noises $\eta_{r},\eta_{i}$ that satisfy $\langle\eta_{r}(t)\eta_{r}(t^{\prime})\rangle=\frac{1}{2}\delta(t-t^{\prime})=\langle\eta_{i}(t)\eta_{i}(t^{\prime})\rangle$ and $\langle\eta_{i}(t^{\prime})\eta_{r}(t)\rangle=0$. Then also $\langle\eta(t)\eta^{*}(t^{\prime})\rangle=\delta(t-t^{\prime})$. Writing $z=x+{\rm i}y$ one has independent equations of motion for $x$ and $y$. The corresponding Fokker-Planck equation for the probability density $P=P(x,y)$ is

[TABLE]

whose stationary distribution is $P_{0}\propto\exp(-\omega(x^{2}+y^{2})/\gamma)$. Alternatively one may use the calculus of complex variables and consider a probability density defined as $Q=Q(z,z^{*})$ which obeys

[TABLE]

The stationary solution is $Q_{0}\propto\exp(-\omega z^{*}z/\gamma)$ which is (of course) equivalent to $P_{0}$ as given above. The following results can be derived by considering separately the real and imaginary parts of $z$ but we use the complex variable representation, which simplifies the analysis.

For biased ensembles of the form (14) with $K=\alpha\int_{0}^{t}{}_{\rm obs}z^{*}(t)z(t)\mathrm{d}t$, the dynamical free energy can be obtained by solving the eigenproblem

[TABLE]

It is easily verified that the eigenfunction with maximal eigenvalue is

[TABLE]

which is valid for $2s\alpha\gamma>-\omega^{2}$ (otherwise the eigenvalues are not bounded above and the dynamical free energy does not exist). The corresponding eigenvalue is

[TABLE]

To obtain the optimal control force one should solve the adjoint eigenproblem

[TABLE]

whose solution is ${\cal F}\propto\exp\left[-\frac{z^{*}z}{2\gamma}\left(\sqrt{\omega^{2}+2s\alpha\gamma}-\omega\right)\right]$. Note that ${\cal F}\propto Q/Q_{0}$ which follows because the underlying equation is reversible (obeys detailed balance). The optimal control potential is $V_{\rm opt}=-2\log{\cal F}$ (up to an arbitrary additive constant) which yields

[TABLE]

Away from transient regions, the distribution of $z$ in the biased ensemble is $P_{\rm ave}(z,z^{*})\propto{\cal F}Q$ so

[TABLE]

For the discussion here, the case of primary interest is when $s\alpha<0$, in which case the control potential $V_{\rm opt}$ has negative curvature and guides the system towards increasingly large values of $z$. As $s\alpha$ tends to $-\omega^{2}/(2\gamma)$ one sees that the variance of $P_{\rm ave}$ diverges. If the original equation (92) was derived by linearisation at small $z$, then this divergence indicates the breakdown of the linear approximation, within the biased ensemble. This is the situation discussed in Sec. 5.1.

Appendix C Convergence of TPS

In order to measure the improvement in sampling that is achieved by guiding (control) forces, we discuss in the main text the asymptotic variance $\sigma^{2}_{\rm TPS}$, see (59). This quantity requires a large amount of TPS data to evaluate it, but does give a reliable estimate of the effort required to obtain an independent sampling from a biased trajectory ensemble. As an alternative, we also consider the autocorrelation function. In the notation of (58) let

[TABLE]

where the average is over many realisations of the TPS algorithm. One sees that $\chi^{\rm TPS}_{m}=\sum_{i,j=1}^{m}\langle\hat{K}_{i}\hat{K}_{j}\rangle-\langle\hat{K}_{i}\rangle\langle\hat{K}_{j}\rangle$ is related to a sum of $C(n)$ over the lag time $n$. Fig. 12 shows results for this correlation function. As in figure 9, one concludes that the sampling is most effective for $c=3$, since the correlations decay most quickly when the control force has this strength. Compared with the asymptotic variance $\sigma^{2}_{\rm TPS}$, results for the autocorrelation function are somewhat easier to obtain in practice. The difficulty is that $\sigma^{2}_{\rm TPS}=\sum_{m=-\infty}^{\infty}C(m)$ has contributions from weak correlations at large $m$: accurate estimation of these (weak) correlations requires very long TPS runs.

Another approach is to consider what fraction of TPS moves are accepted, and how this is affected by the guiding forces. In general, TPS acceptance rates are not reliable as indicators of convergence. For example short shifting moves lead to slow decorrelation of the trajectories, while longer trajectories may decorrelate the trajectory more quickly, even if the acceptance probability is somewhat lower. Hence, if a control force leads to acceptance of longer shifting moves then this can still improve sampling, even at the cost of a lower overall acceptance rate.

Despite these limitations, there is useful information available by monitoring TPS acceptance rates. For TPS with control forces in place, it follows from (23) that a proposed trajectory is accepted with probability

[TABLE]

where $\Delta K$ is the change in activity between the original and proposed trajectory, and similarly $\Delta{\cal A}^{\rm sym}$ is the change in ${\cal A}^{\rm sym}$. For TPS with the optimal control potential then ${\cal A}^{\rm sym}_{\rm MC}-s\Delta K=\psi(s)$ for every trajectory so the acceptance probability is unity. That is

[TABLE]

Joint probability density functions for accepted values of $\Delta{\cal A}^{\rm sym}_{\rm MC}$ and $s\Delta K$ are shown in Fig. 13. The relationship (103) is indicated. There are two effects at play here. For control forces that are close to optimal, the distribution concentrates close to (103). On the other hand, larger control forces tend to suppress the total acceptance, because the forces are not optimal. The most efficient sampling occurs in an intermediate regime. In this case, we find that the the optimal regime is when the typical values of $s\Delta K$ and $\Delta{\cal A}^{\rm sym}$ are of similar sizes.

Appendix D An alternative measure of dynamical activity

Previous work has considered large deviations of the dynamical activity in this system [30, 31], but using a different measure of activity, which is defined in terms of squared particle displacements. One separates the time interval $[0,t_{\rm obs}]$ into $S$ segments, each of length $\Delta t=t_{\rm obs}/S$. Then define

[TABLE]

where $t_{j}=j\Delta t$ and $\Delta\bar{x}_{j}$ is the displacement of the centre of mass of all particles, between times $t_{j-1}$ and $t_{j}$. The activity $K$ of (11) depends on the characteristic length $a$, while $K_{\mathrm{msd}}$ depends on the parameter $\Delta t$. To obtain a corresponding length one may define $a_{\rm msd}=\sqrt{2D_{0}\Delta t}$ where $D_{0}$ is the single-particle diffusion constant. For a direct comparison between $K$ and $K_{\rm msd}$ it is natural to take $a_{\rm msd}\approx a$ since this means that both activity measures are sensitive to motion on the same length scales.

Analogous to (14) we define a biased ensemble with a bias parameter $s_{2}$, as

[TABLE]

Also define $\lambda_{2}=s_{2}L_{\rm r}^{2}/D_{0}$, analogous to (27), and $k_{\rm msd}=K_{\rm msd}/(Nt_{\rm obs})$. With these definitions, Fig. 14 shows that the ensemble of (105) has the same qualitative features as biasing by $K$. Specifically, Fig. 14(a) is analogous to Fig. 4(a) and Fig. 14(b) is analogous to Fig. 5(a). The data collapses when plotting these results as a function of $\lambda_{2}$, consistent with the MFT predictions.