Mechanical analysis of a dynamical phase transition for particles in a channel

We analyse biased ensembles of trajectories for a two-dimensional system of particles, evolving by Langevin dynamics in a channel geometry. This bias controls the degree of particle clustering. On biasing to large clustering, we observe a dynamical phase transition where the particles break symmetry and accumulate at one of the walls. We analyse the mechanical properties of this symmetry-broken state using the Irving-Kirkwood stress tensor. The biased ensemble is characterised by body forces which originate in random thermal noises, but have finite averages in the presence of the bias. We discuss the connection of these forces to Doob's transform and optimal control theory.

We analyse biased ensembles of trajectories for a two-dimensional system of particles, evolving by Langevin dynamics in a channel geometry. This bias controls the degree of particle clustering. On biasing to large clustering, we observe a dynamical phase transition where the particles break symmetry and accumulate at one of the walls. We analyse the mechanical properties of this symmetry-broken state using the Irving-Kirkwood stress tensor. The biased ensemble is characterised by body forces which originate in random thermal noises, but have finite averages in the presence of the bias. We discuss the connection of these forces to Doob's transform and optimal control theory.
In this work, we use numerical simulations to analyse a dynamical phase transition in a system of interacting particles, moving in a two-dimensional channel. The particles follow Langevin dynamics (with inertia); the ensembles are biased by a measurement of particle clustering, which plays a similar role to the dynamical activity considered in previous work [11,28,40]. On biasing to large clustering, the system undergoes a dynamical phase transition and spontaneously breaks symmetry, leading to an increased density near one of the walls of the channel.
In contrast to exclusion processes and other lattice models, this system can be analysed at a mechanical level, including the interparticle forces and the local stress tensor [41]. We exploit this approach to gain phys-ical insight into the phase transition. In the symmetrybroken state, we find a stress gradient that extends into the bulk of the system; this must be balanced by a suitable body force, to sustain the asymmetric density. We show that this force comes from thermal noise forces, which develop non-zero average values in response to the bias. Physically, this shows that large deviation events with increased clustering involve non-typical instances of the thermal noise forces, which push particles towards the walls of the channel, where they become localized. This result reveals the physical origin of the Doob forces, and it also provides a way to measure them directly, which is not possible with standard methods. We emphasize that our analysis follows directly from the (Newtonian) equations of motion: this mechanical perspective should therefore be applicable to a broad array of rare events and large deviations.
In the following, we define the model and relevant observables in Sec. II. Results are presented in Sec. III, including the existence of the dynamical phase transition, the behaviour of the stress, and the connection to fluctuating hydrodynamics. We conclude in Sec. IV by summarising the consequences of these results, and their relevance for future work.

A. System
We consider N particles with mass m and diameter l 0 , interacting by a Weeks-Chandler-Andersen (WCA) potential [42] of strength , in a box of size L x × L. It has periodic boundaries in the y-direction, and walls at Particle i has position r i = (x i , y i ) and momentum p i = mṙ i . Its equation of motion iṡ where U int and V w are the potential energies for particleparticle and particle-wall interactions respectively (see below); alsox is a unit vector in the x-direction, γ is the frictional damping rate, T is the temperature, and the η i are independent unit white noises with η µ i η ν j = δ ij δ µν δ(t − t ), where Greek indices indicate Cartesian components. To maintain a compact notation, T denotes the thermal energy, that is T = k B T ph where T ph is the physical temperature and k B is Boltzmann's constant. It is also common to use γ to denote a friction constant: our frictional damping force is mγṙ i so the friction constant is mγ in our notation. The interaction potential between particle i and the walls is of truncated Lennard-Jones type: where Θ is the Heaviside function, |x i − x w | is the distance from the particle to the nearest wall, the cutoff is l cut = 3l 0 /2, and α is chosen to ensure continuity of the potential at the cutoff (α = 0.080). We take L x = L + 2∆L, where ∆L = l 0 (2 1/6 − (1/4)) accounts for the excluded volume of the walls. Specifically, we vary N and L together, while fixing the (dimensionless) particle concentrationρ = N l 2 0 /L 2 and also ∆L. This choice ensures that the local particle concentration remains constant in the bulk of the system (away from the walls).
The interaction potential between particles i and j is a Weeks-Chandler-Andersen (WCA) potential [42]: where r ij is the distance between the particles and l WCA = 2 1/6 l 0 is the usual WCA cutoff. Hence the potential energy associated with these interactions is In addition to the dimensionless concentrationρ defined above and the number of particles N , a natural set of dimensionless control parameters of the system is obtained by identifying v 0 = T /m as the thermal velocity, and rescaling all variables in terms of l 0 , v 0 , T . One obtains two dimensionless parameters: The parameter˜ determines the strength of interactions, whileγ determines how strongly damped is the particle motion. Since D 0 = T /(mγ), we can also writẽ γ = (l 0 v 0 /D 0 ), so the damping determines the ratio of the thermal velocity to the diffusion constant (as usual). The simulation time step is δt, it is fixed by taking the dimensionless parameter (v 0 δt/l 0 ) = 0.002.
Throughout this work we takeρ = 0.48, a moderate density where particles interactions are significant, but the system is not crowded enough to cause slow dynamics or crystallisation. We set˜ = 1 (we expect results to depend weakly on this parameter) andγ = 10, which corresponds to strong damping. Previous work on biased ensembles focussed on the overdamped limit. Here we consider finite damping so that momenta are welldefined, this is convenient for the mechanical analysis. However, for the largeγ that we consider, we expect the physical behaviour to be similar to the overdamped limit. When presenting numerical results, we use nondimensional units in which l 0 , m, T are set to unity (so v 0 = 1 also).

