Sampling rare events across dynamical phase transitions

TL;DR

This paper reviews the use of cloning Monte Carlo methods to study dynamical phase transitions in stochastic lattice gases, highlighting techniques, challenges, and phenomenology of rare event characterization.

Contribution

It provides a detailed overview of applying cloning Monte Carlo techniques to analyze dynamical phase transitions in particle systems, including practical measurement tips.

Findings

Effective measurement of order parameters for DPTs

Finite-size effects impact DPT measurements

Cloning Monte Carlo methods enable direct observation of rare events

Abstract

Interacting particle systems with many degrees of freedom may undergo phase transitions to sustain atypical fluctuations of dynamical observables such as the current or the activity. This leads in some cases to symmetry-broken space-time trajectories which enhance the probability of such events due to the emergence of ordered structures. Despite their conceptual and practical importance, these dynamical phase transitions (DPTs) at the trajectory level are difficult to characterize due to the low probability of their occurrence. However, during the last decade advanced computational techniques have been developed to measure rare events in simulations of many-particle systems that allow for the first time the direct observation and characterization of these DPTs. Here we review the application of a particular rare-event simulation technique, based on cloning Monte Carlo methods, to characterize DPTs in paradigmatic stochastic lattice gases. In particular, we describe in detail some tricks and tips of the trade, paying special attention to the measurement of order parameters capturing the physics of the different DPTs, as well as to the finite-size effects (both in the system size and number of clones) that affect the measurements. Overall, we provide a consistent picture of the phenomenology associated with DPTs and their measurement.

Figures

Equations

$$P_t(Q)\asymp \exp [+tF(Q/t))],$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}\ln P_t(Q=\hat{q}t)=F(\hat{q}).$$

$$Z_t(\lambda )=\sum _Q{P}_t(Q){\text{e}}^{\lambda Q},$$

$$\theta (\lambda )=\underset{t\to \infty }{\lim }\frac{1}{t}\ln Z_t(\lambda ),$$

$$\theta (\lambda )=\underset{\hat{q}}{\max }\left[F(\hat{q})+\lambda \hat{q}\right].$$

$$j(x,t)=-D(\rho ){\partial }_x\rho (x,t)+\sigma (\rho )E+\xi (x,t).$$

$${\cal I}_{\tau}[\rho,j]=\int_{0}^{\tau}dt\int_{0}^{1}dx\frac{\displaystyle\Big{(}j+D(\rho)\partial_{x}\rho-\sigma(\rho)E\Big{)}^{2}}{\displaystyle 2\sigma(\rho)}\,.$$

$$G(q)=-\underset{\tau \to \infty }{\lim }\frac{1}{\tau }\underset{\{\rho ,j{\}}_0^{\tau }}{\min }{}^{\ast }{\mathcal{I}}_{\tau }(\rho ,j),$$

$$G(q)=-\underset{\rho (x)}{\min }{\int }_0^{1}dx\frac{{\left[q+D(\rho ){\partial }_x\rho -\sigma (\rho )E\right]}^2}{2\sigma (\rho )}.$$

$$\mu (\lambda )=\underset{q}{\max }\left[G(q)+\lambda q\right].$$

$${\mathbb{W}}_{C_0,C_1,\dots ,C_t}=W_{C_t{C}_{t-1}}\dots W_{C_1{C}_0}P(C_0),$$

$$P_t(Q)=\sum _{C_t..C_0}{W}_{C_t{C}_{t-1}}..W_{C_1{C}_0}P(C_0)\delta \left(Q-\sum _{k=0}^{t-1}{q}_{C_{k+1}{C}_k}\right),$$

$$Z_t(\lambda )=\sum _{C_t..C_0}{W}_{C_t{C}_{t-1}}..W_{C_1{C}_0}P(C_0){\text{e}}^{\lambda {\sum }_{k=0}^{t-1}{q}_{C_{k+1}{C}_k}}.$$

$${\tilde{W}}_{C^{\prime }C}\equiv {\text{e}}^{\lambda q_{C^{\prime }C}}{W}_{C^{\prime }C},$$

$$Z_t(\lambda )=\sum _{C_t\dots C_0}{\tilde{W}}_{C_t{C}_{t-1}}\dots {\tilde{W}}_{C_1{C}_0}P(C_0).$$

$$Y_C\equiv \sum _{C^{\prime }}{\tilde{W}}_{C^{\prime }C},$$

$$W_{C^{\prime }C}^{\prime }\equiv \frac{1}{Y_C}{\tilde{W}}_{C^{\prime }C}=\frac{{\text{e}}^{\lambda q_{C^{\prime }C}}}{Y_C}{W}_{C^{\prime }C},$$

$$Z_t(\lambda )=\sum _{C_t\dots C_0}{W}_{C_t{C}_{t-1}}^{\prime }{Y}_{C_{t-1}}\dots W_{C_1{C}_0}^{\prime }{Y}_{C_0}P(C_0).$$

$$\mathcal{N}(C_0,\dots ,C_t;t)=W_{C_t{C}_{t-1}}^{\prime }{Y}_{C_{t-1}}\mathcal{N}(C_0,\dots ,C_{t-1};t-1),$$

$$\mathcal{N}(C_0,\dots ,C_t;t)=W_{C_t{C}_{t-1}}^{\prime }{Y}_{C_{t-1}}\dots W_{C_1{C}_0}^{\prime }{Y}_{C_0}\mathcal{N}(C_0;0),$$

$$⟨\mathcal{N}(t)⟩$$

$$Z_t(\lambda )=\frac{⟨\mathcal{N}(t)⟩}{N_c}=\frac{⟨\mathcal{N}(t)⟩}{⟨\mathcal{N}(t-1)⟩}\frac{⟨\mathcal{N}(t-1)⟩}{⟨\mathcal{N}(t-2)⟩}\dots \frac{⟨\mathcal{N}(1)⟩}{N_c}.$$

$$\theta (\lambda )=\underset{t\to \infty }{\lim }\frac{1}{t}\sum _{k=1}^t\ln \left(\frac{⟨\mathcal{N}(k)⟩}{⟨\mathcal{N}(k-1)⟩}\right),$$

$$\theta (\lambda )\approx -\frac{1}{t}\sum _{k=1}^t\ln X_k\text{for }t\gg 1,$$

$$Y_t^{\text{max}}\equiv \max (Y_{C_t[1]},\dots ,Y_{C_t[N_c]}),$$

$$P_{\lambda ,<}^{\text{max}}(y)=[1-P_{\lambda }^{>}(y){]}^{N_c}.$$

$$P_{\lambda ,<}^{\text{max}}(y)\simeq \exp [-N_c{P}_{\lambda }^{>}(y)].$$

$$P_{\lambda }^{>}[y_{\lambda }^{\ast }(p)]=\frac{1}{N_c}\ln \left(\frac{1}{p}\right).$$

$$P_{{\lambda }_c(p)}^{>}(N_c)=\frac{1}{N_c}\ln \left(\frac{1}{p}\right).$$

$$P_t(C;Q)=\sum _{C^{\prime }}{W}_{C{C}^{\prime }}{P}_{t-1}(C^{\prime };Q-q_{C{C}^{\prime }}).$$

Full text

Large dynamical fluctuations are realizations of the dynamics sustained during a long period of time which deviate very much from their average value. Despite being very unlikely to occur, these fluctuations appear in many different systems carrying a large impact. Examples range from oceanic rogue waves, or chemical reaction kinetics to climate changes or stock market crashes. In the context of nonequilibrium systems, where we cannot derive the macroscopic properties from the Boltzmann-Gibbs distribution, the study of these fluctuations has led to important advances such as the fluctuation theorems and the formulation of a macroscopic fluctuation theory for driven diffusive systems. However, one of the main challenges to make further progress is precisely that large fluctuations are difficult to be observed as their probability is extremely low. Thus, much effort is currently devoted to implement efficient algorithms allowing for the measurement of such atypical events. In this paper we review the application of a computational technique, based on population dynamics, to measure large fluctuations of the current sustained in driven diffusive systems. We focus our attention on the application of the algorithm to the study of the so-called dynamical phase transitions, which appear as a change in the trajectories of the system in order to maximize the probability of sustaining a large fluctuation. Transitions of this kind are unveiled for several paradigmatic stochastic lattice gases after introducing the correct order parameters.

I Introduction

Phase transitions appear ubiquitously in nature, from cosmological scales to the quantum world of elementary particles. Consequently, the theory of critical phenomena has become one of the cornerstones of modern theoretical physics Binney et al. (1992); Zinn-Justin (2002). According to the modern classification (which is based on the classical Ehrenfest scheme) most phase transitions can be broadly divided into two categories: first-order or discontinuous, and second-order or continuous Binney et al. (1992); Zinn-Justin (2002), although other types do exist which challenge this classification. In a typical second-order phase transition, some type of order emerges continuously at a critical point of a control parameter. This peculiar change is captured by a so-called order parameter, and typically signals the spontaneous breaking of a symmetry and an associated non-analyticity of the relevant thermodynamic potential. Conversely, first-order transitions are characterized by an abrupt jump in the order parameter related to a kink in the thermodynamic potential, which leads to a coexistence between different phases Binney et al. (1992); Zinn-Justin (2002).

In recent years these ideas have been generalized to the realm of fluctuations, where dynamical phase transitions have been identified in many different systems Bertini et al. (2005); Bodineau and Derrida (2005); Harris et al. (2005); Bertini et al. (2006); Bodineau and Derrida (2007); Lecomte et al. (2007a, b); Garrahan et al. (2007); Bodineau et al. (2008); Garrahan et al. (2009); Hedges et al. (2009); Chandler and Garrahan (2010); Garrahan and Lesanovsky (2010); Hurtado and Garrido (2011); Garrahan et al. (2011); Pitard et al. (2011); Genway et al. (2012); Ates et al. (2012); Speck et al. (2012); Pérez-Espigares et al. (2013); Harris et al. (2013); Villavicencio-Sanchez et al. (2014); Lesanovsky et al. (2013); Hurtado et al. (2014); Vaikuntanathan et al. (2014); Manzano and Hurtado (2014); Jack et al. (2015); Shpielberg and Akkermans (2016); Zarfaty and Meerson ; Nyawo and Touchette (2016); Manzano and Kyoseva (2016); Lazarescu (2017); Brandner et al. (2017); Karevski and Schütz (2017); Carollo et al. (2017); Baek et al. (2017); Tizón-Escamilla et al. (2017a); Shpielberg (2017); Pinchaipat et al. (2017); Abou et al. (2018); Manzano and Hurtado (2018); Baek et al. (2018); Shpielberg et al. (2018); Pérez-Espigares et al. (2018a, b); Chleboun et al. (2018); Klymko et al. (2018); Whitelam (2018); Vroylandt and Verley (2018); Rotondo et al. (2018); Buča et al. (2019); Doyon and J. (2019). But what is dynamical about dynamical phase transitions? In contrast to standard critical phenomena, which occur at the configurational level when varying a control parameter such as temperature or magnetic field, DPTs appear in trajectory space when conditioning a system of interest to sustain an unlikely value of dynamical observables such as the time-integrated current or the activity (which are key magnitudes when studying respectively nonequilibrium systems and amorphous solids). The different dynamical phases that appear correspond to different types of trajectories adopted by the system during these rare events. Interestingly, some dynamical phases and their corresponding trajectories turn out to be far more probable than anticipated due to the emergence of ordered structures such as traveling waves Bodineau and Derrida (2005); Hurtado and Garrido (2011); Pérez-Espigares et al. (2013); Karevski and Schütz (2017), condensates Harris et al. (2005, 2013); Chleboun et al. (2018) or hyperuniform states Jack et al. (2015); Carollo et al. (2017, 2018a). Another hallmark of a DPT is the appearance of non-analyticities and Lee-Yang singularities Yang and Lee (1952); Arndt (2000); Blythe and Evans (2002); Dammer et al. (2002); Blythe and Evans (2003); Flindt and Garrahan (2013); Hickey et al. (2014); Brandner et al. (2017) in the so-called large deviation function (LDF) which controls the probability of fluctuations. This is a finding of crucial importance particularly in nonequilibrium physics, as these LDFs play a role akin to the equilibrium thermodynamic potentials for nonequilibrium systems, where no bottom-up approach exists yet connecting microscopic dynamics with macroscopic properties Bertini et al. (2015); Derrida ; Barré et al. (2018).

