Dynamical criticality in open systems: non-perturbative physics, microscopic origin and direct observation

Driven diffusive systems may undergo phase transitions to sustain atypical values of the current. This leads in some cases to symmetry-broken space-time trajectories which enhance the probability of such fluctuations. Here we shed light on both the macroscopic large deviation properties and the microscopic origin of such spontaneous symmetry breaking in the open weakly asymmetric exclusion process. By studying the joint fluctuations of the current and a collective order parameter, we uncover the full dynamical phase diagram for arbitrary boundary driving, which is reminiscent of a $\mathbb{Z}_2$ symmetry-breaking transition. The associated joint large deviation function becomes non-convex below the critical point, where a Maxwell-like violation of the additivity principle is observed. At the microscopic level, the dynamical phase transition is linked to an emerging degeneracy of the ground state of the microscopic generator, from which the optimal trajectories in the symmetry-broken phase follow. In addition, we observe this new symmetry-breaking phenomenon in extensive rare-event simulations, confirming our macroscopic and microscopic results.

Introduction.-The discovery of dynamical phase transitions (DPTs) in the fluctuations of nonequilibrium systems has attracted much attention in recent years . In contrast with standard critical phenomena [33,34], which occur at the configurational level, DPTs appear in trajectory space when conditioning the system to sustain an unlikely value of dynamical observables such as the time-integrated current [1,4,27,[35][36][37]. DPTs thus manifest as a peculiar change in the properties of trajectories responsible for such rare events, making these trajectories far more probable than anticipated due to the emergence of ordered structures such as traveling waves [2,9,11,21], condensates [3,12,29] or hyperuniform states [15,22,38]. In all these cases, the hallmark of the DPT is the appearance of a singularity in the so-called large deviation function (LDF), which controls the probability of fluctuations and plays the role of a thermodynamic potential for nonequilibrium systems [35,39,40]. DPTs play a key role to understand the physics of different systems, from glass formers [7,8,[41][42][43][44][45][46] to micromasers and superconducting transistors [47,48], and applications such as DPT-based quantum thermal switches [49][50][51]. Moreover, by making rare events typical with the use of Doob's transform [52][53][54] or optimal fields [40], one may exploit DPTs to engineer and control nonequilibrium systems with a desired statistics on demand [55].
In the context of diffusive systems, DPTs in current statistics have been throughly studied for periodic settings [2,4,9,11,24], in which the broken symmetry is time translational invariance, giving rise to a violation of the so-called additivity principle via travelingwave profiles [9,56]. Nevertheless, it has not been until very recently that other kind of symmetry-breaking scenarios (involving e.g. particle-hole symmetry) have been predicted for open systems [23], i.e. in contact with boundary reservoirs. In particular, a perturbative Landau theory restricted to zero or small boundary gradient has been recently put forward [23,26] which predicts 1 st -and 2 nd -order DPTs in some diffusive media. Key questions remain unanswered, however, such as the direct numerical observation of this DPT, its microscopic origin, the non-perturbative physics beyond the critical point, or its existence under strong boundary driving, the latter being one of the most challenging problems in nonequilibrium physics.
In this work we address these questions in a paradigmatic diffusive system, the open one-dimensional (1d) weakly asymmetric simple exclusion process (WASEP) [57,58]. In particular, by studying the joint fluctuations of the current q and a novel collective order parameter defined by total mass (m), we unveil analytically the full dynamical phase diagram for arbitrary boundary gradients, see Fig. 1. A DPT is observed at a critical current |q c | for any boundary driving symmetric around the density 1/2, i.e. for ρ R = 1 − ρ L (with ρ L and ρ R the left and right reservoir densities, respectively), where the joint mass- current LDF G(m, q) becomes non-convex (see Fig. 2). This signals the breaking of the particle-hole (PH) symmetry present in the governing action but no longer in the optimal trajectories associated to these atypical fluctuations: for |q| < |q c | coexisting low-and high-mass trajectories appear with broken PH-symmetry. An asymmetric boundary gradient favors one of the mass branches, deepening the associated minimum in G(m, q). Interestingly, in the regime where G(m, q) is non-convex, instantonlike time-dependent trajectories connecting the two local minima become optimal, demonstrating dynamical coexistence between the different symmetry-broken phases and signaling a violation of the additivity principle in open systems [2,9,11,56,[59][60][61][62]. A spectral analysis of the microscopic dynamical generator of the WASEP shows that the DPT is triggered by an emerging degeneracy of the associated ground state, from which one can compute the density profiles of the symmetry-broken phase. We provide also the first direct observation of this phenomenon through extensive rare-event simulations [63][64][65][66][67][68]. This work opens the door to studying DPTs in more complex scenarios, as e.g. open high-dimensional systems with multiple conservation laws, and represents a step forward in connecting current fluctuations with metastability and standard critical phenomena.
Model.-The WASEP belongs to a broad class of driven diffusive systems of fundamental interest [35,57,58]. Microscopically it consists of a 1d lattice of L sites, each of which may be empty or occupied by one particle at most. Particles hop randomly to empty neighboring left (right) sites at a rate 1 2 e −E/L ( 1 2 e E/L ), with E an external field. In addition, particles are injected and removed at the leftmost (rightmost) site at rates α and γ (δ and β), respectively, yielding in the diffusive limit boundary particle densities of ρ L = α/(α + γ) and ρ R = δ/(β + δ). At the mesoscopic level, driven diffusive systems like WASEP are characterized by a density field ρ(x, t) which obeys a stochastic equation [69] with D(ρ) and σ(ρ) the diffusivity and mobility coefficients, which for WASEP are D(ρ) = 1/2 and stands for the fluctuating current, and ξ is a Gaussian white noise, with ξ = 0 and ξ( , which accounts for microscopic fluctuations at the mesoscopic level. The density at the boundaries is fixed to ρ(0, t) = ρ L and ρ(1, t) = ρ R ∀t. DPT in the thermodynamics of currents.-When driven by E = 0 and/or ρ L = ρ R , the system relaxes to a nonequilibrium steady state characterized by an average current q and a non-trivial density profile ρ st (x) [70]. Moreover, we can associate to any trajectory {ρ(x, t), j(x, t)} τ 0 an empirical current q = τ −1 τ 0 dt 1 0 dx j(x, t). In the following we show how from the structure of the probability of this current, P (q), we can predict the existence of DPTs associated with  spontaneous symmetry breaking.
The probability P ({ρ, j} τ 0 ) of any trajectory can be computed from Eq. (1) via a path integral formalism [35,36,40], and scales in the large-size limit as P ({ρ, j} τ 0 ) ∼ exp{−L I τ [ρ, j]}, with an action [40] The probability P ({ρ, j} τ 0 ) represents the ensemble of space-time trajectories, from which one can obtain the statistics of any observable depending on {ρ, j} τ 0 . In particular the probability of a given current q can be obtained by minimizing the action functional (2) over all trajectories sustaining such current. This yields in the long-time limit P (q) ∼ exp{−τ LG(q)}, with G(q) = lim τ →∞ 1 τ min * {ρ,j} τ 0 I τ (ρ, j) the current LDF, and * meaning that the minimization must be compatible with the prescribed constraints (q, ρ L,R ). The optimal trajectories ρ q (x, t) and j q (x, t) solution of this variational problem are then those adopted by the system in order to maintain the current q over a long period of time, and turn out to be time-independent in many cases (a conjecture known as additivity principle [56]).
Just as in standard critical phenomena, the action (2) contains the symmetries which are eventually broken.
For WASEP with ρ R = 1 − ρ L it is easy to check that the action (2) is invariant under the transformation ρ → 1−ρ, x → 1 − x, referred to as PH symmetry (resulting from the symmetry of σ(ρ) around ρ = 1/2). The optimal density profile ρ q (x) typically inherits this PH symmetry, mapping onto itself under the above transformation. However, as detailed in the Supp. Mat. [70], for currents below a critical threshold (|q| ≤ |q c |) and large enough E, two different (but equally) optimal profiles ρ ± q (x) appear such that ρ ± q (x) → 1 − ρ ∓ q (1 − x), see inset to Fig. 1, giving rise to a second-order singularity in the current LDF. This spontaneous PH symmetry breaking can be easily understood [23,26] by noting that, in order to sustain a low-current fluctuation, the system can react by either crowding with particles hence hindering motion, or rather emptying the lattice to minimize particle flow. Both tendencies break the action PH symmetry, eventually triggering the DPT.
Order parameter fluctuations.-To better understand this DPT, we study the joint fluctuations of the current and an appropriate global order parameter for the transition, much in the spirit of the paradigmatic Ising model of standard critical behavior [33]. A natural choice for this order parameter is the total mass in the system, which clearly characterizes the DPT in this case but also in more complex scenarios. Indeed, as shown in Fig. 1, the typical mass during a current fluctuation, m q ≡ 1 0 dxρ q (x), exhibits a behavior strongly reminiscent of a standard Z 2 phase transition, capturing the PH symmetry breaking. Defining the empirical mass for a trajectory as m = τ −1 τ 0 dt 1 0 dxρ(x, t), the probability of observing a joint mass-current fluctuation for long times and large system sizes scales as P (m, q) ∼ exp{−τ LG(m, q)}, with G(m, q) = lim τ →∞ 1 τ min * {ρ,j} τ 0 I τ (ρ, j) being the mass-current LDF, such that G(q) = min m G(m, q) = G(m q , q). Within the additivity hypothesis [36,56,61] G(m, q) = min with the optimal profile ρ m,q (x) subject to the constraint m = 1 0 dxρ m,q (x) as well as to fixed boundary conditions. The mass constraint can be implemented using a Lagrange multiplier, and we solve analitycally the resulting problem in terms of elliptic integrals and Jacobi elliptic functions, see [70]. We note that the ρ m,q (x) so obtained can be classified attending to their extrema. Fig. 2 illustrates our results for strong boundary gradients, well beyond the linear nonequilibrium regime. In particular, for PH-symmetric boundaries (ρ R = 1 − ρ L ), the conditional mass-current LDF G(m|q) ≡ G(m, q) − G(q) exhibits a peculiar change of behavior at a critical current |q c |, see panel 2.a: while for |q| > |q c | the LDF G(m|q) displays a single minimum at m q = 1/2, with an associated PH-symmetric optimal profile ( Fig.  2.b), for |q| < |q c | two equivalent minima m ± q appear in G(m|q), each one associated with a PH-symmetry-broken optimal profile ρ ± q (x), see Fig. 2.c, such that . The emergence of this nonconvex regime in G(m|q) signals a 2 nd -order DPT to a PH-symmetry-broken dynamical phase. On the other hand, for PH-asymmetric boundaries (ρ R = 1 − ρ L ), the governing action (2) is no longer PH-symmetric: the asymmetry favors one of the mass branches and the associated G(m|q) displays a single global minimum ∀q and an unique optimal profile (see Fig. 2.d-f), explaining why no DPT is observed in this case [61]. Still, G(m|q) becomes non-convex for low enough currents, and for weak gradient asymmetry metastable-like local minima in G(m|q) may appear [70].
Maxwell construction and additivity violation.-A natural question is whether time-dependent optimal trajectories exist which improve the additivity principle minimizers. The emergence of a non-convex regime in G(m|q) for |q| < |q c | suggests a Maxwell-like instantonic solution in this region [26,32,71]. In particular, as we show in [70], for PH-symmetric boundaries, fixed |q| < |q c | and m ∈ (m − q , m + q ), a trajectory which jumps smoothly (in a finite time) , improves the additivity principle solution, yielding a straight Maxwell-like construction . This corresponds to a dynamical coexistence of the different symmetry-broken phases for |q| < |q c |, as expected for a 1 st -order DPT, see Fig. 1. Similar solutions exist for PH-asymmetric boundaries in regimes where G(m|q) is non-convex, leading to metastable dynamical coexistence, and we note that the role of the instanton around |q| ≈ |q c | can be affected by how the L → ∞ and τ → ∞ limits are taken [26].
Microscopic results: Spectral analysis.-Next we focus on the microscopic understanding of the symmetrybreaking DPT for current statistics. At the microscopic level, a configuration of the 1d WASEP is given by C = {n k } k=1,...,L , where n k = 0, 1 is the occupation number of the lattice's k th site. Within the quantum Hamiltonian formalism for the master equation [72], each configuration is represented as a vector in a Hilbert space, |C = L k=1 (n k , 1−n k ) T , with T denoting transposition. The complete information about the system is contained in a vector |P = (P (C 1 ), P (C 2 ), ...) T = i P (C i ) |C i , with P (C i ) representing the probabilities of the different configurations C i . This probability vector evolves according to the master equation ∂ t |P = W |P , where W defines the Markov generator of the dynamics. Such generator can be tilted W µ,λ [36,39,70] to bias the original stochastic dynamics in order to favor large (low) mass for µ < 0 (µ > 0) and large (low) currents for λ > 0 (λ < 0), with µ and λ the conjugate parameters to the microscopic mass and current observables, respectively. The connection between the biased dynamics and the large deviation properties of our system is established through the largest eigenvalue of W µ,λ [39,73]. Such eigenvalue, denoted by θ 0 (µ, λ), is nothing but the cu- We now consider exact numerical diagonalization of W µ,λ for a particular case of PH-symmetric boundaries and no mass bias (µ = 0). Fig. 3.a shows that the diffusively-scaled spectral gap, L 2 [θ 0 (0, λ) − θ 1 (0, λ)], with θ 1 (0, λ) the next-to-leading eigenvalue of W 0,λ , tends to zero as L increases in a region λ − [23,70]. This means that the 2 nd -order DPT in current statistics unveiled above at the macroscopic level corresponds to an emerging degeneracy of the ground state of W µ,λ (i.e. that corresponding to the leading eigenvalue), in which the sub-leading eigenvalue coalesces with the leading one. Moreover, by varying µ for λ = −E (equiv. q = 0) a remarkable 1 st -order-like behavior associated with a kink of θ 0 (µ, λ = −E) at µ = 0 is found, see inset to Fig. 3.b, consistent with the non-convex behavior of G(m|q = 0) found macroscopically and the associated dynamical coexistence of the two mass branches. Indeed, the numerical inverse Legendre transform of θ 0 (µ, λ = −E) converges to the convex envelope or Maxwell construction of the macroscopic prediction for G(m|q = 0), see Fig. 3.b.
The eigenspace associated to θ 0 (µ, λ) contains the microscopic information about the typical trajectories responsible for a given fluctuation (as parametrized by λ and µ). In this way, the emergence of a degeneracy as L increases points out to the appearance of two competing (symmetry-broken) states. For large but finite L, the spectral gap is small but non-zero and the eigenspace of θ 1 (µ, λ) defines a long-lived metastable state [74][75][76][77]. Using Doob's transform as a tool [53,55], one can show that any state in the degenerate (metastable) manifold is then given by a probability vector |P c MS =L 0 (|R 0 + c |R 1 ) [70]. Here |R i (|L i ) is the right (left) eigenvector associated with θ i (µ, λ) (i = 0, 1), andL 0 is a diagonal matrix whose elements (L 0 ) ii are the i th entries of |L 0 . Moreover, c ∈ [c 1 , c 2 ] is a constant with c 1 (c 2 ) the smallest (largest) entry of the vector L 1 |L −1 0 . Interestingly, our microscopic approach shows that the high-and low-mass states in the symmetry-broken phase then correspond to the states |P c1,c2

MS
, from which the average density profile in each phase can be computed. Fig. 4 shows the profiles so obtained from the exact numerical diagonalization of W µ,λ for two different gradients, and the convergence to the macroscopic prediction as L increases is clear.
Direct observation of the DPT.-So far, we have obtained clear indications of a symmetry-breaking DPT both from a macroscopic approach and a microscopic (spectral) analysis. The question remains as to whether this phenomenon is observable in simulations, which allow to reach larger system sizes. To address this we have performed extensive rare event simulations using the cloning Monte Carlo method [63,64,66,78] to study current statistics in the open 1d WASEP. Starting from random initial configurations, we have measured the optimal density profiles adopted by the system to sustain a highly atypical current, namely q = 0, using a population of 10 4 clones and L = 40. To capture the possible symmetry breaking, we average separately profiles with a total mass above and below 1/2. Fig. 4.a shows the result for ρ L = ρ R = 0.5, while Fig. 4.b displays data for ρ L = 0.8 and ρ R = 0.2 (in both cases E = 4). The measured high-and low-mass optimal profiles again converge towards the macroscopic predictions, strongly supporting our results on the PH-symmetry-breaking scenario.
Conclusions.-We have analyzed from a hydrodynamic, microscopic and computational point of view a 2 nd -order DPT in the current statistis of a paradigmatic driven diffusive system, the open 1d WASEP, unveiling the full dynamical phase diagram for arbitrary current fluctuations and boundary driving. For that we have investigated the joint fluctuations of the current and a collective order parameter, the total mass in the system, finding that the associated LDF becomes non-convex for low enough currents. Microscopically, we link the observed DPT with an emerging degeneracy of the ground state of the tilted dynamical generator, from which the macroscopic optimal profiles can be computed. Our predictions are confirmed by the observation of this DPT phenomenon for the first time in rare event simulations.
We thank Vivien Lecomte for insightful discussions. The research leading to these results has received funding from the EPSRC Grant No. EP/M014266/1 and the Spanish Ministry MINECO project FIS2017-84256-P. C.P.E. acknowledges the funding received from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Cofund Programme Athenea3I Grant Agreement No. 754446. C.P.E, P.I.H and J.P.G acknowledge as well the hospitality and support of the International Centre for Theoretical Sciences (ICTS) in Bangalore (India), where part of this work was developed during the program Large deviation theory in statistical physics: Recent advances and future challenges (Code: ICTS/Prog-ldt/2017/8). We are also grateful for access to the University of Nottingham High Performance Computing Facility, and for the computational resources and assistance provided by PRO-TEUS, the super-computing center of iC1 in Granada, Spain.