B. Time-integrated clustering, and biased ensembles
We analyse dynamical trajectories over an observation time t obs . Our results are controlled by the hydrodynamic behavior of the system, so we define a dimensionless observation time τ obs = t obs /τ L where τ L = L 2 /D 0 is the hydrodynamic time scale, in which D 0 = T /(mγ) is the (bare) particle diffusivity. To measure clustering between particles, write r ij = |r i − r j | and define This function interpolates between a value of order unity when particles i and j are very close, and zero when their distance exceeds the cutoff distance l WCA . The total clustering within a trajectory is measured by integrating over time and summing over all pairs of particles: This quantity has been rendered dimensionless through normalisation by the hydrodynamic time scale. We consider biased ensembles of trajectories in which the average of any observable quantity A is given by where λ is a dimensionless biasing parameter and · 0 represents an average in the equilibrium state of the system. To understand the physical meaning of the bias, consider the average (dimensionless) clustering per particle in this ensemble: For large τ obs , the biased ensembles reproduce the mechanism for (rare) fluctuations of C(τ obs ) [29][30][31]34]. The value of λ encodes the size of the fluctuation, according to C(τ obs ) ≈ N τ obs C(λ). Positive λ corresponds to increased clustering. The factor of L 2 in (7) means that λ biases the hydrodynamic behavior, but has a weak effect on individual particle dynamics [31,40], see Sec. III D for further details. Comparing with previous work for ensembles biased by the dynamical activity [28,30], we expect large clustering to correspond to low dynamical activity [43].

A. Dynamical phase transition
We use transition path sampling (TPS) [44] to sample these biased ensembles. We take τ obs = 0.252, which is sufficiently large to reveal the phase transition, as we will see below. As in [40], we include auxiliary forces in the dynamics to improve sampling, similar to [36]. Fig  1 illustrates the behaviour of C(λ) as λ is increased from zero, for several system sizes, with fixedρ, τ obs . There is a change in behaviour for λ = λ c ≈ 450, which becomes increasingly pronounced for larger systems. This is a dynamical phase transition where symmetry is spontaneously broken (see Fig. 2, discussed below). The large numerical value of λ c is attributable to the small values of C in the unbiased system; the steepness of the WCA potential means that particle separations are never much less than l WCA . Fig. 2 illustrates the spontaneous symmetry breaking for λ > λ c : particles accumulate near one wall of the system. This can be be quantified by the average local density, defined as ρ λ (r) = i δ(r − r i ) λ . Accumulation at either wall is equally likely: as in equilibrium phase transitions, it is convenient to break the symmetry in the numerical computation, to obtain a clear signature of the symmetry-broken state. This is achieved by using asymmetric auxiliary forces in the TPS method, as discussed in Sec. III C, below.
The main features of this transition can be explained within macroscopic fluctuation theory (MFT) [9]. Following [28,40,45,46], one assumes that the clustering is a function of the local density, in which case the hydrodynamic behavior of biased ensembles can be predicted by minimization of a dynamical action that depends on the density alone. This indicates that (i) C(λ) should approach a scaling function in the limit of large system size N (consistent with Fig. 1); and (ii) the density profile should respond most strongly on the largest relevant length scale, which is the system size. This is consistent with Fig. 2 (and with Fig. 3, below). See Sec. III D for additional information on this hydrodynamic analysis.