Symmetry-breaking DPTs are particularly interesting, and their analogy with standard critical phenomena is intriguing. In a standard second-order critical point, a (continuous or discrete) symmetry of the governing action is eventually broken, meaning that the system ground state beyond the critical point has less symmetries than the original action. Symmetry is however recovered by the appearance of different (symmetry-broken) ground states, which map onto each other under the symmetry operator. Remarkably, a similar picture arises in DPTs, but this time at the level of trajectories. To better understand this point, note that the probability of observing a given long-time fluctuation of a dynamical observable is dominated by the probability of the most probable trajectory (or optimal path) leading to such fluctuation Hurtado et al. (2014); Bertini et al. (2015). This defines a sort of dynamical ground state for each fluctuation, i.e. its optimal path, whose properties are crucial to shed light on the physics of the problem of interest Pérez-Espigares et al. (2015, 2016); Tizón-Escamilla et al. (2017b). The action-like functional describing the statistical weight of paths in phase space may have some symmetries (as e.g. time-translation invariance, particle-hole exchange, etc.) which are typically inherited by the associated optimal paths. However, at a second-order-like DPT, the symmetry of the trajectory action is broken: optimal paths do not share the symmetry of their action, but symmetry is restored by the appearance of degenerate optimal trajectories linked by the symmetry transformation. This analogy can be further exploited to obtain deep insights into symmetry-breaking DPTs.

In addition to their conceptual importance, DPTs play also a key role to understand the physics of different systems, from glass formers Garrahan et al. (2007, 2009); Hedges et al. (2009); Chandler and Garrahan (2010); Pitard et al. (2011); Speck et al. (2012); Pinchaipat et al. (2017); Abou et al. (2018) to superconducting transistors and micromasers Garrahan et al. (2011); Genway et al. (2012). There have been also recent applications of DPTs to design quantum thermal switches Manzano and Hurtado (2014); Manzano and Kyoseva (2016); Manzano and Hurtado (2018), i.e. quantum devices where the heat current flowing between hot and cold reservoirs can be completely blocked, modulated or turned on at will. Furthermore, another interesting possibility opens up by noting that rare events can be turned to typical with the use of Doob’s h-transform Doob (1957); Jack and Sollich (2010); Chetrite and Touchette (2015a, b) or external fields with optimal dissipation Bertini et al. (2015). This can be then used to exploit existing DPTs to engineer and control complex systems with a desired statistics on demand Carollo et al. (2018b), a possibility which is being currently explored.

Despite their relevance, observing and characterizing DPTs is a challenging task, the reason being that the spontaneous emergence of large (rare) fluctuations in many-body systems is generally unlikely. However, during the last few years, two new powerful and general methods have appeared to investigate fluctuating behavior in many-particle systems that are broadening our understanding of DPTs. On one hand, at the theoretical level, a macroscopic fluctuation theory (MFT) has been formulated Bertini et al. (2015) which offers variational equations to rationalize dynamical fluctuations and the associated LDFs in interacting many-particle systems arbitrarily far from equilibrium, starting from their fluctuating hydrodynamic description and a few transport coefficients. This deep and rich theoretical scheme won’t be the focus of this paper, though we will refer to some of its predictions; reviews on this formalism can be found elsewhere Bertini et al. (2015); Derrida ; Hurtado et al. (2014). A second tool which has reinvigorated the study of DPTs is the development of advanced computational methods to directly measure LDFs and the associated optimal paths in simulations of interacting many-particle systems Giardinà et al. (2006); Lecomte and Tailleur ; Tailleur and Lecomte (2009); Giardinà et al. (2011). These numerical methods amount to modifying the system dynamics so that the rare events responsible for a large deviation become no longer rare, and involve the simultaneous evolution of multiple copies or clones of the system of interest, which replicate or die in time according to their statistical weight, a technique based on the Diffusion Monte Carlo method of quantum mechanics Anderson (1975). The application of these new tools to simple models, particularly stochastic lattice gases, is providing intriguing evidences of the existence of rich and fundamental structures in the fluctuating behavior of nonequilibrium systems, which emerge mainly via DPTs, crucial to crack this long-unsolved problem.

The aim of this paper consists in reviewing the application of the cloning Monte Carlo method to understand the physics behind a number of dynamical phase transitions of theoretical interest. In particular, we will describe in detail some tricks of the trade, with an emphasis on the definition and measurement of order parameters capturing the physics of the different DPTs, the characterization of the optimal paths responsible for a fluctuation, as well as the finite-size effects (both in the system size and number of clones) that affect the measurements. In order to do so, we first provide a brief description of the cloning Monte Carlo method in Section III, both in its discrete- and continuous-time versions, together with an analysis of the effect of the finite number of clones on the large deviation function estimators. Once the main computational tools have been introduced, we set out to describe some intriguing results for DPTs in the current statistics of two paradigmatic models of transport, namely the Kipnis-Marchioro-Presutti (KMP) model of heat transport Kipnis et al. (1982) and the weakly asymmetric simple exclusion process (WASEP) Derrida (1998). In Section IV we study the spontaneous breaking of time-translation symmetry at the trajectory level in periodic systems. We do so in the $1d$ KMP model, where a DPT into a dynamical phase dominated by ballistic energy packets has been found, and in $1d$ WASEP, where jammed density-wave trajectories dominate low-current fluctuations. In Section V we explore numerically the important role of dimensionality on dynamical phase transitions by investigating vector current statistics in $2d$ WASEP. Interestingly, the complex interplay among an external field, the possible system anisotropy, and vector currents in $d>1$ leads to a rich phase diagram at the fluctuating level, with different symmetry-broken fluctuation phases separated by lines of first- and second-order DPTs. This remarkable competition between different dynamical phases is due to the appearance of a structured vector field coupled to the current, a key feature of high-dimensional (realistic) systems. Section VI is devoted to study a different type of symmetry-breaking phenomenon at the trajectory level which appears in open systems, i.e. coupled to boundary reservoirs which may drive the system out of equilibrium by imposing an external gradient (of e.g. density or temperature). In this case the symmetry that is broken at the DPT is the particle-hole symmetry –a $\mathbb{Z}_{2}$ discrete symmetry–, and we show that the transition persists in the presence of arbitrarily strong (but symmetric) boundary gradients. Finally, Section VII discusses the results presented in this review from a general point of view, providing also some outlook on the work that remains to be done in a near future.

II The statistical physics of trajectories

As mentioned earlier, we will be interested in this paper on DPTs emerging in the statistics of a dynamical observable such as the space&time-averaged current. However, the following large deviation formalism applies as well to any time-integrated observable. In order not to clutter our notation we particularize our discussion to one-dimensional ($1d$) systems (unless otherwise stated), though extensions to arbitrary dimension and vector currents are straightforward Hurtado et al. (2014). With this aim in mind, we hence consider the statistical physics of an ensemble of trajectories conditioned to a given total current $Q$ integrated over a long time $t$. This trajectory ensemble is fully characterized by the probability $P_{t}(Q)$ of all trajectories of duration $t$ with total current $Q$. In most cases of interest, this probability can be shown to obey a large-deviation principle for long times $t$ Bertini et al. (2015); Derrida ; Hurtado et al. (2014), i.e. $P_{t}(Q)$ scales in this limit as

[TABLE]

where the symbol ”$\asymp$” represents asymptotic logarithmic equality, i.e.

[TABLE]

The function $F(\hat{q})$ in Eq. (1) above, with $\hat{q}=Q/t$, defines the large deviation function (LDF) of the current. This LDF is a measure of the (exponential) rate at which the probability of observing an empirical current $\hat{q}$ –appreciably different from its steady-state value $\langle\hat{q}\rangle$– decays as $t$ increases. Note that this implies that $F(\langle\hat{q}\rangle)=0$. In the spirit of ensemble theory of equilibrium statistical mechanics, one can also characterize the system in terms of a dynamical partition function

[TABLE]

or equivalently by the associated dynamical free energy (dFE)

[TABLE]

which is nothing but the Legendre transform of the current LDF, namely

[TABLE]

The intensive parameter $\lambda$ is conjugated to the time-extensive current $Q$. This relation is equivalent to the connection between temperature and energy in equilibrium systems. However, and unlike temperature, the parameter $\lambda$ is non-physical and cannot be directly manipulated is experiments, a main difficulty when studying DPTs which can be however circumvented using the active interpretation of fluctuation formulas Bertini et al. (2015). In any case, fixing $\lambda$ to a constant value is equivalent to conditioning the system of interest to have a time-averaged (intensive) current $\hat{q}_{\lambda}\equiv Q_{\lambda}/t=\partial_{\lambda}\theta(\lambda)$, so by varying $\lambda$ one may change the associated current and therefore move from one dynamical phase to another.

Most of the systems whose dynamical fluctuations we are interested in can be described at the mesoscopic level by a locally-conserved density field $\rho(x,t)$ which evolves in time according to a fluctuating hydrodynamic equation. This can be seen as a continuity equation $\partial_{t}\rho+\partial_{x}j=0$ coupling the local density field $\rho(x,t)$ with a fluctuating local current $j(x,t)$. This current field typically obeys Fick’s (or Fourier’s) law, and includes a noise term which captures all the fast degrees of freedom which have been re-summed in the coarse-graining procedure leading to this mesoscopic description, namely

[TABLE]

Here $D(\rho)$ and $\sigma(\rho)$ are the diffusivity and mobility transport coefficients, respectively, and $E$ is a possible external field applied on the system of interest. The stochastic field $\xi$ is a Gaussian white noise, with $\langle\xi\rangle=0$ and $\langle\xi(x,t)\xi(x^{\prime},t^{\prime})\rangle=L^{-1}\sigma(\rho)\delta(x-x^{\prime})\delta(t-t^{\prime})$, where $L$ is the size of the system. Moreover, this fluctuating hydrodynamic description must be supplemented by appropriate boundary conditions, which can be either periodic or open (see below). Interestingly, we can associate with any trajectory $\{\rho(x,t),j(x,t)\}_{0}^{\tau}$ of duration $\tau$ in mesoscopic phase space an empirical space&time-averaged current $q=\tau^{-1}\int_{0}^{\tau}dt\int_{0}^{1}dx~{}j(x,t)$. Due to the diffusive scaling when going from the microscopic description to a mesoscopic stochastic field theory, the relation between the macroscopic current $q$ and the microscopic one $\hat{q}$ can be shown to be $q=L\hat{q}$ Derrida . The probability $P(\{\rho,j\}_{0}^{\tau})$ of any trajectory can be computed using a path integral formalism Hurtado et al. (2014); Bertini et al. (2015); Derrida , and scales in the large-size limit as $P(\{\rho,j\}_{0}^{\tau})\asymp\exp\{-L\,{\cal I}_{\tau}[\rho,j]\}$, with an action Bertini et al. (2015)

[TABLE]

This probability measure represents the ensemble of space-time trajectories at this mesoscopic level of description. The probability of a given current $q$ can be now obtained by minimizing the action functional (7) over all trajectories sustaining such current, leading in the long-time limit to $P(q)\asymp e^{\tau LG(q)}$, with a current LDF

[TABLE]

where ∗ means that the minimization procedure must be compatible with the prescribed constraints ($q$, boundary conditions, $\partial_{t}\rho+\partial_{x}j=0$, etc). The optimal trajectories $\rho_{q}(x,t)$ and $j_{q}(x,t)$ solution of this variational problem define the path the system follows to sustain a current $q$ over a long period of time, and turn out to be time-independent in many cases (a conjecture known as additivity principle Bodineau and Derrida (2004)). In this case, the current LDF simplifies to

[TABLE]

As we will analyze below, sometimes this time-independent solution becomes unstable for large enough current fluctuations under periodic boundary conditions, a DPT leading to a more complex (and now time-dependent) traveling wave optimal trajectory $\rho_{q}(x,t)=\omega_{q}(x-vt)$ characterized by a non-trivial velocity $v$. More details on this MFT problem can be found elsewhere Bodineau and Derrida (2005); Hurtado and Garrido (2011); Pérez-Espigares et al. (2013); Hurtado et al. (2014).

Note that both in the general case and in the time-independent approximation, equivalent variational problems can be formulated for the dynamical free energy $\mu(\lambda)$ Bertini et al. (2015); Derrida ; Hurtado et al. (2014), which reads

[TABLE]

The relation between the microscopic and macroscopic LDFs and their corresponding dynamical free energies are $F(\hat{q})=L^{-1}G(q)$ and $\theta(\lambda)=L^{-1}\mu(\lambda)$ in $d=1$ Hurtado et al. (2014); Derrida . Dynamical phase transitions correspond to singularities in the dFE $\mu(\lambda)$ or equivalently in the current LDF $G(q)$. These singularities (which typically are first- or second-order, as explained in the previous section) are accompanied by peculiar changes in the most probable trajectories or optimal paths associated with these fluctuations. For instance, a broad current interval may be dominated by time-independent optimal trajectories, i.e. by the additivity principle solution Bodineau and Derrida (2004); Hurtado and Garrido (2009a), while in some cases it is known that this additivity conjecture breaks down at some critical current, beyond which time-dependent traveling-wave-like optimal paths dominate the variational problem. This singular change can be rationalized as a second-order DPT where time-translation symmetry is broken Bodineau and Derrida (2005); Bertini et al. (2005, 2006); Hurtado and Garrido (2011); Pérez-Espigares et al. (2013); Hurtado et al. (2014). Similarly, systems with particle-hole symmetry may exhibit regimes of current fluctuations where the dominant trajectory (or dynamical ground state as termed in the introduction) breaks such particle-hole invariance Baek et al. (2017, 2018); Pérez-Espigares et al. (2018b).

III Sampling rare events with cloning Monte Carlo

We now turn our attention to the measurement and characterization of rare events in simulations of stochastic many-particle systems. This is of course a difficult task as these rare events are exponentially unlikely (as measured by their associated large deviation function), and hence typically hard to observe in experiments or simulations. Different computational strategies have been proposed to solve this issue, ranging from transition path sampling techniques Dellago et al. (2002); Bolhuis et al. (2002) or density matrix renormalization group methods Gorissen et al. (2009); Gorissen and Vanderzande (2012); Gorissen et al. (2012) to cloning Monte Carlo simulations Giardinà et al. (2006); Lecomte and Tailleur ; Tailleur and Lecomte (2009); Giardinà et al. (2011). In this paper we will focus on the latter method, which is best suited to measure in nonequilibrium situations the probability of a large deviation for time-extensive observables such as the current or the activity. We will describe two different versions of this general method, namely as applied to stochastic many-particle systems evolving with discrete- and continuous-time dynamics. Computing the probability of a rare event is equivalent to calculating the leading spectral properties of a tilted or deformed stochastic matrix which no longer conserves probability. In this way, the key idea consists in reading this tilted matrix as a new dynamical generator which makes the rare event typical, interpreting the non-conservation of probability as an effective population dynamics where different clones or copies of the system reproduce and die according to their statistical weights in the tilted dynamics.

Before proceeding, we note that the material in this section has been described in more detail elsewhere Giardinà et al. (2006); Lecomte and Tailleur ; Tailleur and Lecomte (2009); Hurtado and Garrido (2009b, 2010); Giardinà et al. (2011); Hurtado et al. (2014), here we summarize the main ideas and results for clarity and completeness purposes.

III.1 Discrete-time cloning algorithm

We hence consider a stochastic many-particle Markovian system in discrete time whose microscopic configuration at time $t$ is given by $C_{t}$. We will be here interested in the statistics of trajectories, which are nothing but sequences of configurations $\{C_{0},C_{1},\ldots,C_{t}\}$. If $W_{C^{\prime}C}$ is the transition probability from configuration $C$ to configuration $C^{\prime}$ in the stochastic model of interest, then by virtue of the Markov property the probability of a path is given by

[TABLE]

where $P(C_{0})$ is an initial distribution for the first configuration. Associated with each jump in configuration space, $C\to C^{\prime}$, one can define the elementary current $q_{C^{\prime}C}$ involved in this microscopic transition. The probability of observing a total time-integrated current $Q$ after a time $t$ is now given by the sum of the probabilities of all paths compatible with such current, and hence can be formally written as

[TABLE]

where the Dirac delta-function implements the constraint on the total current. Working with a global constraint on the total current as in Eq. (12) is typically difficult, so it seems convenient to change ensemble in the spirit of equilibrium statistical mechanics and work with the moment-generating function of the current distribution, i.e. the dynamical partition function $Z_{t}(\lambda)=\sum_{Q}P_{t}(Q)e^{\lambda Q}$ of previous section, see Eq. (3). Using (12) in the definition of $Z_{t}(\lambda)$, we find that

[TABLE]

This expression suggests defining a modified dynamics

[TABLE]

so the dynamical partition function can be now written as

[TABLE]

However, the transition matrix with elements $\tilde{W}_{C^{\prime}C}$ does not define a proper stochastic matrix, as the modified dynamics is not normalized, $\sum_{C^{\prime}}\tilde{W}_{C^{\prime}C}\neq 1$. Introducing now the exit rates

[TABLE]

we can define now a bona fide normalized dynamics

[TABLE]

which is now correctly normalized. The dynamical partition function now reads

[TABLE]

This sum over paths can be read from a computational point of view in terms of a population dynamics, which combines a cloning or replication step proportional to the exit rate of a configuration, $Y_{C}$, followed by a configuration jump governed by the modified normalized dynamics $W^{\prime}_{C^{\prime}C}$. The cloning step can be realized by substituting a particular configuration $C$ by an integer number of copies (0, 1, 2, …) with average $Y_{C}$, which then evolve independently according to the stochastic matrix $W^{\prime}_{C^{\prime}C}$. This population dynamics typically gives rise to an exponential increase or reduction in the number of copies of the system. Indeed, for a given trajectory $\{C_{0},\ldots,C_{t}\}$, the number ${\cal N}(C_{0},\ldots,C_{t};t)$ of copies at time $t$ obeys a simple recurrence relation

[TABLE]

so iterating we arrive at

[TABLE]

where ${\cal N}(C_{0};0)=N_{c}P(C_{0})$, with $N_{c}$ the initial total number of copies. The average $\langle{\cal N}(t)\rangle$ of the total number of clones at time $t$ is obtained by summing over all trajectories of such duration, and is proportional to the dynamical partition function $Z_{t}(\lambda)$, see Eq. (18) above,

[TABLE]

where we have used in the last step the exponential scaling of the dynamical partition function at long times, see Eq. (4), which demonstrates the exponential scaling of the population of clones as time evolves. In this way,

[TABLE]

The last equality trivially splits the total growth of the population into the different increments at each time step. This allows us to deal numerically with the exponential explosion (or implosion) of the total population. Indeed, we find for the dynamical free energy, see Eq. (4),

[TABLE]

with the definition $\langle{\cal N}(0)\rangle\equiv N_{c}$.

We are now in position to exploit the previous ideas to sample the sum over paths that define the dynamical partition function $Z_{t}(\lambda)$ using Monte Carlo techniques, in particular inspired by the Diffusion Monte Carlo method of quantum mechanics Anderson (1975). This sampling can be realized by considering a large ensemble of $N_{c}\gg 1$ copies of the system of interest (also dubbed clones, and hence the name of cloning Monte Carlo methods) which evolve sequentially according to the following steps Giardinà et al. (2006, 2011):

1. (a)

Each copy evolves independently as dictated by the modified normalized dynamics $W^{\prime}_{C^{\prime}C}$, which favors certain local jumps depending on the current $q_{C^{\prime}C}$ involved, see Eq. (17). 2. (b)

Each copy $m\in[1,N_{c}]$ (in configuration $C_{t}[m]$ at time $t$) is cloned with rate $Y_{C_{t}[m]}$, i.e. if $\lfloor x\rfloor$ represents the integer part of $x$, we generate a number $K_{C_{t}[m]}=\lfloor Y_{C_{t}[m]}\rfloor+1$ of identical clones with probability $Y_{C_{t}[m]}-\lfloor Y_{C_{t}[m]}\rfloor$, or $K_{C_{t}[m]}=\lfloor Y_{C_{t}[m]}\rfloor$ otherwise. This cloning step includes the possibility of a copy leaving no offspring if $K_{C_{t}[m]}=0$.

This procedure gives rise to a total of $N^{\prime}_{c}(t)=\sum_{m=1}^{N_{c}}K_{C_{t}[m]}$ copies after cloning all of the original $N_{c}$ copies. To avoid numerical issues with the exponential explosion or implosion of the total population, we add a third, population-control step to the method described above:

1. (c)

Once all copies have been evolved and replicated as dictated by their exit rates, the total number of copies $N^{\prime}_{c}(t)$ is sent back to $N_{c}$ by an uniform cloning rate $X_{t}=N_{c}/N^{\prime}_{c}(t)$.

Fig. 1 sketches the cloning Monte Carlo algorithm. When the total number of clones $N_{c}$ is large enough (and this is a key point in the algorithm), the global cloning factors $X_{t}$ are a good estimator of the average population increments in each time step, $X_{t}\approx\langle{\cal N}(t-1)\rangle/\langle{\cal N}(t)\rangle$, and hence we obtain an estimator of the dynamical free energy

[TABLE]

an expression which is expected to be exact in the limit of infinite number of clones, $N_{c}\to\infty$.

III.2 Effects of a finite population of clones on the estimation of dynamical free energies

Of course, in a real simulation we have at our disposal a hopefully large but still finite number of copies $N_{c}$, and this limitation may drive the algorithm into troubles. A first problem is apparent: the cloning method will typically fail whenever the largest exit rate $Y_{C_{t}[m]}$ among the set of $N_{c}$ copies at a given time becomes of the order of $N_{c}$ itself Hurtado and Garrido (2009b). If this is the case, configuration $C_{t}[m]$ will overpopulate all the other copies after the replication step (see Fig. 1 and steps (b)-(c) in the algorithm of previous section), hence introducing a bias in the Monte Carlo sampling.

The emergence of this problem depends on the magnitude of the current fluctuation, or equivalently on the value of conjugated parameter $\lambda$. The idea now is to estimate the critical $\lambda_{c}$ beyond which this main source of error becomes dominant, using tools from extreme value statistics Sornette (2006). In order to to do so, we first define for a fixed parameter $\lambda$ the maximum exit rate among the set of $N_{c}$ copies at a given time $t$,

[TABLE]

and consider the probability $P_{\lambda,<}^{\text{max}}(y)$ that $Y_{t}^{\text{max}}$ is smaller than a threshold $y$. If $P_{\lambda}^{>}(y)$ is the probability that the exit rate of one clone is larger than $y$, and we assume statistical independence among the different clones, then $P_{\lambda,<}^{\text{max}}(y)$ can be simply written as

[TABLE]

In the limit of large $N_{c}$ the main contribution to $P_{\lambda,<}^{\text{max}}(y)$ is dominated by the tails of $P_{\lambda}^{>}(y)$, where this probability distribution is small, so in this large-$N_{c}$ limit we have

[TABLE]

Consider now the following question: what is the value $y^{*}_{\lambda}(p)$ of the maximum which will not be exceeded with probability $p$? By definition, this value can be obtained from $p=P_{\lambda,<}^{\text{max}}[y^{*}_{\lambda}(p)]$, which according to Eq. (27) leads to the following expression

[TABLE]

Now, for obvious reasons the maximum exit rate allowed by the algorithm is $N_{c}$ itself, and this limit leads to a critical value of the current-conjugated parameter $\lambda$ for a given confidence limit $p$ beyond which sampling problems in the cloning Monte Carlo algorithm will become evident Hurtado and Garrido (2009b). To obtain this critical parameter $\lambda_{c}(p)$ (beyond which the maximum that will not be exceeded with probability $p$ is larger than $N_{c}$) we set $y^{*}_{\lambda_{c}}(p)=N_{c}$ in the previous identity, arriving at

[TABLE]

This condition signals (with confidence level $p$) the onset of the systematic bias due to the finite number of clones in simulations, and defines the critical current-conjugated parameter $\lambda_{c}(p)$ beyond which this happens. To further proceed, we need the probability $P_{\lambda}^{>}(y)$ that the exit rate of a clone is larger than $y$, which implies some knowledge of the statistics at the end of a rare event Hurtado and Garrido (2009b).

We hence consider the probability $P_{t}(C;Q)$ that the system is in configuration $C$ at time $t$ with a total time-integrated current $Q$. This probability obeys a simple recurrence relation (a sort of master equation) of the form

[TABLE]

Iterating in time the previous relation, we arrive at

[TABLE]

where we have used the path probability $\mathbb{W}_{C_{0}\ldots C}=W_{CC_{t-1}}\ldots W_{C_{1}C_{0}}P(C_{0})$ for brevity, see also Eq. (11). Note that $P_{t}(Q)=\sum_{C}P_{t}(C;Q)$ as expected, see Eq. (12) above. In this way, the probability measure of a configuration at the end of a large deviation event of current $\hat{q}=Q/t$, or endtime statistics, can be written as

[TABLE]

Introducing now the moment-generating function of the distribution $P_{t}(C;Q)$,

[TABLE]

such that $Z_{t}(\lambda)=\sum_{C}Z_{t}(C;\lambda)$, see Eq. (15), one can show that for long times $t\gg 1$ Hurtado et al. (2014); Hurtado and Garrido (2010)

[TABLE]

where ${\hat{q}}_{\lambda}=\partial_{\lambda}\theta(\lambda)$ is the conjugated current to parameter $\lambda$, which maximizes the Legendre transform (5). A direct inspection of Eq. (III.2) at the light of the cloning Monte Carlo method of the previous section, see also Eq. (21), shows that $P_{\lambda}^{\text{end}}(C)$ is proportional to the number of clones at a given time in configuration $C$ out of the total number of clones $N_{c}$.

Once $P_{\lambda}^{\text{end}}(C)$ has been defined, we can write the probability that a given clone has an exit rate $Y$ as

[TABLE]

with $Y_{C}=\sum_{C^{\prime}}\tilde{W}_{C^{\prime}C}$. In this way

[TABLE]