B. Stress and force balance
So far, the analysis of this dynamical phase transition has similar features to previous work on simpler systems [25,40,45]. We now focus on mechanical properties of this biased ensemble, which have not been considered before, to our knowledge. We write the equation of motion (1) asṗ where is the body force on particle i. These are forces that do not conserve the total particle momentum P = Physically: the interparticle forces obey Newton's 3rd law so they can't change the total particle momentum, but the body forces can inject momentum from the particles' environment.] By analogy with ρ λ , define the average body-force density as F λ (r) = i f i δ(r − r i ) λ . Cartesian components of F λ are denoted by F α λ with α ∈ {x, y}, we use a similar notation for other vectors and tensors. The biased ensembles of trajectories that we consider are homogeneous in time, which means that they must have balanced mechanical forces. Specifically, body forces must be balanced by stress gradients: where Π λ is the local stress tensor (a position-dependent 2 × 2 matrix), and ∇ β indicates differentiation with respect to r β . We use the procedure of Irving and Kirkwood [41] to compute the stress in the biased ensemble, see Appendix A for details. Within that framework, we emphasise that (12) can be derived directly from the equations of motion of the system, there is no assumption that the system be near-equilibrium, only that it is stationary. Hence the formalism also applies in the biased ensemble.
The biased states retain translational invariance along the y-direction, so (12) reduces to For the unbiased dynamics, Eq (11) implies that F x 0 (x) = −V w (x)ρ 0 (x). Away from the wall, this quantity vanishes, so Π xx 0 is independent of x [by (13)]: in fact −Π xx 0 equals the pressure of the bulk fluid. This situation is shown in black in Fig. 3: the stress is constant in the bulk, while body forces generate gradients near the wall. The biased ensemble is time-reversal symmetric which means that p i δ(r−r i ) λ = 0, so the friction term in (11) does not contribute to the average body force. Hence, For λ > λ c , one sees from Fig. 3 that there is a stress gradient that extends into the bulk of the channel, so Eq. 13 requires a non-zero body force F x λ there. But V w = 0 in the bulk, so a non-zero body-force in (14) requires that the averaged Langevin noises must be non-zero, in the biased ensemble. This result may be counter-intuitive: it occurs because the biased ensemble changes the probabilities of different trajectories, so it also changes the probabilities of particular realisations of the thermal forces η i . This leads to a position-dependent body force within the biased ensemble. Indeed, some such force must be present, to maintain the asymmetric density profiles in Fig. 2.