where $H(x)$ is the Heaviside step function. Therefore the knowledge of endtime statistics during a large deviation event allows for an estimation of the range of validity of the cloning Monte Carlo method Hurtado and Garrido (2009b). This calculation can be made completely explicit in particular models where analytical expressions for the endtime distribution $P_{\lambda}^{\text{end}}(C)$ can be obtained Hurtado and Garrido (2009b). In these cases it has been found that the probability distribution of exit rates exhibits power law tails, i.e. $P_{\lambda}(Y)\sim Y^{-\alpha(\lambda)}$ in the limit of large $Y$ which dominates the onset of the bias, with an exponent $\alpha(\lambda)$ which is typically a rational function of $\lambda$. This power law behavior is then inherited by the cumulative distribution $P_{\lambda}^{>}(y)\sim y^{-[\alpha(\lambda)-1]}$, and this can be used together with Eq (28) to arrive at an expression for the critical value of the current-conjugated parameter $\lambda_{c}(p)$ beyond which the bias due to the finite number of clones of the simulation appears, i.e.

[TABLE]

The accuracy of this expression, which bounds the range of validity of the cloning method, has been corroborated in numerical experiments of simple stochastic lattice gases Hurtado and Garrido (2009b). Interestingly, the logarithmic dependence on the number of clones $N_{c}$ in the previous equation suggests that an appreciable increase of the range of validity of the algorithm in $\lambda$-space demands an exponential increase in the number of clones.

III.3 Statistics and averages during a rare event

In previous sections we have shown how to sample the dynamical free energy associated with the current (i.e. the Legendre transform of the current LDF) using a technique known as the cloning Monte Carlo method. Our aim in this paper is to apply this method to understand the physics of dynamical phase transitions, which appear as singularities in the current LDF or the dFE and are accompanied by peculiar changes in the trajectories associated with rare events, which may include symmetry-breaking phenomena as captured by certain order parameter. These features can be only determined if we can measure observables of interest (as e.g. an order parameter, or a density profile) during a large deviation, and this implies the characterization of the statistics of configurations in the course of a rare event.

In order to define this statistics, we first introduce the probability $P_{t}(C_{t^{\prime}},t^{\prime};Q)$ that the system of interest was in configuration $C_{t^{\prime}}$ at time $t^{\prime}$ when at time $t$ the total integrated current is $Q$. Most importantly, time scales are such that $1\ll t^{\prime}\ll t$, so all times involved are large enough for the asymptotic large-deviation regime to be reached. To write down a formula for this probability, we just have to realize that this probability is nothing but the likelihood of all paths which traverse configuration $C_{t^{\prime}}$ at time $t^{\prime}$ such that the total accumulated current is $Q$. This reasoning leads to the following expression

[TABLE]

Note that the sum over paths leaves fixed the configuration $C_{t^{\prime}}$ at time $t^{\prime}$. The probability of an arbitrary configuration $C$ during a large deviation event of the current ${\hat{q}}=Q/t$ (also known as midtime statistics) can be now written as

[TABLE]

and does not depend on the duration $t$ of the large deviation event and the time of measurement $t^{\prime}$ as far as $1\ll t^{\prime}\ll t$. As before, if we define now the moment-generating function of the above distribution,

[TABLE]

it can be shown easily that the probability weight of configuration $C_{t^{\prime}}$ at intermediate time $t^{\prime}$ in a large deviation event of current ${\hat{q}}=Q/t$ can be written as Hurtado et al. (2014); Hurtado and Garrido (2010)

[TABLE]

for long times such that $1\ll t^{\prime}\ll t$, and where ${\hat{q}}$ and $\lambda$ are Legendre conjugated duals, ${\hat{q}}={\hat{q}}_{\lambda}=\partial_{\lambda}\theta(\lambda)$. To measure the typical value of an observable during a large deviation event, we have to define averages over the midtime statistics introduced above. In this way, if we denote our observable as $A(C)$, its average value during a large deviation event of the current with conjugated parameter $\lambda$ (i.e. with a current value of $\hat{q}_{\lambda}$) can be written as

[TABLE]

This expression can be now conveniently read from a computational point of view to obtain an algorithm to measure averages of observables of interest during a large deviation event Tailleur and Lecomte (2009). In particular, the idea is to run the cloning algorithm as described in previous sections, so that it generates trajectories typical of a current fluctuation ${\hat{q}}_{\lambda}$. The difference now is that, when the time evolution reaches a chosen intermediate time $t^{\prime}$, such that $1\ll t^{\prime}\ll t$, the value of the observable $A$ associated with the actual configuration is attached to each clone, see Eq. (40). In this way, whenever a clone is replicated, it carries a record of its value of observable $A$ at time $t^{\prime}$. The midtime average of the observable of interest is now obtained from the average of the values of $A$ attached to each surviving clone after a long time interval has passed since the measurement time $t^{\prime}$ Tailleur and Lecomte (2009). Since midtime statistics does not depend on $t$ or $t^{\prime}$ as long as $1\ll t^{\prime}\ll t$, to better sample the midtime average it is convenient to nest in time several measurement-average pairs (characterized by a lag time $t_{\text{lag}}=t-t^{\prime}$), as illustrated in Fig. 2, leaving a reasonable waiting time $t_{\text{wait}}$ to avoid possible correlations between consecutive measurements. Note that since usually $t_{\text{wait}}<t_{\text{lag}}$, several measurements must be carried out and transported in parallel before averaging.

This midtime average, which includes the effect of the population cloning (and hence the exit rated $Y(C)$ along the different trajectories) is clearly different from the naive average of the observable $A(C)$ among the different clones at a given time, which yields the endtime statistics associated with a large deviation event introduced in the previous section. It can be shown Hurtado et al. (2014); Hurtado and Garrido (2010); Jack and Sollich (2010), based on a spectral analysis of the deformed matrix with elements $\tilde{W}_{C^{\prime}C}$, that endtime and midtime statistics during a large deviation event are different. Indeed midtime statistics captures the physics of the system during the rare event, while the endtime statistics only characterize transient behavior which appears both at the end and at the beginning of a rare event Bertini et al. (2015); Derrida ; Jack and Sollich (2010).

III.4 Continuous-time cloning algorithm

The previous time-discrete dynamics that allows for the measurement of time-integrated rare events can be readily extended to continuous-time dynamics. The evolution of the probability $P_{t}(C)$ of being in configuration $C$ at time $t$ is then governed by the master equation

[TABLE]

where $W_{CC^{\prime}}$ is now the transition rate from $C^{\prime}$ to $C$ and $R(C)=\sum_{C^{\prime}}W_{C^{\prime}C}$ is the escape rate from configuration $C$. On the other hand, the probability of being in configuration $C$ having a time-integrated current $Q$ up to time $t$ is $P_{t}(C;Q)$, which evolves according to

[TABLE]

As in the discrete-time case, it is more convenient to work with the Laplace transform of $P_{t}(C;Q)$

[TABLE]

whose evolution is given, by virtue of Eq. (41), by

[TABLE]

Here $\tilde{W}_{CC^{\prime}}=e^{\lambda q_{CC^{\prime}}}W_{CC^{\prime}}$ stands for modified transition rates, which can be much larger than the original $W_{CC^{\prime}}$ for positive currents when $\lambda>0$ and for negative currents when $\lambda<0$. Summing over every configuration $C$ we get the dynamical partition function

[TABLE]

Thus, given $\lambda$ we could get $\theta(\lambda)$ by virtue of (43) just by simulating the evolution of $Z(C;\lambda)$ by means of (42) up to time $t$ for many different realizations. However, this is not straightforward as Eq. (42) is not a stochastic equation since $R(C)\neq\tilde{R}(C)\equiv\sum_{C^{\prime}}\tilde{W}_{C^{\prime}C}$. Nevertheless, it can be rewritten as a stochastic evolution which alternates with a cloning process

[TABLE]

with $\gamma(C)=\tilde{R}(C)-R(C)$. This dynamics can be computationally carried out by considering $N_{c}$ copies of the system which evolve as follows Tailleur and Lecomte (2009) (see Fig. 3).

1. (a)

Set the time to $t_{\alpha}$, with $t_{\alpha}$ being the time of the first clone $\alpha$ to evolve. 2. (b)

Change the configuration of copy $\alpha$ from configuration $C_{\alpha}$ to $C^{\prime}_{\alpha}$ with probability $\tilde{W}_{C^{\prime}_{\alpha}C_{\alpha}}/\tilde{R}(C_{\alpha})$. Compute the sojourn time $\Delta t$ of copy $\alpha$ in configuration $C^{\prime}_{\alpha}$ until the next jump, distributed according to the Poisson law $P(\Delta t)=\tilde{R}(C^{\prime}_{\alpha})e^{-\Delta t\tilde{R}(C^{\prime}_{\alpha})}$. 3. (c)

Clone configuration $C^{\prime}_{\alpha}$ with rate $Y(C^{\prime}_{\alpha})=e^{\Delta t\gamma(C^{\prime}_{\alpha})}$, i.e. generate a number $K_{\alpha}=\lfloor Y(C^{\prime}_{\alpha})\rfloor+1$ of identical clones with probability $Y(C^{\prime}_{\alpha})-\lfloor Y(C^{\prime}_{\alpha})\rfloor$; or $K_{\alpha}=\lfloor Y(C^{\prime}_{\alpha})\rfloor$ otherwise. If $K_{\alpha}=0$, configuration $C^{\prime}_{\alpha}$ is erased.

After a long time $t$, this cloning protocol changes the initial number of copies, resulting in an exponential growth or decrease of the population given by Eq. (43). However, as we have to deal with a finite number of clones, we resize homogeneously the population to the original number of copies after each evolution:

1. (d)

If $K_{\alpha}=0$ then one randomly chosen copy $\beta\neq\alpha$ is replicated, while if $K_{\alpha}>1$ then $K_{\alpha}-1$ copies, among the total $N_{c}+K_{\alpha}-1$, are randomly chosen and erased.

The resizing factor $X_{1}=N_{c}/(N_{c}+K_{\alpha}-1)$ is then stored at each time and the dynamical free energy can be computed as

[TABLE]

where $M$ is the total number of jumps up to time $t$. Fig. 3 sketches this procedure.

It is worth mentioning that one could alternatively analyse the above continuous-time system by simulating each clone for a fixed time $\Delta t$ and then using the discrete-time cloning method with $\Delta t$ as time step, see Brewer et al. (2018) for details. However, the statistical errors strongly depend on $\Delta t$, but these can be considerably mitigated with a more efficient method of producing the new population by reducing the deletion of clones as described in Brewer et al. (2018).

III.5 Recent improvements of the Monte Carlo cloning algorithm

The cloning method described above in its discrete and continuous versions, despite being very successful, suffers from the fact that an exponential number of clones is needed to properly simulate trajectories associated with extremely rare events –corresponding to large values of $\lambda$–, as shown by (35), or in the vicinity of dynamical phase transitions. This has motivated recent improvements of the method aiming to mitigate such finite-population effects. One possibility recently explored consists in introducing an interpolation technique for the large deviation function based on the analysis of the systematic errors of a birth-death process Nemoto et al. (2017a). An extension of such technique has led to a more efficient version of the cloning algorithm in continuous time which significantly improves the large deviation function estimators in the long time and large number of clones limit Guevara Hidalgo et al. (2017). In addition, some rigorous bounds on convergence properties of the algorithm have been established Angeli et al. (2018, 2019).

Another upgrade of the cloning method has been developed by incorporating to the population dynamics a controlling force which is determined by iterating a measurement-and-feedback scheme Nemoto et al. (2016). This method, recently extended to Markov jump processes Nemoto et al. (2017b), is independent of the choice of the force in the limit of large number of clones. However, for a finite population of clones, it drastically improves the accuracy of the results, being necessary to simulate just one clone for the optimal force, as in that case the cloning factor $\gamma(C)$ in Eq.(44) is constant, i.e. independent of the configuration, and equal to $\theta(\lambda)$. Thus, the dynamics associated with a rare fluctuation is stochastic and the rare trajectories of the original system correspond to the typical trajectories of the system under the action of the optimal force. The latter can be exactly obtained in terms of the eigenvector of the tilted generator associated with its largest eigenvalue $\theta(\lambda)$, but this is a very challenging task for many-particle systems, both numerically and analytically. Remarkably, the effective dynamics under the optimal force corresponds to the Doob’s h-transform Jack and Sollich (2010); Chetrite and Touchette (2015a, b), which transforms the non-physical tilted dynamics into a proper stochastic dynamics. Nevertheless, despite the exact optimal force is difficult to obtain, combining an iterative estimation with the population dynamics has shown to be highly effective in finite-size scaling analysis of the first order dynamical phase transition present in the Fredrickson-Andersen model Nemoto et al. (2017b).