C. Optimal-control representation and Doob's transform
In fact, these body forces are already familiar from theories of biased ensembles. As described above, the effects of the bias λ can be reproduced by adding control forces to the equations of motion, leading to an auxiliary model [31,33,34] as in Doob's transform. The equation of motion for the auxiliary model iṡ where φ i is a control force on particle i that depends (in general) on the positions and momenta of all particles in the system, andη i is a white noise with the same statistical properties as η i . Denote averages in the unbiased steady state of this auxiliary model by A aux . Then the Doob transform yields [33,34] that for a large class of observables A, in the limit of large t obs . The control forces φ i depend on λ, so A aux does too.
Note that the auxiliary model does not capture transient behavior of the biased ensemble for times t ≈ 0 and t ≈ t obs [30,34]. This restricts the the class of observables A to those for which transient effects have a negligible contribution for large t obs . All observables considered here are within this class.
In the steady state of the auxiliary model, the φ i act as body forces. The stress tensors of the biased ensemble and the auxiliary model are equal [by (16)] which means that their body forces must also be equal [by (12)]. Equating these forces yields (17) where we used that η i δ(r − r i ) aux = 0. The left hand side of (17) is the contribution of the φ i to the body force of the auxiliary model, and the right hand side is the contribution of the noise forces to the body force in the biased ensemble. We denote the left hand side of (17) by ρ λ (r)φ ave (r), so that φ ave (r) is the average control force on a particle at r.
Eq. (17) means that the control forces that appear in Doob's transform have a mechanical interpretation: they generate the biased noise forces identified in (14). To infer the control forces themselves, combine (12,14,17) to obtain That is, the (averaged) control force can be estimated from measurements of the density ρ λ and the stress Π λ , within the biased ensemble. Of course, the Doob force φ i may depend on the positions of all particles, while measuring the stress gradient only gives the average force φ ave (r). Eq. (17) could be generalised by adding additional delta functions inside the averages, such that Doob forces can be related to conditional averages of η i . However, we concentrate here on the information that is available from the stress tensor, since this is more easily measured than the noise forces themselves.
Integrating (18), we express the control force as φ x ave (x) = −V λ (x)/ρ λ (x) where we define the Doob stress, This is the part of the stress that originates from the noise forces. It is estimated numerically in Fig. 4. To approximate φ, we fit the stress as This function is zero near the walls, but negative in the bulk (corresponding to a force towards the left wall). Close inspection of Fig. 4 shows that the resulting fit to the stress is not perfect in the regions close to the wall, where layering takes place. However, it does capture the behaviour in bulk, and hence the hydrodynamic response to the bias.
A sample fit is shown in Fig. 4. Since the stress gradient occurs on the scale of the system size L, this control force is of order 1/L. This is another indication of the hydrodynamic response to the bias, which appears when a large number of particles each receives a weak bias, leading to a macroscopic response [31].
The estimated control force φ est can be inserted into (15) as φ i =xφ est (x i ), wherex is a unit vector in the x-direction. This gives an auxiliary model whose density profile is similar to the biased ensemble -this is not the exact Doob dynamics but it can be used within the TPS algorithm, to propose new trajectories that are more likely to be accepted. Using this method and systematically refining φ est greatly improves numerical performance of the algorithm, similar to [22,[36][37][38]40]. Exactly this method was used to obtain the data shown here -without such control forces, the simulations for larger systems and larger biases would have been computationally intractable. In this sense, insights from the mechanical analysis can be exploited to improve numerical methods.
Note that including control forces within TPS does not change the result of the numerical computation, which always samples the ensemble defined by (8), as long as t obs is large enough: this feature was discussed extensively in [36] in the context of population dynamics algorithms, and also in [40] for TPS. To understand it, one should think of TPS as a Monte Carlo (MC) procedure for trajectories of the system. Including auxiliary forces amounts to a different set of proposed MC updates, which does not change the ensemble being sampled [Eq. (8)], but may improve the acceptance rate. To verify that the target ensemble is not affected by the control forces, we checked that for λ < λ c then the sampled density profile is symmetric, even when asymmetric auxiliary forces are used. For λ > λ c , the auxiliary forces still do not affect the sampled distribution, but they do improve the acceptance rate within the broken-symmetry state. (Without such forces, the TPS algorithm tends to propose trajectory updates where the particles move away from the walls. These updates tend to be rejected, leading to inefficient sampling.)