A similar approach has been devised in Ref. Ray et al. (2018), where the control strategy is performed by approximating the eigenvectors by means of guided distribution functions. In particular, the eigenvectors are approximated by product states of clusters of different number of sites. In this case, the iterative procedure to determine the eigenvectors is carried out outside the dynamics itself, in contrast to the approach of Refs. Nemoto et al. (2016, 2017b). Analogous methods have been proposed in the context of the umbrella sampling technique Oakes et al. (2018), in which a reference dynamics is introduced to perform importance sampling of the rare trajectories of interest Klymko et al. (2018). In addition, besides the improvement in the copying and selection method pointed out at the end of the last section, further progress has been made regarding the parallel implementation of the cloning algorithm in Brewer et al. (2018), as well as its generalization to non-Markovian systems Cavallaro and Harris (2016). Finally, it is worth mentioning that another recent algorithm using risk-sensitive and feedback control methods has been introduced to estimate the effective dynamics of a rare event from a single trajectory Ferré and Touchette (2018).

IV Time-translational symmetry breaking at the trajectory level

Once equipped with the general computational tools described in previous sections to sample rare current fluctuations, we set out to study its application to explore different dynamical phase transitions which appear at the trajectory level. In particular, in this section we will use the cloning Monte Carlo method to characterize the DPT into a dynamical phase with broken time-translation invariance, appearing in two paradigmatic models of transport for nonequilibrium physics, namely the Kipnis-Marchioro-Presutti (KMP) model of heat conduction Kipnis et al. (1982) and the weakly asymmetric simple exclusion process (WASEP) Derrida (1998).

IV.1 Ballistic energy packets in a model of heat conduction

As we have already pointed out in the introduction, DPTs arise in the large size and/or long time limit and manifest as a non-analyticity in the LDF of time-averaged observables such as the current. Just as in standard critical phenomena, DPTs can be either discontinuous (first-order) or continuous (second-order). The latter are associated with a symmetry breaking phenomenon, i.e. the optimal path leading to a rare fluctuation breaks a symmetry of the corresponding large deviation functional (8). In this section we shall focus on an equilibrium model of energy diffusion featuring a continuous DPT which breaks the time-translational invariance.

In 1982 Kipnis et al. (1982), C. Kipnis, C. Marchioro and E.Presutti introduced a simple stochastic model to understand in rigorous terms energy transport in systems with many degrees of freedom. This model, denoted hereafter as KMP model, has been paramount for the study of nonequilibrium phenomena since its formulation, leading to a number of new ideas and breakthroughs in this field. In particular, the authors of Kipnis et al. (1982) were able to prove rigorously, starting from its microscopic Markovian dynamics, that the $1d$ KMP model in contact with two boundary thermal baths at different temperatures obeys Fourier’s law. This is a key result in mathematical physics, as the microscopic foundations of Fourier’s law still remain unknown in more realistic systems Bonetto et al. (2000).

The KMP model consists of a $1d$ lattice of $L$ sites, see Fig. 4(a)-(b), though its definition can be generalized to different lattices in arbitrary dimension. The system’s microscopic configuration is defined by a set $C\equiv\{e_{i},i=1,\ldots,L\}$, where $e_{i}\in\mathbb{R}^{+}$ is the energy of site $i$. The dynamics is stochastic, and time can be discrete or continuous. In an elementary step, a pair of nearest neighbor sites $(i,i+1)$ is chosen at random, and the pair total energy is randomly redistributed locally so that the pair total energy is conserved in the interaction. In this way $(e_{i},e_{i+1})\to[e^{\prime}_{i}(p),e^{\prime}_{i+1}(p)]$ with

[TABLE]

where $p$ is a uniform random number $p\in[0,1]$, so that $(e_{i}+e_{i+1})=[e^{\prime}_{i}(p)+e^{\prime}_{i+1}(p)]$ $\forall p$. With the aim of studying energy transport under a temperature gradient, KMP Kipnis et al. (1982) coupled the model so defined to two boundary thermal baths at different temperatures, see Fig. 4(a). However, for our purposes here it is convenient to study a microcanonical version of the KMP model, i.e. defined on a closed $1d$ lattice under periodic boundary conditions, see Fig. 4(b). In this case, the system is isolated and the total energy per particle, $\rho_{0}\equiv L^{-1}\sum_{i=1}^{L}e_{i}$, is conserved in the evolution. We note that in both cases the model can be generalized to include an external field $E$ driving the energy field is some preferred direction, though we will mostly focus on the $E=0$ KMP case below.

The macroscopic limit of this model is taken under diffusive scaling, with $x=i/L$ the macroscopic spatial variable and $t=m/L^{2}$ the macroscopic time ($m$ is the microscopic time variable), so the macroscopic energy density is defined as the average energy on each site; $\rho(x,t)=\langle e_{i=Lx}\rangle_{m=tL^{2}}$. In Kipnis et al. (1982) it was shown by means of the duality technique Carinci et al. (2013) that, at the macroscopic level, this model obeys Fourier’s law with a diffusivity $D(\rho)=1/2$. As a consequence, when in contact with reservoirs at different temperatures $T_{\text{L}}\neq T_{\text{R}}$ at the left (L) and right (R) boundaries, the system reaches in the long time limit a stationary state characterised by a linear energy profile and a constant energy current, i.e.

[TABLE]

On the other hand, when boundary conditions are periodic, i.e. in the absence of contact with thermal reservoirs, the KMP steady state is homogeneous, with a constant energy profile and no net current

[TABLE]

Moreover, the KMP mobility transport coefficient can be shown to be $\sigma(\rho)=\rho^{2}$ Prados et al. (2012); Hurtado et al. (2014), and from these two transport coefficients, $D(\rho)$ and $\sigma(\rho)$, the macroscopic fluctuation theory (MFT) Bertini et al. (2015) of Section §II offers precise predictions on the fluctuating behavior of the KMP model, both under a temperature gradient Hurtado and Garrido (2009b, 2010) and under periodic boundary conditions Hurtado and Garrido (2011).

For the case considered here, i.e. current statistics in the KMP model under periodic boundary conditions, MFT predicts a dynamical phase transition from (a) a phase for small current fluctuations around the average $q_{st}=0$ where the Additivity Principle holds so the optimal profiles solution of the MFT problem are time-independent, $\rho_{q}(x,t)=\rho_{0}$ and $j_{q}(x,t)=q$ (see Section §II), to (b) another dynamical phase for large current fluctuations where time-dependent optimal profiles emerge. Indeed, one can argue that small current fluctuations around the average result from the sum of local and uncorrelated jumps of the energy field, hence leading to Gaussian current statistics in agreement with the central limit theorem. However, for large enough currents beyond a critical threshold $|q_{c}|$, with

[TABLE]

the system finds optimal to pack energy into a coherent traveling wave $\rho_{q}(x,t)=\omega_{q}(x-vt)$, moving at constant velocity $v$, which facilitates the rare event. Fig. 4(c) depicts the dynamical phase diagram for this second-order DPT, and includes the effect of a possible external field $E$ on the critical current Hurtado et al. (2014). For the particular case $E=0$ studied here, the appearance of ballistic energy packets for large current fluctuations is remarkable as it happens in an isolated (microcanonical) equilibrium system, spontaneously breaking time-translational symmetry in $1d$. This is a beautiful example of the general observation that symmetry-breaking phase transitions forbidden in $1d$ equilibrium steady states can, however, emerge in the statistics of trajectories.

Our purpose in this section consists in characterizing the DPT in the KMP model using the cloning Monte Carlo method described in previous sections, providing some tips and tricks of the trade in the computation of large deviation functions, as well as details on order parameters for the transition and a characterization of the emergent structures. For simplicity we focus below on the discrete-time version of the cloning method as applied to the KMP model. If $C\equiv\{e_{i},i=1,\ldots,L\}$ is the system configuration at a given time and

[TABLE]

is the configuration resulting from $C$ after an interaction has occurred at pair $(k,k+1)$ with an exchange parameter $p$, see Eqs. (45)-(46), then the transition probability between these two configurations is simply $W_{C_{k}^{(p)},C}=L^{-1}$, as the interaction probability of all $L$ pairs is the same and $p$ is an homogeneous random number in the unit interval. We now define the energy current involved in this transition as the energy traveling to the right divided by the total number of pairs, namely

[TABLE]

so the modified transition probability $\tilde{W}_{C_{k}^{(p)},C}$ defined in the cloning Monte Carlo method, see Eq. (14), can be written as

[TABLE]

The associated exit rate $Y_{C}$, see Eq. (16), hence reads

[TABLE]

so the normalized modified dynamics to study current statistics in the $1d$ periodic KMP model is just

[TABLE]

see Section §III.1. Note that this modified dynamics weights each possible transition according to the current involved in the step and the intensive parameter $\lambda$ conjugated to the current. This new stochastic dynamics can now be used in conjunction with the clone population replication/pruning scheme of §III.1 to sample the current large-deviations statistics in this model.

A first evidence of the DPT phenomenon described above for the KMP model is expected to appear at the level of individual trajectories of the energy field for current fluctuations below and beyond the critical current $|q_{c}|$. With this aim we simulated the KMP model on a $1d$ ring lattice using the cloning method for two different values of the conjugated parameter $\lambda$, such that in one case the current $q=q_{\lambda}$ associated with this $\lambda$ was below the critical threshold, $q<|q_{c}|$, while in the other case $q>|q_{c}|$. Figs. 5(a)-(b) show the resulting typical trajectories in both cases as a spatiotemporal raster plot of the energy field. In particular, Fig. 5(a) shows that the system can sustain a net current fluctuation $q<|q_{c}|$ in the absence of any macroscopic structure, i.e. with an homogeneous average density profile. However, Fig. 5(b) clearly demonstrates that in order to sustain a large energy current ($|q|>|q_{c}|$) the system accumulates energy in the form of packets that propagate ballistically (with a constant velocity on average) across the lattice. By contrast, trajectories producing moderate currents do not feature any spatial structure.

A more quantitative characterization of these energy packets for $|q|>|q_{c}|$ can be obtained by measuring their average shape. Due to the packet’s ballistic motion and the lattice periodicity, a naive averaging of density profiles along multiple rare-event trajectories would miss the packet structure, leading instead to homogeneous density profiles. In order to correctly measure the traveling-wave energy profile and not to blur away its spatial structure, we performed averages of the density field around its instantaneous center of mass. Indeed, due to the periodicity of the $1d$ lattice we can consider the system as a ring embedded in a two-dimensional space, and this allows us to assign an angle $\theta_{i}=2\pi i/L$ to each site $i\in[1,L]$ in the lattice. In this way we define now the angular position of the energy field center of mass for a given configuration $C\equiv\{e_{i},i=1,\ldots,L\}$ as $\theta_{\text{CM}}=\tan^{-1}(y_{\text{CM}}/x_{\text{CM}})$, where

[TABLE]

Hence, every time we have to average configurations during a large deviation current event, we rotate the system by an angle $\theta_{\text{CM}}$ before averaging so that the angular center of mass always sits in the origin. In this way we correctly capture the spatial structure of the traveling-wave profile, which would otherwise fade away. Some caution is needed however: this measuring method leads to spurious weakly-structured energy profiles for subcritical current fluctuations $|q|<|q_{c}|$, as averaging random homogeneous energy profiles around their random center of mass results in a weak but non-trivial spatial structure. This potential problem can be easily fixed by noting that this spurious profile is of course independent of the current $q$ for $|q|<|q_{c}|$, so it can be easily subtracted. Fig. 5(c) shows the measured average optimal energy profile responsible for a current fluctuation $q_{\lambda}$ obtained with this technique as a function of the conjugate parameter $|\lambda|$. The appearance of a critical value $\lambda=\lambda_{c}$, or equivalently a critical current $q_{c}$, where the optimal profile changes from homogeneous to structured clearly signals the onset of the DPT into the traveling-wave dynamical phase. These measurements were performed for a system with $L=32$ sites using a large number of clones $N_{c}=10^{4}$.

As described in Section §III.2, these measurements are strongly affected by finite size effects both on the system sizes and the number of clones employed in the sampling, so a careful analysis of how $L$ and $N_{c}$ affect the measurements is always needed. Indeed, in the inset to Fig. 6 we display the measured energy packet profile for a current fluctuation in the traveling-wave phase ($\lambda=4>\lambda_{c}$) for increasing system sizes $L$. We observe how the measured profiles suffer from the strong finite-size effects anticipated above, but for $L=32$ there is already a very good agreement with the optimal energy profiles predicted by MFT Hurtado and Garrido (2011); Hurtado et al. (2014). A similar analysis can be performed as a function of the number of clones $N_{c}$, and this suggests that $N_{c}=10^{4}$ suffices in this case to sample accurately current statistics and the associated structures across the DPT. As we will see below, these finite-population effects increase in higher dimensions, demanding the simulation of very large numbers of clones in parallel.

Another key feature of the supercritical energy packets is their propagation velocity. To measure this observable, we define a quasi-instantaneous velocity for a given trajectory as the average slope of the packet peak trajectory in a space-time diagram for a short time interval, see Fig. 5(b). In particular, by fitting the motion of the center of mass during small time intervals $\Delta t$ to a ballistic law, $r(t+\Delta t)-r(t)=vt$, see Fig. 5(b), and using this observable to obtain midtime averages using the sampling scheme of Section §III.3, we obtained the data plotted in the main panel of Fig. 6. In particular, we plot in this graph the velocity of the center of mass motion as a function of the bias parameter $\lambda$ for different system sizes, as well as the macroscopic theoretical prediction based on MFT calculations Hurtado and Garrido (2011). An excellent agreement is found already for $L=32$, while clear finite-size corrections are observed below this size, as in other observables above. Note that while the center of mass velocity is linear in $\lambda$ in the subcritical region (i.e. this velocity is just proportional to the current), for currents beyond the critical threshold, or equivalently $|\lambda|>|\lambda_{c}|$, the velocity $v(\lambda)$ becomes a nonlinear monotonously increasing function of $\lambda$, another trait of the DPT.

Finally, the cloning algorithm permits to measure not only the typical profiles giving rise to a large dynamical fluctuation but also its dynamical free energy $\mu(\lambda)$, which is nothing but the Legendre transform of the current LDF $G(q)$, with $\lambda$ being the parameter conjugate to the current $q$. In Fig.7(a) we present the theoretical $\mu(\lambda)$ along with the numerical values measured for different system sizes and $N_{c}=10^{4}$. We clearly observe how for currents below the critical one ($|\lambda|<|\lambda_{c}|$) the Gaussian solution –associated with the time-independent homogeneous profiles– is the minimizing one, while for currents above the critical threshold ($|\lambda|>|\lambda_{c}|$) the dFE branch associated with traveling-wave profiles is the optimal one [see inset to Fig. 7(a)]. Note that in virtue of the Legendre transform connecting $G(q)$ and $\mu(\lambda)$ given by (10), we get that $P_{\rm MFT}(q)>P_{\rm flat}(q)$ for $|q|>|q_{c}|$ implies that $\mu_{\rm MFT}(\lambda)>\mu_{\rm flat}(\lambda)$ for $|\lambda|>|\lambda_{c}|$, as observed in Fig. 7(a). This highlights the fact that large current fluctuations are far more probable –via the emergence of complex spatio-temporal structures– than what is predicted naively by Gaussian statistics. The transition between both regimes corresponds to a second-order DPT as $\mu(\lambda)$ presents a discontinuity in its second derivative, i.e. $\lim_{|\lambda|\to|\lambda_{c}|^{-}}\partial^{2}_{\lambda}\mu(\lambda)\neq\lim_{|\lambda|\to|\lambda_{c}|^{+}}\partial^{2}_{\lambda}\mu(\lambda)$. This is displayed in Fig. 7(b) as obtained from MFT Hurtado and Garrido (2011), where we also present its first derivative $\mu^{\prime}(\lambda)$ corresponding to the current $q_{\lambda}$ associated with each $\lambda$.

IV.2 Jammed trajectories in a model of particle diffusion

Another archetypical stochastic lattice model which displays an interesting DPT is the weakly asymmetric simple exclusion process (WASEP) Gärtner (1987); De Masi et al. (1989), a schematic model of particle diffusion. This model exhibits a DPT in its current statistics which is mathematically very similar to that of the KMP model described above Bodineau and Derrida (2005); Pérez-Espigares et al. (2013). However, as we will see below, the DPT in WASEP has a markedly different physical interpretation, which again gives valuable insights into the range of possibilities that an apparently simple system has at its disposal in order to maximize the probability of a given fluctuation. Moreover, understanding this DPT in $1d$ will help us to better understand the important effect of dimension on the current statistics of driven system, which will be the focus of the Section §V below.

The WASEP is defined on a $1d$-lattice of $L$ sites, each of which can be either empty or occupied at most by one particle, see Figs. 8(a)-(b). As for the KMP model, the WASEP can be defined in arbitrary dimensions and different types of lattices (indeed we will study a $2d$ case in next section). A microscopic configuration is defined as $C\equiv\{n_{i},i=1,\ldots L\}$, where $n_{i}=0,1$ is the occupation number of site $i$ at a given time, which as before can be either discrete or continuous. Every particle can randomly jump to one of its nearest neighboring sites, provided this is empty. The jump rate is $p^{\pm}=pe^{\pm E/L}$, for jumps to the right ($+$) or left ($-$) neighbor, with $p$ an overall jump rate which sets the microscopic time scale, and $E$ an external field biasing the motion along a given direction. The term weakly in WASEP denotes the fact that the rate asymmetry, $p^{+}-p^{-}$, scales as $L^{-1}$ for large $L$, allowing the external field to compete on equal footing with the particle diffusion process. As in the KMP model of previous section, we may choose to (a) couple the WASEP model to two particle reservoirs at the boundaries, which may work at different chemical potentials (or densities) thus inducing an external gradient, or rather (b) impose periodic boundary conditions so the total number of particles in the system $M\equiv\rho_{0}L=\sum_{i=1}^{L}n_{i}$ is conserved during the evolution. To implement boundary reservoirs, particles are injected to the first (last) site at rate $\alpha$ ($\delta$) provided it is empty and removed from the first (last) site at rate $\gamma$ ($\beta$) provided it is occupied [see Fig. 8(a)]. This mechanism mimics the effect of boundary reservoirs at densities $\rho_{\text{L}}=\alpha/(\alpha+\gamma)$ and $\rho_{\text{R}}=\delta/(\delta+\beta)$. On the other hand, removing the extremal reservoirs and introducing periodic boundary conditions we obtain the microcanonical version of the WASEP that will be the focus of study in this section. The stationary state under periodic boundaries has a homogeneous density profile $\rho_{st}(x)=\rho_{0}=M/L$ equal to the average density in the system, and a net current $q_{st}(E)=\sigma(\rho_{0})E=2p\rho_{0}(1-\rho_{0})E$ due to the presence of the weak external field.

Averaging over trajectories on each site and applying the diffusive scaling limit, we can define a macroscopic density field for WASEP as $\rho(x,t)=\langle n_{i=Lx}\rangle_{m=tL^{2}}$, with $m$ the microscopic time. In this mesoscopic limit, the system evolution is captured by a fluctuating hydrodynamics equation, see Eq. (6) in Section §II, defined in terms of two transport coefficients, the diffusivity $D(\rho)$ and the mobility $\sigma(\rho)$, which for WASEP are $D(\rho)=p$ and $\sigma(\rho)=2p\rho(1-\rho)$ Spohn (2012). Using these two transport coefficients, the macroscopic fluctuation theory of Section §II predicts a DPT in the current statistics of the microcanonical $1d$ WASEP to a dynamical phase dominated by traveling-wave density profiles, with a mathematical structure very similar to that of the DPT in the KMP model but with several key differences. For instance, the critical current where the DPT appears has the same general form as in the KMP model, see Eq. (50), namely

[TABLE]

However, and most importantly, while $\sigma^{\prime\prime}(\rho_{0})>0$ for the KMP model, for the WASEP we have $\sigma(\rho)=2p\rho(1-\rho)$ so that $\sigma^{\prime\prime}(\rho)=-4p<0$. This key difference, which originates in the particle exclusion –a main feature of the WASEP–, immediately introduces a critical value for the external field,

[TABLE]

below which no DPT is possible as the critical current $|q_{c}|$ becomes complex. Fig. 8(c) shows the resulting dynamical phase diagram for current statistics in $1d$ WASEP as derived from MFT, see Eq. (58). In this way we can write

[TABLE]

so clearly $|q_{c}|<|q_{st}|$ for WASEP under external fields of magnitude larger than the critical one. This is another important difference between the KMP and WASEP DPTs. In particular, for the KMP model there is no current in the steady state (in the absence of external field) and the DPT appears for large enough currents, where the system facilitates the fluctuation by localizing energy in packets that travel ballistically across the system. On the other hand, and in stark contrast, the DPT in WASEP appears for current below the average stationary current. In this case the dynamic phase transition corresponds to the emergence of a macroscopic jammed state which hinders transport of particles to facilitate a current fluctuation well below the average.

To characterize this sort of jamming transition at the trajectory level, we performed Monte Carlo cloning simulations to unveil the current statistics of the $1d$ periodic WASEP Pérez-Espigares et al. (2013) for an average density $\rho_{0}=0.3$ under a supercritical external field $E=10>E_{c}$, and for increasing system sizes $L\in[8,64]$. Fig. 9 summarizes our results, that we obtained using $N_{c}=2\times 10^{4}$ clones evolving in parallel. In particular, Figs. 9(a)-(b) depict spatiotemporal raster plots of trajectories obtained in cloning simulations for currents larger and smaller than the critical current $|q_{c}|$. The trajectory associated with a current fluctuation $|q|>|q_{c}|$ is homogeneous and apparently random, with no clear structure, while for $|q|<|q_{c}|$ particles are localized in a jammed region which moves at constant velocity across the system. The average shape of the jammed density profile was measured in detail using midtime sampling methods, see Figs. 9(c)-(d) and Section §III.3. Although finite-size effects are strong as expected, a neat convergence towards the macroscopic MFT predictions is observed already for sizes $L=64$ when using a sufficiently large number of clones $N_{c}=2\times 10^{4}$. Similar finite-size effects appear in our measurements of the dynamical free energy $\mu(\lambda)$ and its derivatives, see Figs. 9(e)-(f), though again an excellent convergence to the MFT result is observed. Note in particular the non-quadratic shape of the dFE $\mu(\lambda)$ in the traveling-wave region, which corresponds to non-Gaussian current fluctuations in this regime Pérez-Espigares et al. (2013), as well as the second-order character of the transition, see $\mu^{\prime\prime}(\lambda)$ in Fig. 9(f).

V Dynamical criticality in high-dimensional driven systems

As we have already discussed in the introduction, the discovery of DPTs in the fluctuations of many-particle systems is an important finding for nonequilibrium physics. The reason is that the large deviation functions controlling the statistics of these fluctuations are the best candidates we have nowadays to generalize the concept of thermodynamic potentials to systems out of equilibrium, where no bottom-up approach exists yet connecting microscopic dynamics with macroscopic properties. In this way, understanding the singularities appearing in LDFs seems relevant to better understand nonequilibrium phenomena. In addition, the emergence of order associated with rare fluctuations implies that these extreme events are far more probable than previously anticipated Lam et al. (2009); Hurtado et al. (2014); Zarfaty and Meerson . However, up to now most works on DPTs have focused on oversimplified $1d$ models as the ones discussed in previous section Hurtado and Garrido (2011); Pérez-Espigares et al. (2013); Hurtado et al. (2014); Vaikuntanathan et al. (2014); Jack et al. (2015); Shpielberg and Akkermans (2016); Zarfaty and Meerson ; Baek et al. (2017); Karevski and Schütz (2017); Garrahan and Lesanovsky (2010); Ates et al. (2012); Lesanovsky et al. (2013); Carollo et al. (2017) or fluctuations of scalar ($1d$) observables in $d>1$ Garrahan et al. (2007, 2009); Hedges et al. (2009); Chandler and Garrahan (2010); Pitard et al. (2011); Speck et al. (2012); Pinchaipat et al. (2017); Abou et al. (2018); Garrahan et al. (2011); Genway et al. (2012); Manzano and Hurtado (2014). The challenge thus remains to understand DPTs in more realistic settings relevant for actual experiments, as e.g DPTs in the fluctuations of fully vectorial observables in $d$ dimensions and how they are affected by the (possible) system anisotropy.

With this idea in mind, we employed the cloning Monte Carlo technique to study current fluctuations in a $2d$ version of the WASEP model presented in the previous section. In particular we define now the model on a $2d$ square lattice of size $N=L\times L$ with periodic boundaries (i.e. with the topology of a torus) where $M\leq N$ particles evolve, so the global density is $\rho_{0}=M/N$. As before, each site may contain at most one particle. These particles can jump stochastically to neighboring empty sites, which now lie along the $\pm{\alpha}$-direction ($\alpha=x,y$), at a rate $p^{\alpha}_{\pm}\equiv\text{exp}[\pm E_{\alpha}/L]/2$, with ${\bf{E}}=(E_{x},E_{y})$ being in this case a vectorial external field. Fig. 10 shows an sketch of the $2d$ WASEP.

Interestingly, for large ${\bf{E}}$ (needed to overcome a critical field value) and the moderate system sizes that we are capable to explore with the cloning algorithm, the field per unit length ${\bf{E}}/L$ is strong enough to induce an effective anisotropy in the medium. This effective anisotropy (a main feature of many real systems) enhances diffusivity and mobility along the strongest field direction, an effect that can be taken into account within macroscopic fluctuation theory via an effective anisotropy parameter $\epsilon$. Indeed, by expanding the microscopic transition rate $p_{\pm}^{\alpha}$ to second order in the field per unit length, i.e.

[TABLE]

it is easy to show using a simple random walk argument that, while the first-order term gives rise to the standard field biasing dynamics along a given direction, the second-order perturbation results in effective differences of diffusivity and mobility along the different field directions. In this way, the transport coefficients characterizing hydrodynamic behavior in the $2d$ WASEP are a diffusivity matrix ${\hat{D}}(\rho)\equiv D(\rho){\hat{\cal A}}$ and a mobility matrix ${\hat{\sigma}}(\rho)=\sigma(\rho){\hat{\cal A}}$, with $D(\rho)=1/2$ and $\sigma(\rho)=\rho(1-\rho)$ as in the $1d$ WASEP of previous section, and ${\hat{\cal A}}$ being a constant diagonal matrix