D. Hydrodynamic response in biased ensemble
We noted in Sec. III A that the field λ biases the hydrodynamic behaviour of the system. This section gives some extra information on this point.
Based on the analogy between biased ensembles and thermodynamics [29][30][31], it is natural to define a measure of clustering that is extensive in space and time: Then the definition of the biased ensemble in (8) is In systems without hydrodynamic modes, it would be natural to use a definition similar to (22), but with λD 0 /L 2 replaced by an intensive field that is often denoted by s. Then one would keep s fixed in a joint limit L, t obs → ∞, which corresponds to applying a bias of order unity to each particle, even as the system size tends to infinity.
In the present context, we fix λ as L, t obs → ∞. Hence one sees from (22) that the intensive field s = λD 0 /L 2 is being reduced as the system size increases. The situation is familiar in systems with hydrodynamic modes [26,40,45], where weak biases on many particles can add up coherently to generate strong responses. This lends a degree of universality to such behaviour, which can be captured by theories like macroscopic fluctuation theory [9]. Similar mechanisms are at work in the system considered here.
We sketch the MFT analysis of the transition (following the analysis of [40] and [28,45] for one-dimensional systems): Rescale in hydrodynamic units asr = r/L andt = D 0 t/L 2 and assume that the clustering can be parameterised in terms of the local density as Q(t obs ) ≈ t obs 0 q(ρ(r,t))drdt for some "clustering density" q(ρ). Then the MFT action for the biased ensemble would be whereJ is the hydrodynamic current that obeys ∂tρ = −∇ ·J, and D, σ are the density-dependent diffusivity and mobility of MFT. (The effect of the walls would be incorporated through hard boundaries, after rescaling to the hydrodynamic scale.) Using time-translation invariance and time-reversal invariance of the biased ensemble, the time integral can be ignored and one hasJ = 0. Finally, the hydrodynamic density would be obtained by minimising the functional In the present context, the functional forms of D, σ, q are not known, so a quantitative analysis is not possible. However, it is clear on general grounds that for q > 0, positive λ will drive the system to an inhomogeneous state, similar to transitions in other contexts [40,45]. Fig. 5 shows results for N = 48, compared with the corresponding results for N = 24, as already shown in Figs. 2 and 3. Ignoring the layering effects near the walls (which are non-hydrodynamic in nature), and focussing on length scales of the order of the system size, one sees a semi-quantitative match of the density and stress gradients that develop as a function of λ. For full consistency with MFT, one would require that the density and stress would be scaling functions of x/L, in the limit of large systems N → ∞. In that case, the stress gradient in bulk would be O(1/L), as discussed in the main text.
In practice, the systems considered are far from the limit N → ∞, so microscopic length scales also affect the results. Still, the observed behaviour is consistent with the hydrodynamic theory.

IV. OUTLOOK
We summarize the results of this work: These twodimensional systems of interacting particles support a dynamical phase transition, associated with large deviations where the particles are clustered more than usual. At this transition, a symmetry is spontaneously broken, leading to particle aggregation at a wall. These transitions are naturally studied via biased ensembles of trajectories, which we have implemented numerically by TPS. In addition, while some similar behaviour is observed in lattice models [25,45], the fact that this system follows Newtonian dynamics means that the balance of mechanical forces can be investigated within the biased ensemble. In particular, the stress tensor of Irving and Kirkwood [41] can be computed within the biased ensemble, and yields useful information, even in systems far from equilibrium. (This property distinguishes the stress from objects like the thermodynamic free energy, which is no longer meaningful in biased ensembles.) This mechanical analysis shows that if particle motion depends on the bias parameter λ, this dependence must have its origin in some underlying forces. These forces appear as non-zero averages for the Langevin noises η i , and as control forces that appear in Doob's transform. Their presence can be inferred from the stress tensor. The physical picture is that the biased ensemble changes the probabilities of trajectories of the system, which changes the statistical properties of the noise. In particular, the bias tends to selects noise realisations that push particles towards one of the walls. This is the cause of the particle aggregation at the walls, and of the broken symmetry.
We emphasise that this mechanical perspective is not specific to the system analysed here, it is valid in any system where the stress tensor can be defined. In particular, we have considered inertial motion, but the same methodology is applicable in the overdamped limit, so large deviations of diffusions [34] can be analysed in the same way. For these reasons, we hope that this approach will provide new opportunities for understanding large deviations in physical systems, and for characterising them numerically.