[TABLE]

which modelizes the system underlying anisotropy in terms of an anisotropy parameter $\epsilon$.

The variational problem of MFT can be solved in this complex scenario, and predicts a rich dynamical phase diagram for the statistics of the space&time-averaged current vector ${\bf{q}}$ for the $2d$ WASEP, with some surprising new physics when compared to the $1d$ case. Fig. 11 shows the predicted dynamical free energy $\mu({\bf{z}})$, with ${\bf{z}}\equiv{\bm{\lambda}}+{\bf{E}}$, for different anisotropies $\epsilon$ as well as the dynamical phase diagram in each case and the typical trajectories in each dynamical phase. In particular, MFT predicts for large enough external fields $|{\bf{E}}|>E_{c}$ a second-order DPT for currents ${\bf{q}}\cdot{\hat{\cal A}}^{-1}{\bf{q}}=\sigma_{0}^{2}\Xi_{c}$ (or equivalently ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}=\Xi_{c}$), with $\sigma_{0}\equiv\sigma(\rho_{0})$ and $\Xi_{c}$ a critical threshold. This critical DPT line separates a homogeneous fluctuation phase with structureless trajectories and Gaussian current statistics characterized by a quadratic dFE $\mu_{\text{G}}({\bf{z}})\equiv({\bf{z}}\cdot\sigma_{0}{\hat{\cal A}}{\bf{z}}-{\bf{E}}\cdot\sigma_{0}{\hat{\cal A}}{\bf{E}})/2$, and a non-Gaussian dynamical phase for small currents, ${\bf{q}}\cdot{\hat{\cal A}}^{-1}{\bf{q}}\leq\sigma_{0}^{2}\Xi_{c}$ or ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}\leq\Xi_{c}$. As in the $1d$ WASEP, coherent jammed states emerge in the non-Gaussian phase in the form of traveling-wave trajectories, thus breaking time-translational symmetry. Such jammed states hamper particle flow enhancing the probability of low-current fluctuations, but rather counter-intuitively these jammed states in $2d$ are surprisingly extended and noncompact. Interestingly, for mild or no anisotropy, $\epsilon<\epsilon_{c}$, different symmetry-broken $1d$ density waves dominate different current vector regimes, see Figs. 11(a)-(b). Lines of first-order DPTs separate both density-wave phases, and dynamical coexistence emerges along these first-order lines between the two traveling-wave phases, see Figs. 11(g)-(i). However, the competition between the two different density-wave phases for low current fluctuations, modulated by the orientation of the current vector, disappears at a critical anisotropy $\epsilon_{c}\approx 0.035$, beyond which a single dynamical density-wave phase dominates the non-Gaussian regime, see Fig. 11(c).

Our aim here is to search for compelling evidences of this complex DPT in numerical simulations of the $2d$ WASEP, and to characterize in detail the emergent order predicted for low currents. To do so, we explored the statistics of the particle vector current in this model using massive cloning Monte Carlo simulations Giardinà et al. (2006); Lecomte and Tailleur ; Giardinà et al. (2011); Tizón-Escamilla et al. (2017a), see Section §III. In particular, we focused on density $\rho_{0}=0.3$ under a strong external field ${\bf{E}}=(10,0)$ along the $\hat{x}$-direction, such that $|{\bf{E}}|>E_{c}$, investigating in detail finite-size effects up to sizes $N=144$. Moreover, as the number of clones needed to observe a given rare event typically grows exponentially with the system size, all the more the rarer the event is Hurtado and Garrido (2009b); Nemoto et al. (2016), we needed to reach the extraordinary number of $N_{c}=5.12\times 10^{5}$ clones evolving in parallel for a long time to correctly sample the tails of the current distribution.

The main panel in Fig. 12 shows the measured $\mu({\bf{z}})$ for $N=144$ and $N_{c}=5.12\times 10^{5}$ as a function of $z=|{\bf{z}}|$ for different current orientations $\phi=\tan^{-1}(z_{y}/z_{x})$, while the inset analyzes both finite-$N$ and finite-$N_{c}$ corrections in our measurements as a function of $z$ for $\phi=0$. In particular, the inset shows the effect of a varying number of clones on the measured $\mu({\bf{z}})$ for the largest $N=144$: while finite-$N_{c}$ corrections are weak for $z\approx E$ (i.e. ${\bm{\lambda}}\approx 0$ or ${\bf{q}}\approx{\bf{q}}_{st}=\sigma_{0}{\hat{\cal A}}{\bf{E}}$, small fluctuation regime), as otherwise expected, these corrections mount up as $z\to 0$ (equivalently ${\bm{\lambda}}\to-{\bf{E}}$ or ${\bf{q}}\to 0$) tending to overestimate the value of $|\mu({\bf{z}})|$. In any case, these corrections scale as $1/N_{c}$, in agreement with Guevara Hidalgo et al. (2017); Nemoto et al. (2017a), so as to observe an excellent convergence toward the macroscopic fluctuation theory prediction (solid line) for large enough $N_{c}$ and the largest $N$ explored, as displayed in the inset to Fig. 12. On the other hand, finite-$N$ corrections measured for the largest $N_{c}=5.12\times 10^{5}$ affect equally our measurements of the dFE irrespective of the value of $z$, reflecting the distance of the microscopic simulation to the predicted behavior based on a macroscopic approach. These finite-$N$ corrections, which underestimate $|\mu({\bf{z}})|$ for large fluctuations ($z\to 0$), decay at good pace as $N$ increases, showing an excellent convergence to the MFT prediction already for $N=144$ and $N_{c}=5.12\times 10^{5}$.

In this way, and in agreement with MFT, the measured dFE is fully compatible with the Gaussian prediction $\mu_{\text{G}}({\bf{z}})$ for ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}\geq\Xi_{c}$. This confirms the idea that small current fluctuations come out from the random superposition of typically uncorrelated, localized jump events which sum up incoherently to yield Gaussian statistics and structureless typical trajectories (see below). Interestingly, our data exhibit a weak dependence of the dFE $\mu({\bf{z}})$ on the angle $\phi$ in this Gaussian regime, a clear hallmark of the effective anisotropy $\epsilon$ mentioned above. Indeed, by fitting the observed data in the region ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}\geq\Xi_{c}$ to the $\epsilon$-dependent quadratic form $\mu_{\text{G}}({\bf{z}})\equiv({\bf{z}}\cdot\sigma_{0}{\hat{\cal A}}{\bf{z}}-{\bf{E}}\cdot\sigma_{0}{\hat{\cal A}}{\bf{E}})/2$ for each angle $\phi$, see inset in Fig. 13, we can estimate the effective value of the anisotropy parameter for $N=144$, obtaining that $\epsilon\approx 0.038$ properly describes the observed weak anisotropy. This effective anisotropy is slightly larger than the critical anisotropy $\epsilon_{c}\approx 0.035$ beyond which a single symmetry-broken phase dominates the non-Gaussian regime, see Fig. 11(c), an observation consistent with additional results at the trajectory level (see below). Furthermore, Fig. 12 clearly shows that the Gaussian fluctuation regime ends up for ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}<\Xi_{c}$, where systematic deviations from the quadratic form $\mu_{\text{G}}({\bf{z}})$ become apparent. This change of behavior, in excellent agreement with MFT predictions, signals the onset of the DPT into a non-Gaussian current fluctuation regime characterized by traveling density-wave trajectories.

The observed crossover in the dFE at ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}=\Xi_{c}$, see Fig. 12, is a strong evidence of the existence of the predicted DPT. In order to confirm this evidence, and to fully characterize the DPT, we now search for an order parameter exhibiting a smooth but apparent change as we cross the critical line. This is the smoking gun of any continuous phase transition such as the DPT here studied Binney et al. (1992). Moreover, the order parameter introduced must distinguish between the competing symmetry-broken phases (i.e. the different jammed density-wave trajectories) which are expected to appear for low current fluctuations. Such order parameter can be defined by performing a tomographic analysis of configurations, i.e. by taking $1d$ slices of our $2d$ system and analyzing density profiles in each slice. In particular we consider a microscopic particle configuration ${\bf n}=\{n_{ij};i,j\in[1,L]\}$ and slice it along one of the principal axes, say $\hat{x}$, defining the configuration of the $j$-slice as

[TABLE]

see e.g Figs. 14(a) and 14(d). To properly take into account the periodic boundaries [i.e. the system torus topology, see Fig. 14(b) and Fig.14(e)], and in analogy with the angular analysis introduced in Section §IV.1 for the $1d$ case, we consider each $j$-slice as a ring of fixed radius embedded in $2d$, so we can assign to each site $i\in[1,L]$ an angle $\theta_{i}=2\pi i/L$, and compute the angular position of the center of mass for the $j$-slice, $\theta_{\text{cm}}^{(j)}$. This is defined as

[TABLE]

where we have introduced the definitions

[TABLE]

with $M_{j}=\sum_{i=1}^{L}n_{ij}$ the total number of particles in this slice. In this way, a small dispersion of the angular centers of mass across the different slices clearly signals the formation of a coherent jam along the $x$-direction and the associated density wave in the orthogonal direction, see Fig. 14(c). On the other hand, a large dispersion of $\theta_{\text{cm}}^{(j)}$ across the different $j\in[1,L]$ is the typical signature of a structureless, homogeneous random configuration, see Fig. 14(d) and Fig.14(f). Therefore we introduce a measure of such dispersion as the variance of the angular centers of mass across the different slices,

[TABLE]

where we have defined

[TABLE]

for any arbitrary local observable $f_{j}$. We finally define the tomographic $x$-coherence as

[TABLE]

where the average $\langle\cdot\rangle_{\bm{\lambda}}$ is taken over the biased ${\bm{\lambda}}$-ensemble, i.e. over all trajectories statistically relevant during a rare event of fixed current-conjugated parameter ${\bm{\lambda}}$. We can define in an equivalent way the tomographic $y$-coherence $\Delta_{y}({\bm{\lambda}})$ to detect particle jams along the $y$-direction, and Fig. 13 shows these two order parameters measured across the DPT as a function of $z=|{\bm{\lambda}}+{\bf{E}}|$.

The results obtained with these two complementary order parameters fully confirm the presence of the DPT, shedding additional light on the symmetry-breaking process. In particular, $\Delta_{x}(z)$ increases steeply for ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}\leq\Xi_{c}$ and all angles $\phi$ of the current vector, while $\Delta_{y}(z)$ remains small and does not change appreciably across the DPT. This phenomenology clearly indicates that a coherent particle jam emerges along the $x$-direction for ${\bf{z}}\cdot{\hat{\cal A}}{\bf{z}}\leq\Xi_{c}$ and all angles, as in the sketch of Fig. 14(a) above. This means that only one of the two possible symmetry-broken density-wave phases appear in our simulations (regardless of the current vector orientation). Such result agrees with the MFT prediction in the supercritical anisotropy regime $\epsilon>\epsilon_{c}$, see Fig. 11(c), and is consistent with the measured effective anisotropy $\epsilon\approx 0.038>\epsilon_{c}$, see inset in Fig. 13. In addition, the behavior of both tomographic coherences $\Delta_{\alpha}$ ($\alpha=x,y$) across the DPT is consistent with the emergence of a traveling density-wave with structure in $1d$ and not in $2d$, as in the latter case both $\Delta_{\alpha}$ should increase upon crossing the critical line. Moreover, the acute but continuous change of $\Delta_{x}({\bf{z}})$ across the DPT is consistent with a second-order transition, in agreement with the MFT prediction.

VI Particle-hole symmetry breaking at the fluctuation level

The dynamical phase transitions in current statistics reported so far take place in periodic settings. In such systems, the additivity principle (which conjectures the time-independence of the dominant trajectories responsible of a fluctuation, see Section §III) is violated via the emergence of traveling-wave density profiles that break the time-translational invariance. For several years it was not clear whether other types of DPTs could appear for open diffusive systems, where energy or particles are injected through reservoirs attached to the boundaries. However, it was recently shown that first- and second-order DPTs might also occur for open systems with an external field through a novel mechanism: particle-hole (PH) symmetry breaking Baek et al. (2017, 2018). This important prediction, based on a perturbative Landau theory derived within MFT, is however restricted to (a) fluctuations around the critical current and (b) equal boundary densities/temperatures (or at most infinitesimally small gradients). Additionally, its microscopic origin is not understood yet and, most importantly, such a DPT has never been observed in numerical experiments, which might offer clues on novel phenomenology far from the critical point not yet explored. To shed light on all these issues, we have thoroughly studied the open $1d$ WASEP using both MFT and cloning Monte Carlo simulations in search of this elusive DPT Pérez-Espigares et al. (2018b).

The open $1d$ WASEP, already described in Section §IV.2, is characterized at the microscopic level by biased bulk jump rates $p^{\pm}=pe^{\pm E/L}$ to the right ($+$) and to the left ($-$), see Fig. 8(a). Taking $p=1/2$ for simplicity, this model is characterized at the macroscopic level by a diffusivity $D(\rho)=1/2$ and a mobility $\sigma(\rho)=\rho(1-\rho)$. For these transport coefficients, the current large deviation function $G(q)$ –given by Eq. (9)– remains invariant under a PH transformation, namely $\rho\to 1-\rho$ and $x\to 1-x$, whenever the densities at the boundaries satisfy the condition $\rho_{\text{L}}=1-\rho_{\text{R}}$. Interestingly, the work in Baek et al. (2017, 2018) predicts that this PH symmetry can be eventually broken when the system is conditioned to maintain atypically low currents in the presence of a large enough field ($|E|>E_{c}=\pi$). More in detail, in order to sustain a current $q$ above a critical threshold $(|q|>|q_{c}|)$, the system adopts an optimal profile $\rho_{q}(x)$ that is PH-symmetric, i.e. $\rho_{q}(x)=1-\rho_{q}(1-x)$, so the optimal profile inherits the symmetry of the governing action, the current LDF (9). This is illustrated by the black density profiles of the inset to Fig. 15, corresponding to the currents given by the black symbols in the main panel. Nevertheless, MFT predicts that for currents below the critical point $(|q|<|q_{c}|)$ the PH symmetry of the current LDF is broken: two coexisting symmetry-broken optimal profiles $\rho_{q}^{\pm}(x)$ appear, such that $\rho_{q}^{\pm}(x)\neq 1-\rho_{q}^{\pm}(1-x)$. These two different profiles are however linked by the PH transformation, i.e. $\rho_{q}^{\pm}(x)=1-\rho_{q}^{\mp}(1-x)$, thus restoring the broken symmetry. This is displayed in the inset to Fig. 15 by the red, $\rho_{q}^{+}(x)$, and blue profiles, $\rho_{q}^{-}(x)$, which are respectively associated with the currents given by the red and blue symbols of the main panel. As argued in Pérez-Espigares et al. (2018b), this symmetry-breaking phenomenon is captured by a global order parameter such as the total mass of the system $m=\int_{0}^{1}\rho(x)dx$. In particular, by studying the joint fluctuations of the current $q$ and this collective order parameter $m$, the dynamical phase diagram for arbitrary boundary gradients (both symmetric and asymmetric) can be unveiled. This is shown in the main panel of Fig. 15, which displays the mass $m_{q}$ of the optimal trajectory responsible for a current fluctuation $q$ for different boundary drivings. Note the $\mathbb{Z}_{2}$-type symmetry-breaking transition appearing at a critical current $|q_{c}|$ for $\rho_{\text{R}}=1-\rho_{\text{L}}$.

The DPT here reported for the open WASEP is, as its periodic counterpart (see Section IV.2), of second-order type and takes place whenever the densities at the boundaries are PH symmetric, $\rho_{\text{L}}=1-\rho_{\text{R}}$. As in standard critical phenomena, much information of the transition can be obtained by studying fluctuations of the order parameter (i.e. the total mass) across the critical point. In this case this means investigating the the joint mass-current LDF $G(m,q)$, or rather the conditional LDF $G(m|q)=G(m,q)-G(q)$. This can be computed from (9) just by minimizing over density profiles with the additional constraint $m=\int_{0}^{1}\rho(x)dx$. Fig. 16 illustrates the MFT prediction for $-G(m|q)$ and different boundary drivings. Note in particular that the optimal mass $m_{q}$ for a current fluctuation $q$, shown in the main panel of Fig. 15, is nothing but the mass minimizing $-G(m|q)$ for each $q$. Fig. 16 also shows the optimal profiles $\rho_{m,q}(x)$, derived from the MFT, associated with each mass $m$ for different currents $q$ (top and bottom rows) and different boundary drivings with $E=4$. For $\rho_{\text{L}}=1-\rho_{\text{R}}$, the joint LDF $-G(m|q)$ moves from a region with one global minimum to a symmetry-broken phase with two global minima as $|q|$ crosses a critical threshold $|q_{c}|$, a behavior typical of a $\mathbb{Z}_{2}$-type transition, and this is accompanied by an emergent degeneration in the optimal profiles for $|q|<|q_{c}|$, $\rho_{q}^{\pm}(x)=1-\rho_{q}^{\mp}(1-x)$, see the top and bottom rows of Figs. 16(a)-(b). On the other hand, for asymmetric boundary drivings, $\rho_{\text{L}}\neq 1-\rho_{\text{R}}$, there is always one single global minimum in $-G(m|q)$ $\forall q$, although a local minimum might appear for low enough currents associated with dynamical metastability, see Fig. 16(c).

A natural question that arises is whether a time-dependent density profile might become a better minimizer of the LDF $G(m|q)$ for some values of $m$. Interestingly, the answer is affirmative in the symmetry-broken phase, as a result of the non-convexity of $G(m|q)$ Pérez-Espigares et al. (2018b). In this case, a dynamical time-coexistence between the upper- and lower-branch optimal profiles $\rho_{q}^{\pm}(x)$ gives rise to the convex envelope of $G(m|q)$ for $m\in(m_{q}^{-},m_{q}^{+})$, with $m_{q}^{\pm}=\int_{0}^{1}\rho_{q}^{\pm}(x)dx$. This means that in order to produce a mass $m\in(m_{q}^{-},m_{q}^{+})$ sustaining a current $q$ during a long time interval $\tau$, the system adopts optimal profile is $\rho_{q}^{+}(x)$ for a fraction of time $\nu\tau$ with $\nu\in[0,1]$ while it maintains $\rho_{q}^{-}(x)$ for the rest of the time $(1-\nu)\tau$. This coexistence in time results in a convex envelope –as a Maxwell-like construction– which reads $G(m|q)=\nu G(m_{q}^{-}|q)+(1-\nu)G(m_{q}^{+}|q)$, thus maximizing the probability $P(m|q)$, as in a standard first-order phase transition Pérez-Espigares et al. (2018b).

So far we have analyzed the DPT at a macroscopic scale, however, in order to have a better understanding of the phenomenon we set out to study its origin from the microscopic dynamics. This approach may be carried out via exact diagonalization of the generator of the dynamics and by means of cloning numerical simulations. The former analysis shows that the DPT here reported corresponds to a degeneracy of the ground state (the one associated with the leading eigenvalue) of the dynamical generator. It turns out that its spectral gap closes as $L$ increases, i.e. the sub-leading eigenvalue coalesces with the leading one. For large but finite $L$, this gap is small but non-zero, thus giving rise to two long-lived metastable states (MS) which converge to the macroscopic symmetry-broken profiles for increasing sizes [see density profiles with square and triangle symbols in Figs. 17(a)-(b)]. Details on the derivation of these MS optimal profiles, obtained using the so-called Doob’s h-transform, can be found in Ref.Pérez-Espigares et al. (2018b).

Despite the physical understanding resulting from the spectral analysis, a direct observation of this phenomenon –occurring in the far tails of the current distribution–, remains to be done. For that reason we have carried out numerical simulations using the cloning Monte Carlo algorithm described in previous sections. This method allows us to reach larger system sizes than with the diagonalization method and to confirm the theoretical predictions for different boundary drivings. Results are shown in Figs. 17(a)-(b), where the optimal profiles in the symmetry broken phase (associated with $q=0$, thus corresponding to $\lambda=-E$ with $E=4$) have been displayed with purple down triangles both for (a) $\rho_{\text{L}}=\rho_{\text{R}}=0.5$ and (b) $\rho_{\text{L}}=0.8$ and $\rho_{\text{R}}=0.2$, using a population of $N_{c}=10^{4}$ clones. It is worth noting that, in order to correctly average the optimal profiles in the cloning simulations, we have to distinguish between those realizations sustaining a mass above $1/2$ and those sustaining a mass below $1/2$. Thus, averaging them separately, in order not to blur away the correct spatial structure of the density profiles, we obtain the results depicted in Fig. 17. Taking together the density profiles obtained from the spectral analysis of the microscopic generator and those obtained from the cloning Monte Carlo algorithm, we observe a clear convergence toward the theoretical MFT (symmetry-broken) prediction, though finite-size effects are still present. This suggest to pursue further cloning simulations introducing the recent advances of Section §III.5 to minimize sampling issues so as to reach larger system sizes with a reasonable number of clones. This will be the focus of future work.

VII Conclusions

In this paper we have reviewed a powerful Monte Carlo technique to sample rare events in many-particle systems. This computational method consists in modifying the system dynamics so that the rare events of interest become no longer rare, and involve the simultaneous evolution of multiple copies or clones of the system, which replicate or die in time according to their statistical weight. In particular, we have described in detail two different versions of this general method, namely as applied to discrete- and continuous-time stochastic many-particle systems. Special emphasis has been put on understanding the role of finite-size effects (both on the number of clones and the system size) on the sampling of the rare events of interest, as well as on the measurement and characterization of average trajectory observables characterizing a given rare event (midtime statistics).

The application of this new computational tool to simple models, particularly stochastic lattice gases, provides intriguing evidences of the existence of rich and fundamental structures in the fluctuating behavior of nonequilibrium systems, which typically emerge via dynamical phase transitions (DPTs). These DPTs, which are the focus of the second part of this work, appear when a system with many degrees of freedom sustains atypical fluctuations of dynamical observables such as the current or the activity. This leads in some cases to symmetry-broken space-time trajectories which enhance the probability of such events due to the emergence of ordered structures.

Despite their importance, and due to their low probability of occurrence, these DPTs are extremely difficult to characterize empirically. This makes the cloning Monte Carlo method an ideal tool to observe and characterize for the first time complex DPTs in many-particle systems. In this way, we have described some tricks of the trade related to the sampling of rare events across dynamical phase transition using the cloning method, with an emphasis on the definition and measurement of order parameters capturing the physics of the different DPTs, and the characterization of the optimal paths responsible for a fluctuation. In particular, we have described the application of the cloning method to uncover different DPTs in the current statistics of two paradigmatic models of transport, namely the Kipnis-Marchioro-Presutti (KMP) model of heat transport and the weakly asymmetric simple exclusion process (WASEP). These models exhibit spontaneous breaking of time-translation symmetry at the trajectory level under periodic boundary conditions, via the appearance of a time-dependent traveling wave. Despite their mathematical similarities, these DPTs have a radically different interpretation: while the DPT in the KMP model appears for large enough currents, where ballistic energy packets emerge to facilitate these fluctuations, in the WASEP the transition kicks in for low currents, where jammed density-wave trajectories dominate.

These DPTs are fully confirmed in numerical cloning experiments, which allow a detailed characterization of the emergence of order across the transition. Moreover, the versatility of the cloning Monte Carlo method has allowed us to explore the role of dimensionality in DPTs. In particular, we have investigated the vector current statistics in the $2d$ WASEP, where the complex interplay among an external field, the possible system anisotropy, and the vectorial character of currents leads to a rich phase diagram at the fluctuating level, with different symmetry-broken fluctuation phases separated by lines of first- and second-order DPTs. These predictions, based on macroscopic fluctuation theory, are fully confirmed in numerical cloning experiments by introducing novel order parameters to characterize the competing ordered phases which emerge for low currents. Finally, we have explored both theoretically and numerically a different type of symmetry-breaking phenomenon at the trajectory level which appears in open systems, i.e. coupled to boundary reservoirs which may impose an external gradient (of e.g. density or temperature). In this case a novel DPT appears which breaks the particle-hole symmetry –a $\mathbb{Z}_{2}$ discrete symmetry–, a transition that persists in the presence of arbitrarily strong (but symmetric) boundary gradients. The cloning method allows to observe and characterize for the first time this transition.

While the DPTs studied here appear in transport models of increasing complexity, we are still far from a complete understanding of the possible dynamical phases which may emerge in the fluctuations of more realistic systems, as e.g. hydrodynamic-type media characterized by several coupled and locally-conserved fields evolving in time in high dimensions. Our results above suggest that the resulting dynamical phase diagrams can be very rich, with multiple dynamical phases competing across different fluctuation regimes. The cloning method reviewed in this paper seems the ideal tool to explore this interesting phenomenology. However, given the scaling of the number of clones with the system size (including its dimension) needed to properly sample a given rare event, it seems likely that the standard cloning method will be insufficient to explore these events, and future work must focus on the application of the recent improvements of this method to compute rare event statistics in more realistic systems.

Acknowledgements.

We thank Robert Jack for insightful discussions. The research leading to these results has received funding from the Spanish Ministry MINECO project FIS2017-84256-P. Additionally, this study has been partially financed by the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), ref. SOMM17/6105/UGR.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. We are also grateful for the computational resources and assistance provided by PROTEUS, the super-computing center of Institute Carlos I in Granada, Spain.