ACKNOWLEDGMENTS
We thank Kris Thijssen, Mike Cates, and Tal Agranov for helpful discussions. We are grateful to the EPSRC for support through a studentship for JD (ref EP/N509620/1) and research funding to RLJ (ref EP/T031247/1). The data supporting this publication are available at https://doi.org/10.17863/CAM.86077.
Appendix A: Irving-Kirkwood stress 1. Definition, and derivation of (13) This section reviews the computation of the IK stress tensor [41]. For a modern presentation, we follow [47] (in particular, the Supplemental Material of that work).
As a preliminary for computing the stress, we consider the particle current. Define the empirical particle density and momentum density aŝ These quantities have implicit time-dependence via the positions and momenta. Then the time derivative of the particle density is where we used the chain rule andṙ i = p i /m (gradients ∇ are with respect to r). Since all particles have equal mass, the particle current may be identified aŝ (r) = (1/m)p(r), and one recognises (A2) as the continuity equation The manipulations so far are familiar from analysis of particle currents. The IK stress [41] is derived by applying similar procedures to the momentum density. The time derivative ofp is available by using the chain rule together with (10): denoting Cartesian components by Greek indices we obtain where F ij is the force on particle i from particle j, whose components are (A5) Eq. (A4) for the momentum is analogous to Eq. (A2) for the density. The analogue of the current in this case will be the stress tensor. To see this, use we use the insight of [41]: our discussion follows [47]. The key point is that (A4) can be rewritten as where we use (throughout this Appendix, but not in the main text) the convention of implicit summation of repeated Greek indices, and the IK stress tensor iŝ is a function which distributes unit weight over a straight line connecting r i and r j . The consistency of (A6,A7) with (A4) is demonstrated in App. A 2, below.
Observe that the WCA interaction is always repulsive [the object in square brackets in (A5) is positive], which means that diagonal elements ofΠ are always negative. This ensures that the (instantaneous local) mechanical pressure tr(−Π/2) is always positive.
We now derive (13) of the main text. Note that (A6) is an exact identity: it was derived from the equations of motion and it holds for every trajectory of the system. This means that it can be used to analyse force balance in the biased ensemble. In particular, we take the average of (A6) within the biased ensemble, to obtain The left hand side is zero (we again exclude transient regimes for t ≈ 0, t obs ). Define also the ensembleaveraged stress tensor as Then (A9) becomes 0 = F µ (r) + ∇ ν Π µν λ (r) which is (13) of the main text. [We used the definition of the bodyforce density F λ = i f i δ(r − r i ) λ .]

Validation of the IK stress formula (A7)
To see that (A6) is equivalent to (A4), we follow [47,48] and note the following property of h: where the gradient is with respect to r. To show this: Observe that for any path from a to b and any function Considering the straight line path from a to b and parameterising by λ, we define r λ = λb + (1 − λ)a. Then where the second line uses ∇f (x) = − f (r)∇δ(r−x)dr (the integral runs over all space, so an integration by parts shows that this applies for any function f ). Using again the definitions of the delta function and of h, this can be rearranged as: This holds for every function f , so the object in square brackets must vanish, this implies (A11). Now take (A11) with a, b = r j , r i and multiply by F µ ij , yielding ∇ ν [(r ν j − r ν i )F µ ij h(r; r i , r j )] = F µ ij δ(r − r i ) + F µ ji δ(r − r j ) (A14) (we used that F ij = −F ji ). Using this result with (A7), one obtains Interchanging dummy indices in the first sum and plugging into (A6), one recovers (A4). We note that the IK stress is not the only choice for a tensor that satisfies (A4), hence local stress is not uniquely defined, see for example [47]. However, the forces that we compute via stress gradients are unique. (The difficulty with the stress itself arises because interparticle forces have finite range; a unique stress can be defined mesoscopically by averaging over scales much larger than the range of the force.)

Stress measurements in the biased ensemble
This section details the numerical procedure that was used to estimate the stress.
For numerical measurements, (A7) is awkward because of the Dirac delta functions. A similar issue arises when estimating the density from simulations: the solution in that case is to measure the density as a histogram.
Specifically, one considers a region Ω of the system and definesN Ω = Ωρ (r)dr (A16) which is equal to the number of particles in Ω. This is easily computed in simulations, and one may estimate the average local density at r as N Ω /|Ω| where Ω is a region centred at r, and |Ω| is the volume of Ω. A similar method can be applied to measure the stress (see for example [49]), we definê where H Ω (r i , r j ) = Ω h(r; r i , r j )dr is the fraction of the line from r i to r j that passes through the volume Ω, and the notation i ∈ Ω indicates that we sum over particles whose positions r i are inside Ω. ThenΠ µν Ω is the local stress, averaged over Ω. This is called the volumeaveraged stress [49].
In practice, we divide the system into square boxes of size l 0 /8 and measure the stress in each box. Results for the stress tensor are then obtained by averaging along the y direction. The definition (A17) ensures that this procedure (measuring the stress in small boxes followed by averaging over the boxes) is equivalent to making the measurement on larger boxes from the outset.