# First-order dynamical phase transition in models of glasses: an approach based on ensembles of histories

Juan P. Garrahan, Robert L Jack, Vivien Lecomte, Estelle Pitard, Kristina

van Duijvendijk, Frédéric van Wijland

► To cite this version:

Juan P. Garrahan, Robert L Jack, Vivien Lecomte, Estelle Pitard, Kristina van Duijvendijk, et al. First-order dynamical phase transition in models of glasses: an approach based on ensembles of histories. Journal of Physics A: Mathematical and Theoretical, 2009, 42 (7), pp.075007. 10.1088/1751-8113/42/7/075007. hal-00364938

HAL Id: hal-00364938

https://hal.science/hal-00364938v1

Submitted on 21 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# First-order dynamical phase transition in models of glasses: an approach based on ensembles of histories

Juan P. Garrahan,¹ Robert L. Jack,²,³ Vivien Lecomte,⁴,⁵ Estelle

Pitard,⁶ Kristina van Duijvendijk,⁴ and Frédéric van Wijland⁴

¹School of Physics and Astronomy, University of Nottingham, Nottingham, NG7 2RD, UK

²Department of Chemistry, University of California, Berkeley, CA 94720-1460

³Department of Physics, University of Bath, Bath BA2 7AY, UK

⁴Laboratoire Matière et Systèmes Complexes (CNRS UMR 7057), Université Paris Diderot, 10 rue Alice Domon et Léonie Duquet, 75205 Paris cedex 13, France

⁵Département de Physique de la Matière Condensée, Université de Genève, 24 quai Ernest-Ansermet, 1211 Genève, Switzerland

⁶Laboratoire des Colloïdes, Verres et Nanomatériaux (CNRS UMR 5587), Université de Montpellier II, place Eugène Bataillon, 34095 Montpellier cedex 5, France

We investigate the dynamics of kinetically constrained models of glass formers by analysing the statistics of trajectories of the dynamics, or histories, using large deviation function methods. We show that, in general, these models exhibit a first-order dynamical transition between active and inactive dynamical phases. We argue that the dynamical heterogeneities displayed by these systems are a manifestation of dynamical first-order phase coexistence. In particular, we calculate dynamical large deviation functions, both analytically and numerically, for the Fredrickson-Andersen model, the East model, and constrained lattice gas models. We also show how large deviation functions can be obtained from a Landau-like theory for dynamical fluctuations. We discuss possibilities for similar dynamical phase-coexistence behaviour in other systems with heterogeneous dynamics.

## I. INTRODUCTION

In this paper we describe in detail a theoretical method for the study of the dynamics of glassy systems [1, 2, 3, 4]. This approach is in essence a statistical mechanics of the trajectories of the dynamics, or *histories*, as it is based on the study of large deviation functions [5]—which can be thought of as generalized free-energies—of dynamic observables. In particular, we use the tools of Ruelle’s thermodynamic formalism [6, 7], as applied to continuous time Markov chains [8], to study kinetically constrained models (KCMs) of glass formers [9]. In a recent letter [10] we showed using these methods that the dynamics of KCMs takes place on a first-order coexistence line between active and inactive dynamical phases, in accordance with previous suggestions [11]. Here we expand significantly on Ref. [10], demonstrating in detail the existence of the first-order dynamical phase transition, and discussing the Landau-like approach [12] that we use to characterise the dynamical phases, and the transition between them. The dynamical transition we find in KCMs [10] is related neither to a thermodynamic transition, nor to a finite temperature (or finite density) dynamical singularity. Our results, therefore, point towards a perspective [11] on glasses which is distinct from other approaches, such as the random first-order transition theory [13, 14, 15, 16], frustration-limited domains [17], or mode-coupling theory [18].

The paper is organized as follows: in Section II we introduce our dynamical tools and the ensemble of histories in which the dynamical phase transition takes place. In Section III we describe the models that we will consider. We show the existence of a dynamical phase transition in Section IV comparing different models and establishing minimal conditions that are sufficient to ensure a dynamical transition. In Section V we discuss the ensemble of histories in detail, considering statistical properties of the active and inactive phases, and a dynamical analogue of phase separation. We summarise our results in Section VI, and consider some open questions.

## II. DYNAMICAL TOOLS: THE s-ENSEMBLE

### A. Motivations

In this article, we are concerned with fluctuations in dynamical observables such as the amount of dynamical activity in a glassy system, integrated over a long time $t$ and over a large (but finite) system. To investigate these fluctuations, we consider statistical properties of the histories followed by the system. Ensembles of histories are central to the thermodynamic formalism developed by Ruelle and coworkers [6] (see [7] for a comprehensive review). While thermodynamics is concerned with probability distributions over configurations of a large system, we will apply
the thermodynamic formalism to probability distributions over histories. We begin by discussing the physical content of the observables that we will consider.

In the Boltzmann-Gibbs theory, the macroscopic features of large systems are characterised by determining the statistical properties (the mean value and fluctuations) of extensive observables, such as the energy or the number of particles. In a microcanonical approach, one considers the properties of a system with fixed total energy E. They are obtained from the counting factor

$$ \Omega(E, N) = \begin{cases} \text{number of configurations} \\ \text{with energy } E \end{cases} \qquad (1) $$

where N represents the size (the volume) of the system. In the large size limit ($N \to \infty$), we define the entropy density $s(e) = \lim_{N \to \infty} \frac{1}{N} \ln \Omega(eN, N)$, which represents the relative weight of configurations with energy density e.

In a dynamical context, we consider histories of the system between an initial time $\tau = 0$ and a final time $\tau = t$. Instead of considering the statistics of the energy E, we will consider an observable A, that is extensive in the observation time t. The dynamical analog of $\Omega(E, N)$ is the probability distribution of this observable

$$ \Omega_{\text{dyn}}(A, t) = \begin{cases} \text{fraction of histories with a given value of the} \\ \text{time-extensive observable } A \end{cases} \qquad (2) $$

On a mathematical level, the choice of the observable A is somewhat arbitrary, although application of the thermodynamic formalism requires that the quantity $\frac{1}{t} \log \Omega_{\text{dyn}}(at, t)$ should have a finite limit for large times t. Subject to this constraint, the choice of the order parameter A is informed by physical insight: we should use an observable that reveals the essential physical processes at work in the system. For example, in non-equilibrium systems in contact with two reservoirs of particles, we might define A as the total particle current: the number of particles transferred from one reservoir to the other between times 0 and t (see, for example, Refs [19, 20]). In the context of glassy phenomena, we consider observables that measure the “activity” or the “complexity” of the history [8, 10, 11].

Returning to the Boltzmann-Gibbs approach, it is useful to define the canonical ensemble through the partition function

$$ Z(\beta, N) = \sum_{E} \Omega(E, N) e^{-\beta E} \qquad (3) $$

which characterises a system at a given temperature $\beta^{-1}$. Within this framework, phase transitions can be identified from singularities in the intensive free energy, $f(\beta) = -\lim_{N \to \infty} \frac{1}{\beta N} \ln Z(\beta, N)$. The dynamical analog of this thermodynamic partition sum is

$$ Z_A(s, t) = \sum_A \Omega_{\text{dyn}}(A, t) e^{-sA} \qquad (4) $$

where we introduced an intensive field s conjugate to A. This field will play a role similar to the inverse temperature $\beta$. The dynamical partition function $Z_A(s, t)$ is the central object of Ruelle's thermodynamic formalism.

We have focused on the correspondence between the thermodynamic limit of large system size ($N \to \infty$) and the long time limit ($t \to \infty$) in Ruelle's formalism. In the following, we will consider systems for which the large time limit is to be taken at fixed system size: in some cases, we will then take a second limit of large system size N. If we consider systems with no thermodynamic phase transitions, then no singular behaviour arises on taking the limit of large N at fixed t. In this case, we expect the limits of large N and large t to commute, but this is clearly not the case in general.

## B. Systems with Markov dynamics: statistics over histories

### 1. Continuous time Markov evolution

We now give more precise definitions of the quantities discussed so far, by reviewing the construction of the ensemble of histories for stochastic systems. We focus on continuous-time Markov dynamics (in this section, we follow [8]). The system is defined by a finite set of configurations $\{C\}$. Its dynamical evolution is defined by the rates $W(C \to C')$ for transitions from configuration C to configuration C'. Thus, the probability $P(C, t)$ of being in configuration C at time t evolves according to a master equation:

$$ \partial_t P(C, t) = -r(C)P(C, t) + \sum_{C'} W(C' \to C)P(C', t) \qquad (5) $$
<figure><img src="image_2.png" /><figcaption>FIG. 1: A history of duration $t$ is defined by a sequence of configurations $C_0 \to \dots \to C_K$ and a sequence of jump times $t_1, \dots, t_K$. Between $t_k$ and $t_{k+1}$, the system stays in configuration $C_k$.</figcaption></figure>

where

$$r(C) = \sum_{C'} W(C \to C') \qquad (6)$$

represents the rate of escape from $C$. Equ.(5) is sufficiently general to describe kinetically constrained models, spin facilitated models, or lattice gases, with $C$ representing the configuration of the whole lattice in each case.

Starting from a configuration $C_0$ at initial time $t=0$, the system will experience a fluctuating number of changes of configuration (“jumps”) between 0 and $t$. We shall refer to the number of jumps as the “activity” and denote it by $K$. A history (or trajectory) consists of a sequence $C_0 \to \dots \to C_K$ of visited configurations, and a sequence of times $t_1, \dots, t_K$ at which the jumps occur (Fig. 1). We stress that for a fixed observation time $t$, the number of jumps is a fluctuating quantity: it depends on the particular history followed by the system between 0 and $t$. We refer to histories with many hops (large $K$) as ‘active’ histories and those with few hops (small $K$) as ‘inactive’.

We use the notation $\langle \mathcal{O} \rangle$ for an average of the observable $\mathcal{O}$, over histories of the system. We consider observables that depend on the entire history of the system, through the configurations visited and the time spent in each: that is, $\mathcal{O} = \mathcal{O}(C_0 \dots C_K, t_1 \dots t_K)$. In general, we have

$$\langle \mathcal{O} \rangle = \sum_K \sum_{C_0 \dots C_K} \int dt_1 \dots dt_K p_0(C_0) \left[ \prod_{k=1}^{K} W(C_{k-1} \to C_k) \right] \exp \left[ -\sum_{k=1}^{K} r(C_k)(t_{k+1} - t_k) \right] \mathcal{O}(C_0 \dots C_K, t_1 \dots t_K) \quad (7)$$

where the limits on the time integrals are $t_1 > 0$, $t_k < t_{k+1}$, and $t_K < t_{K+1} \equiv t$; we use $p_0(C_0)$ to denote the probability distribution of the initial configuration $C_0$. We use a compact notation for averages of this form:

$$\langle \mathcal{O} \rangle = \sum_{\text{hist}} \text{Prob}[\text{hist}] \mathcal{O}[\text{hist}]. \qquad (8)$$

where Prob[hist] plays the role of a probability density in the space of histories.

## 2. Time-extensive observables

Having defined our system and its histories, we now turn to the choice of the time-extensive observable $A$. A simple choice of this observable will be the activity $K$. Each time the system changes configuration $C \to C'$ the activity $K$ is incremented: $K \to K + 1$. More generally, we can consider an observable $A$ that is incremented at each jump, with the increment $\alpha(C, C')$ depending on the configurations before and after the jump. That is, for a given history with $K$ changes of configurations

$$A[\text{hist}] = \sum_{k=0}^{K-1} \alpha(C_k, C_{k+1}) \qquad (9)$$

Again, we note that if $\alpha(C, C') = 1$ then $A$ is the activity $K$.

To construct the dynamical partition sum, we start with a ‘microcanonical’ approach, classifying trajectories by their values of $A$. We generalise the probability $P(C, t)$, defining $P(C, A, t)$ as the probability of being in configuration $C$ at time $t$, having measured a value $A$ of the time-extensive variable between 0 and $t$. Its evolution in time is given by the master equation

$$\partial_t P(C, A, t) = \sum_{C'} W(C' \to C) P(C', A - \alpha(C', C), t) - r(C) P(C, A, t) \qquad (10)$$
Thus, the probability of measuring a value A for the observable A in a history of length t is

$$ \Omega_{\text{dyn}}(A, t) = \sum_{C} P(C, A, t) \qquad (11) $$

which we identify as the quantity introduced in (2).

### 3. Canonical description: evolution in the s-ensemble

We have defined the distribution $\Omega_{dyn}(A, t)$ that is the analog of the microcanonical counting factor $\Omega(E, N)$. We now introduce the analog for the canonical (Boltzmann-Gibbs) ensemble, parameterized by a field s. This involves a modification to the statistical weight of each history:

$$ \text{Prob}[\text{hist}] \rightarrow \text{Prob}[\text{hist}]e^{-sA[\text{hist}]} \qquad (12) $$

Thus, in the 's-ensemble', averages of observables $\mathcal{O}$ are given by

$$ \langle \mathcal{O} \rangle_s = \frac{1}{Z_A(s, t)} \sum_{\text{hist}} \mathcal{O}[\text{hist}] \text{Prob}[\text{hist}] e^{-sA[\text{hist}]} = \frac{\langle \mathcal{O} e^{-sA} \rangle}{\langle e^{-sA} \rangle} \qquad (13) $$

where

$$ Z_A(s, t) = \sum_{\text{hist}} \text{Prob}[\text{hist}]e^{-sA[\text{hist}]} = \langle e^{-sA} \rangle \qquad (14) $$

is the dynamical partition function, introduced in (4). (The subscript A of $Z_A$ serves as a reminder that the field s is conjugate to A.)

Averages in the ensemble with $s = 0$ correspond to the steady state averages of $\mathcal{O}$. A priori, this is the only physically accessible ensemble. Positive or negative values of s favor histories with non-typical values of A. For our purposes, working in the s-ensemble is simpler than considering ensembles with fixed values of A. We take the Laplace transform of $P(C, A, t)$ with respect to A:

$$ \hat{P}_A(C, s, t) = \sum_A e^{-sA} P(C, A, t) \qquad (15) $$

From (10), the equation of motion for $\hat{P}_A(C, s, t)$ is

$$ \partial_t \hat{P}_A(C, s, t) = \sum_{C'} e^{-s\alpha(C', C)} W(C' \rightarrow C) \hat{P}_A(C', s, t) - r(C) \hat{P}_A(C, s, t), \qquad (16) $$

or, in an operator notation, $\partial_t \hat{P}_A = W_A \hat{P}_A$, where $W_A$ operates in the space of configurations $\{C\}$. Its matrix elements are

$$ (W_A)_{C,C'} = W(C' \rightarrow C)e^{-s\alpha(C',C)} - r(C)\delta_{C,C'}. \qquad (17) $$

Some properties of the operator W are discussed in appendix A: Equ. (A2) states that $\hat{P}_A(C, s, t)$ behaves in the large time limit as $\hat{P}_A(C, s, t) \sim R_0(C, s)e^{t\psi_A(s)}$ where $\psi_A(s)$ is the largest eigenvalue of $W_A$ and $R_0(C, s)$ is the associated eigenvector. Thus, for large times,

$$ Z_A(s, t) = \sum_C \hat{P}_A(C, s, t) \sim e^{t\psi_A(s)}, \qquad (18) $$

and we will refer to $\psi_A(s)$ as (the negative of) the dynamical free energy per unit time. Summing Eq. (15) over $C$, probability conservation implies $Z_A(0, t) = 1$, so that $\psi(0) = 0$ for all stochastic systems.
### 4. Large deviation functions

In the Boltzmann-Gibbs theory, entropy and free energy are related through a Legendre transform (as can be seen from (3) or [21]) which provides a link between microcanonical and canonical ensembles. We have already defined the function $\psi_K(s)$, which is the dynamical analog of the free energy density $f(\beta)$. The dynamical analog of the entropy density $s(e)$ is

$$ \pi(a) = \lim_{t \to \infty} \frac{1}{t} \log \Omega_{\text{dyn}}(at, t) \qquad (19) $$

which determines the large-$t$ scaling of the probability of observing a value at for the observable A.

For large times, the sum in (4) is dominated by the maximum of $\Omega_{\text{dyn}}(A, t)$, so that $\pi(a)$ and $\psi_A(s)$ are related through a Legendre transform:

$$ \psi_A(s) = \max_a (\pi(a) - sa) \qquad (20) $$

If the function $\pi(a)$ satisfies $\pi''(a) \le 0$, it can be obtained from the inverse transform

$$ \pi(a) = \min_s (\psi_A(s) + sa) \qquad (21) $$

Physically, the quantity $\pi(a)$ describes the large fluctuations of A. It is maximal at the most probable value of a, which is the mean value of $A/t$, in the limit of large time $t$. Gaussian fluctuations of $A/t$ are described by the quadratic approximation of $\pi(a)$ around its maximum. Expanding $\pi(a)$ beyond quadratic order gives information about non-Gaussian fluctuations of $A/t$, which are referred to as large deviations [5]. Alternatively, one may characterise these fluctuations through $\psi_A(s)$, since the cumulants of A are obtained from the derivatives of $\psi_A(s)$ through $\lim_{t \to \infty} \frac{1}{t} \langle A^p \rangle_c = (-1)^p \frac{\partial^p \psi_A(s)}{ds^p} \Big|_{s=0}$, where, as usual, $\langle A^p \rangle_c$ is the p-th cumulant of A.

### 5. Time-extensive observables varying continuously in time

In addition to time-extensive order parameters of the form given in (9), we also consider those of the form

$$ B[\text{hist}] = \sum_{k=0}^{K} (t_{k+1} - t_k) b(C_K) = \int_{0}^{t} dt' b(C(t')), \qquad (22) $$

where we introduced a configuration-dependent observable $b(C)$. In the sum over $k$, we define $t_0 = 0$ and $t_{K+1} = t$ so that the time spent in configuration $C_k$ is simply $t_{k+1} - t_k$. In the integral representation, we have represented the trajectory by a function $C(t')$ which takes the value $C_k$ for $t_k < t' < t_{k+1}$. The time-integrated energy of the system is an observable of the form $B$, in which case $b(C)$ is simply the energy of configuration $C$. Then, defining $P(C, B, t)$ by analogy with $P(C, A, t)$, we have

$$ \partial_t P(C, B, t) = \sum_{C'} W(C' \to C) P(C', B, t) - r(C) P(C, B, t) - b(C) \frac{\partial}{\partial B} P(C, B, t) \qquad (23) $$

We define an s-ensemble associated with the observable B through

$$ \text{Prob}[\text{hist}] \to \text{Prob}[\text{hist}]e^{-sB[\text{hist}]} \qquad (24) $$

Then, repeating the analysis of Section II B 3, the analog of $\psi_A(s)$ is $\phi_B(s) = \lim_{t \to \infty} \frac{1}{t} \ln\langle e^{-sB} \rangle$. This quantity is equal to the maximal eigenvalue of an operator $W_B$, whose elements are

$$ (W_B)_{C,C'} = W(C' \to C) - [r(C) + sb(C)] \delta_{C,C'} \qquad (25) $$

In the following, we concentrate our study on time-extensive variables of type A. Some connections between s-ensembles parameterized by observables of types A and B discussed in Appendix B
## 6. Variational approach for $\psi_A(s)$

The models considered in this work have dynamics which obey detailed balance with respect to an equilibrium distribution $P_{\text{eq}}(C)$: that is, $P_{\text{eq}}(C)W(C \to C') = P_{\text{eq}}(C')W(C' \to C)$. This allows us to derive a variational bound on the dynamical free energy $\psi_K(s)$. To achieve this, we symmetrise the evolution operator $\tilde{W}_K$, defining $\tilde{W}_K$ through the similarity transformation $(\tilde{W}_K)_{C,C'} = P_{\text{eq}}^{-1/2}(C)(\tilde{W}_K)_{C,C'}P_{\text{eq}}^{1/2}(C')$. Hence,

$$ (\tilde{W}_K)_{C'C} = e^{-s} [W(C \to C')W(C' \to C)]^{\frac{1}{2}} - r(C)\delta_{CC'} = (\tilde{W}_K)_{CC'} \quad (26) $$

Since $\tilde{W}_K$ and $W_K$ are related by a similarity transformation, their eigenspectra are identical. We therefore use a variational principle (valid for any symmetric operator) to determine their common maximal eigenvalue:

$$ \psi_K(s) = \max_{\{V(C)\}} \frac{\sum_{C,C'} V(C)(\tilde{W}_K)_{C,C'} V(C')}{\sum_C V(C)^2} = \max_{|V\rangle} \frac{\langle V | \tilde{W}_K | V \rangle}{\langle V | V \rangle} \quad (27) $$

At $s=0$, the maximum is achieved for $V(C) = P_{\text{eq}}(C)^{1/2}$, and $\psi(0) = 0$, as required.

Interestingly, the quantity to be maximised in (27) has a physical interpretation. For any history of the system, the fraction of time spent in each configuration $C$ defines a quantity known as the experimental measure. As we discuss in appendix C, Donsker-Varadhan theory relates the probability of observing a particular experimental measure to an expectation value of the form $\langle V|\tilde{W}_K|V \rangle$. In Section V, we will use these results to investigate fluctuations in the $s$-ensemble.

# III. MODELS AND ORDER PARAMETERS

## A. Kinetically constrained models: FA, East, TLG and KA models

Kinetically constrained models [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34] are simple lattice models of glasses which can account for a large range of dynamical phenomena associated to the glass transition. This includes: super-Arrhenius temperature dependence of timescales, non-exponential relaxation, spatially heterogeneous dynamics, transport decoupling, and aging and memory effects. The thermodynamic properties of KCMs are very simple, and their non-trivial features arise from dynamical rules which forbid or favor some transitions, while maintaining detailed balance with respect to a trivial equilibrium distribution over configurations. For a review on KCMs see [9].

We first consider models with binary spins $n_i = 0, 1$ where $i = 1, \dots, N$ are the sites of a lattice. In spin-facilitated models, sites with $n_i = 1$ represent excitations, which promote local activity. The models evolve by single spin-flips, which occur with rates

$$ W(n_i \to 1 - n_i) = C_i(\{n_j\}) \frac{e^{\beta(n_i-1)}}{1 + e^{-\beta}} \quad (28) $$

where $\beta$ is the inverse temperature, and the kinetic constraint enters through the function $C_i(\{n_j\})$, which is a function of the neighbors $n_j$ of $i$, but does not itself depend on $n_i$. In this case, it is simple to verify that the model obeys detailed balance with respect to the equilibrium distribution

$$ P_{\text{eq}}(\{n_i\}) = \prod_i \frac{e^{-\beta n_i}}{1 + e^{-\beta}} \quad (29) $$

and that the excitation density is $c \equiv \langle n_i \rangle = (1 + e^{\beta})^{-1}$.

In the one-spin facilitated Fredrickson-Andersen (FA) model [9, 22], $C_i = 1$ if any of the nearest neighbors $j$ of $i$ are in the excited state, $n_j = 1$; otherwise $C_i = 0$. We also consider the three-dimensional variant of the East model [23, 27, 35] in which $C_i = 1$ for site $i = (x, y, z)$ if at least one of the sites $(x-1, y, z), (x, y-1, z)$ or $(x, y, z-1)$ is in the excited state; otherwise $C_i = 0$.

In addition, we consider lattice gas models [9, 24, 25], in which particles move from site to site, with at most one particle per site. Sites which are occupied have $n_i = 1$, and unoccupied sites have $n_i = 0$. Particles move between sites $i$ and $j$ with rate $C_{ij}(\{n_k\})$ so that the model has a conserved density $\rho = N^{-1} \sum_i n_i$. The rate $C_{ij}(\{n_k\})$ is non-zero only for nearest neighbor sites $i$ and $j$, and it is independent of $n_i$ and $n_j$. Thus, equilibrium state has a trivial distribution: all configurations with density $\rho$ have equal probability. As an example of such a model, we
consider the two-vacancy facilitated triangular lattice gas, or (2)-TLG [25], which is defined on a triangular lattice, with a constraint $C_{ij}$ which is equal to unity if the two common nearest neighbors of sites $i$ and $j$ are vacant, and zero otherwise. Similarly, the (2,2) variant of the Kob-Andersen (KA) lattice model [24] is defined on a square lattice, with $C_{ij} = 1$ if at least one neighbour $k \neq j$ of site $i$ has $n_k = 0$ and at least one neighbour $k' \neq i$ of site $j$ has $n_{k'} = 0$. Otherwise $C_{ij} = 0$.

## B. Reducibility of KCMs and sums over histories

The construction of the s-ensemble in Section II assumed that the system of interest has a single steady state to which it converges in the long-time limit. For finite-sized stochastic systems, this convergence is ensured as long as the dynamics are irreducible [9]: that is, it must be possible for every configuration of the system to be reached from every other configuration. For KCMs, this is not the case in general. For example, in the FA model, there are no transitions either into or out of the configuration with no excited sites ($n_i = 0$ for all $i$). For the models considered here, these states are usually considered to be irrelevant because they have a contribution to the Gibbs measure that vanishes exponentially in the thermodynamic limit, at all temperatures $T > 0$.

However, when considering large deviations, these states may become relevant. In order to ensure convergence to a single steady state, we define our unbiased measure over histories as in (7), with a distribution of initial conditions $p_0(C_0)$ that is non-zero only for configurations in the largest irreducible partition of the dynamics. That is, we do not allow the system to occupy configurations that cannot be reached from representative configurations taken from the relevant Gibbs ensemble. For the FA model and East models, this simply means that the system may not occupy the configuration which has no excited sites. Practically, this means (for example) that the maximum in the variational expression (27) should be taken with the constraint that $V(C)$ is finite only for configurations in the largest irreducible partition.

Instead of restricting initial conditions in this way, one could instead consider large deviations in a subsystem of size $N$ that is embedded in a larger system of size $N' \gg N$: this was the approach taken in [11]. As usual, we expect these approaches to be equivalent in the limit of large system size $N$.

## C. Kinetically constrained models : bosonic and mean-field variants

It is convenient to define a 'bosonic' variant of the one-spin facilitated FA model [36, 37], in which $n_i$ may be any integer greater than or equal to zero. We take

$$W(n_i \rightarrow n_i + 1) = C_i(\{n_j\})e^{-\beta}, \quad W(n_i \rightarrow n_i - 1) = C_i(\{n_j\})n_i \quad (30)$$

where $C_i$ is again independent of $n_i$ so that

$$P_{\text{eq}}(\{n_i\}) = \prod_i \frac{c^{n_i}}{n_i!} e^{-c} \quad (31)$$

where $c \equiv \langle n_i \rangle = e^{-\beta}$. For the bosonic FA model in finite dimension, we take $C_i(\{n_j\}) = \sum_{\langle j \rangle} n_j$, where the sum is over the nearest neighbours $j$ of site $i$.

For this bosonic model, it is convenient to use the Doi-Peliti representation [38]. We define bosonic operators $a_i$, $a_i^\dagger$ and $\hat{n}_i = a_i^\dagger a_i$, with $[a_i, a_j^\dagger] = \delta_{ij}$ and a vacuum state $|0\rangle$ for which $a_i|0\rangle = 0$ for all $i$. The Doi-Peliti representation of the operator $\mathcal{W}$ is defined by $(\mathcal{W})_{CC'} = \langle 0 | \left[ \prod_i \frac{a_i^{n_i}}{n_i!} \right] \mathcal{W} \prod_i (a_i^\dagger)^{n_i'} |0\rangle$ where the configurations $C$ and $C'$ have occupations $\{n_i\}$ and $\{n_i'\}$ respectively. In the s-ensemble, we are interested in the operator $\mathcal{W}_A$ defined in (17). In the case where the observable $A$ is the activity $K$, we have

$$\mathcal{W}_K^{(\text{FA})} = \sum_{\langle ij \rangle} \left[ e^{-s}(ca_i^\dagger + a_i)\hat{n}_j + e^{-s}(ca_j^\dagger + a_j)\hat{n}_i - 2(\hat{n}_i\hat{n}_j + c) \right] \quad (32)$$

where the sum runs over (distinct) pairs of nearest neighbours.

In addition, it is useful to consider a mean-field variant of the FA model, in which the facilitation function of site $i$ depends symmetrically on the state of all sites. That is, $W(n_i \rightarrow n_i + 1) = N^{-1} \sum_j n_j e^{-\beta}$, and $W(n_i + 1 \rightarrow n_i) = N^{-1} \sum_j n_j n_i$, which satisfy detailed balance with respect to (31). In the Doi-Peliti representation, the master
operator is simply

$$ \mathbb{W}_K^{(\text{FA,mf})} = (2N)^{-1} \sum_{ij} \left[ e^{-s} a_i^\dagger (ca_j^\dagger + a_j) a_i + e^{-s} a_j^\dagger (ca_i^\dagger + a_i) a_j - 2(a_i^\dagger a_j^\dagger a_j a_i + c) \right] \quad (33) $$

Due to the symmetry between sites, the properties of the model can be obtained from a single co-ordinate: the total number of excitations $n_{tot} = \sum_i n_i$, whose equilibrium distribution is Poissonian with mean $cN$. In this co-ordinate, the master-like equation (16) has a closed form, and the matrix elements of the relevant operator are simply

$$ (\mathbb{W}_K^{(\text{FA,mf})})_{n'_{\text{tot}}, n_{\text{tot}}} = cn_{\text{tot}}(e^{-s}\delta_{n_{\text{tot}}+1, n'_{\text{tot}}} - \delta_{n_{\text{tot}}, n'_{\text{tot}}}) + \frac{n_{\text{tot}}}{N}(n_{\text{tot}}-1)(e^{-s}\delta_{n_{\text{tot}}-1, n'_{\text{tot}}} - \delta_{n_{\text{tot}}, n'_{\text{tot}}}) \quad (34) $$

## D. The A-model and the AA model

It will be useful to compare the FA model with two other models, which we call the A and AA models. These names are motivated by the schematic representations of their fundamental processes, as $A \leftrightarrow \emptyset$ and $A+A \leftrightarrow \emptyset$. Here we have used an alternative notation to avoid confusion with the observable A used to define the s-ensemble.

We define the A-model and its bosonic variant by removing the kinetic constraints from the FA model: that is, $C_i(\{n_j\}) = 1$, independent of the state of the system. In this model, excitations are created and destroyed singly, independent of site. The A-model has the same equilibrium distribution as the FA model, but its large deviations can be solved exactly.

We also compare the FA model with a model in which particles appear and annihilate (AA) in pairs [39]. This so-called AA model is related to a variant of the FA model, through a similarity transformation that connects their master operators [37]. The AA model is defined for binary spins $n_i = 0, 1$. In this model, the excitations move between adjacent sites with rate $D$, and appear and annihilate in pairs with rates $k$ and $k'$ respectively. Schematically, we write

$$ 0_i 1_j \stackrel{D}{\leftrightarrow} 0_i 0_j, \quad 1_i 1_j \stackrel{k}{\rightarrow} 0_i 0_j, \quad 0_i 0_j \stackrel{k'}{\rightarrow} 1_i 1_j \quad (35) $$

for neighbouring sites $i$ and $j$. The equilibrium state of this model is of the form (29), with $e^{-\beta} = \sqrt{k'/k}$. It is also convenient to consider bosonic and mean-field variants of this model, defined analogously to their FA counterparts. In the bosonic AA model, we generalise to $n_i \ge 0$, using rates

$$ (n_i, n_j) \stackrel{D n_i}{\rightarrow} (n_i - 1, n_j + 1), \quad (n_i, n_j) \stackrel{k n_i n_j}{\rightarrow} (n_i - 1, n_j - 1), \quad (n_i, n_j) \stackrel{k'}{\rightarrow} (n_i + 1, n_j + 1) \quad (36) $$

which obey detailed balance with respect to (31) with $k' = ke^{-2\beta}$ as before. In the Doi-Peliti formalism, we have

$$ \mathbb{W}_K^{(\text{AA})} = \sum_{\langle ij \rangle} e^{-s} [k' a_i^\dagger a_j^\dagger + k a_i a_j + D(a_i^\dagger a_j + a_j^\dagger a_i)] - [k \hat{n}_i \hat{n}_j + ck' + D(\hat{n}_i + \hat{n}_j)] \quad (37) $$

In the mean field variant of the AA model, diffusion occurs between all pairs of sites ($i \neq j$), with rate $(D/N)$; pair creation and annihilation processes occur with rates $(k/N)n_i n_j$ and $(k'/N)$ for all pairs of sites $i \neq j$; and we also allow for on-site pair creation and annihilation: $n_i \rightarrow n_i \pm 2$ with rates $(k'/N)$ and $(k/N)n_i(n_i - 1)$. In the Doi-Peliti representation, the master operator is

$$ \mathbb{W}_K^{(\text{AA,mf})} = \frac{1}{2N} \sum_{ij} \left[ e^{-s} (k' a_i^\dagger a_j^\dagger + k a_i a_j) - (k a_i^\dagger a_j^\dagger a_j a_i + ck') \right] + \frac{D}{2N} \sum_{i \neq j} \left[ e^{-s} (a_i^\dagger a_j + a_j^\dagger a_i) - (\hat{n}_i + \hat{n}_j) \right] \quad (38) $$

For finite systems, the restriction to $i \neq j$ in the diffusion term means that the master-like equation cannot be written in terms of the single co-ordinate $n_{tot}$, except at $s = 0$. However, in the limit of large-$N$, this single co-ordinate is sufficient, and the master-like operator for this co-ordinate reduces to

$$ (\mathbb{W}_K^{(\text{AA,mf})})_{n'_{\text{tot}}, n_{\text{tot}}} = k' N (z \delta_{n_{\text{tot}}+2, n'_{\text{tot}}} - \delta_{n_{\text{tot}}, n'_{\text{tot}}}) + k \frac{n_{\text{tot}}}{N} (n_{\text{tot}} - 1)(z \delta_{n_{\text{tot}}-2, n'_{\text{tot}}} - \delta_{n_{\text{tot}}, n'_{\text{tot}}}) + D(z-1)n_{\text{tot}} \quad (39) $$

with $z = e^{-s}$.

## E. Relevant observables

We now discuss the observables that we will use to define the s-ensemble, and those that we will use to characterise trajectories within that ensemble.
1. The activity K and the complexity $Q_+$

We have already defined the activity $K$, which counts the number of changes of configuration in a dynamical trajectory. In the context of dynamically heterogeneous systems such as glass-formers, the local activity can be used to distinguish mobile and immobile regions of the system. The large deviations of the extensive activity $K$ are used to characterise trajectories which are more or less mobile than average.

We note that $K$ is of the form given in (9) with $\alpha(C', C) = 1$, so the properties of the relevant $s$-ensemble are encoded in the operator

$$ (\mathbb{W}_K)_{C,C'} = e^{-s} W(C' \rightarrow C) - r(C) \delta_{C,C'} \quad (40) $$

Systems with dynamical heterogeneities are likely to present a wide distribution of very different histories. One way of characterizing this diversity is provided by the dynamical complexity of the histories [11, 40, 41]. In the context of dynamical system theory, this quantity is called the Kolmogorov-Sinai entropy [6]. It provides one with the information content of the history and is defined as the logarithm of the probability of the history. As discussed in [8], the appropriate generalisation of this approach to systems with Markov dynamics is to consider the entropy associated with the measure over sequences of configurations $C_0 \rightarrow \dots \rightarrow C_K$ [8]. This amounts to performing a coarse-graining in time: it means that the information associated with the time intervals between changes of configuration is ignored when calculating the complexity. The definition of the dynamical complexity is

$$ Q_+ = \sum_{k=0}^{K-1} \ln \frac{W(C_k \rightarrow C_{k+1})}{r(C_k)}, \quad (41) $$

which is of the form given in (9). Thus, we define a dynamical partition sum

$$ Z_+(s, t) = \langle e^{-sQ_+} \rangle. \quad (42) $$

The corresponding dynamical free energy is $\psi_+(s) = \lim_{t \to \infty} \frac{1}{t} \ln Z_+(s, t)$ which corresponds to the topological pressure of dynamical system theory. The analog of the Kolmogorov-Sinai entropy $h_{KS}$ is

$$ h_{KS} = - \lim_{t \to \infty} \frac{1}{t} \langle Q_+ \rangle = \frac{d}{ds} \psi_+(s) \quad (43) $$

which provides a measurement of the dynamical complexity of the histories in the steady state. In the examples of glass formers we will study below, the dynamical ensembles given by $K$ and $Q_+$ are qualitatively similar: we concentrate on the activity $K$ for simplicity.

2. Fluctuation theorem in the $s$-ensemble

The Gallavotti–Cohen relation holds also in the quasi-stationary state at fixed value of $K$ (or $s$), and therefore the fluctuation–dissipation theorem holds there as well. In order to see this, we parallel the reasoning presented by Lebowitz and Spohn in [42], and we construct the operator governing the dynamics not only at fixed value of the activity $K$ but also at fixed value of the entropy current $Q_S = \sum_{n=0}^{K-1} \ln \frac{W(C_n \rightarrow C_{n+1})}{W(C_{n+1} \rightarrow C_n)}$, which, in terms of the variables $s$ and $\lambda$ conjugate to the activity $K$ and the entropy current $Q_S$ respectively, leads to the following pseudo–evolution operator,

$$ (\mathbb{W}(s, \lambda))_{C,C'} = e^{-s} W(C' \rightarrow C)^{1-\lambda} W(C' \rightarrow C)^{\lambda} - r(C) \delta_{C,C'} \quad (44) $$

whose property $\mathbb{W}(s, \lambda)^{\dagger} = \mathbb{W}(s, 1-\lambda)$ ensures that its largest eigenvalue $\psi$ verifies $\psi(s, \lambda) = \psi(s, 1-\lambda)$. For system with particle conservation, and subject to a field driving the system out of equilibrium, we note that the entropy current $Q_S$ is directly proportional to the total current of particles flowing through the system [42]. In that case, the generalized symmetry $\psi(s, \lambda) = \psi(s, 1-\lambda)$ implies a fluctuation-dissipation like relation in the $s$-ensemble.

3. Order parameters within the $s$-ensemble

As well as using the observables $K$ and $Q_+$ to define $s$-ensembles through (12), we also characterise the $s$-ensemble by using two other order parameters. For spin-facilitated models, we consider the mean excitation density:

$$ \rho_K(s) = \lim_{t \to \infty} \frac{1}{N} \left\langle \int_0^t d\tau \sum_i n_i(\tau) \right\rangle_s \quad (45) $$
For lattice gas models, the particle density is specified by the initial conditions, so we require a different order parameter. The average activity is given by $\frac{1}{Nt}\langle K \rangle_s$. One can also consider the mean escape rate $r(C)$ which depends only on configurations of the system. Again, we time-average this quantity along the trajectories, and divide by the system size $N$, defining

$$r_K(s) = \lim_{t \to \infty} \frac{1}{Nt} \left\langle \int_0^t d\tau r(C(\tau)) \right\rangle_s \quad (46)$$

# IV. DYNAMICAL TRANSITIONS IN MODELS OF GLASS-FORMERS

## A. Existence of a transition in KCMs: variational bounds

It is clear from their equilibrium distributions $P_{eq}(C)$ that KCMs have no phase transitions at any finite temperature. That is, their thermodynamic free energies are analytic functions of temperature (or chemical potential). However, we now show that in the limit of large time $t$ and large system size $N$, the dynamical free energy density $N^{-1}\psi_K(s)$ has a singularity at $s = 0$. To be precise, the dynamical free energy has a discontinuous first derivative with respect to $s$, so we interpret this singularity as a dynamical analog of a first-order phase transition.

The proof of such a transition is based on the escape rates $r(C)$ from the configurations of the model. We establish two bounds on $\psi(s)$. Firstly, the number of jumps $K$ is non-negative, so Eq. (14) implies that $Z_K(s,t)$ is a non-increasing function of $s$. Thus $\psi_K(s)$ is also non-increasing. Further, $\psi_K(0) = 0$, so we have

$$\psi_K(s) \le 0, \quad s \ge 0 \quad (47)$$

Secondly, we can use the variational result (27) with $V(C_1) = 1$ for just one configuration $C_1$, and $V(C) = 0$ otherwise to establish

$$\psi(s) \ge -\min_C[r(C)] \quad (48)$$

for all $s$. For our purposes, the most important property of the kinetically constrained models defined above is that they have

$$\lim_{N \to \infty} N^{-1} \min_C r(C) = 0 \quad (49)$$

This can be established by explicit construction. In the FA and East models, we simply consider a configuration containing exactly one excitation, which has escape rate $2dc$ in the FA case and $dc$ in the East model, where $d$ is the spatial dimension (in the bosonic variants, these rates are $2de^{-\beta}$ and $de^{-\beta}$). In the (2)-TLG, all of the particles in the model can be arranged in a single compact cluster, in which all but six of the particles are unable to move: this configuration has $r(C) = 6$. For the KA model, a similar construction leads to configurations with $r(C) = 4$. Thus, combining (47-49), we have

$$\lim_{N \to \infty} \frac{1}{N} \psi_K(s) = 0, \quad s \ge 0 \quad (50)$$

Recalling that $\langle K \rangle = t(d/ds)\psi_K(s)$, We define the mean activity per site per unit time as

$$K(s) = \lim_{N \to \infty} \frac{1}{N} \frac{d\psi_K(s)}{ds} \quad (51)$$

and we can see that

$$K(s) = 0, \quad s > 0 \quad (52)$$

Further, from Eq. (B9), we have $K(0) = t\langle r \rangle = t \sum_C P_{eq}(C)r(C)$. Since the distributions $P_{eq}(C)$ have simple forms in kinetically constrained models, this quantity can be calculated explicitly: the limit $K_{eq} = \lim_{N \to \infty} N^{-1} \sum_C P_{eq}(C)r(C)$ is finite and positive for all the models that we consider. Finally, it follows from (14) that $K(s)$ is non-increasing, so that

$$K(s) \ge K_{eq}, \quad s \le 0 \quad (53)$$

with $K_{eq}$ finite. Eqs. (52) and (53) establish the discontinuity of $K(s)$ at $s = 0$: in the limit of large system size, the dynamical free energy has a discontinuous first derivative which we refer to as a first-order dynamical phase transition. We have established the existence of such a transition in the FA, East and (2)-TLG models, in all dimensions and for all finite temperatures [and for all finite densities $\rho$ in the (2)-TLG]. That is, the simple phase diagram shown in fig 2 is generic to all of these models.
<figure><img src="image_3.png" /><figcaption>FIG. 2: (Left): Generic 'dynamic phase diagram' for spin-facilitated KCMs such as the FA and East models. There is a dynamical phase coexistence boundary at $s = 0$, for all finite temperatures. The boundary ends in a dynamical critical point at $s = T = 0$. For the (2)-TLG and KA models, the picture is identical if the temperature $T$ is replaced by the fraction of vacant sites, $1 - \rho$. (Center): Variational estimates for the activity per site $\mathcal{K}(s)$, in the bosonic FA model, at $c = 0.25$. For $s > 0$, the result $\mathcal{K}(s) = 0$ is exact in all dimensions, from Eq. (52). For $s < 0$, the dashed line shows the lower bound obtained from (53), while the solid line is the variational estimate $\mathcal{K}_{\text{var}}(s) = -\frac{1}{N}(\text{d}/\text{ds})\psi_{\text{var}}(s)$, obtained from (58). As discussed in the text, the solid line gives the exact result for the mean-field variant of the FA model. (Right) Again, we show the exact result $\psi_K(s) = 0$ for $s > 0$, together with the variational lower bounds (58) (solid line, exact for the mean-field variant) and (48) (dashed line).</figcaption></figure>

## B. Variational free energy for the excitation density $\rho_K(s)$

The analysis given above establishes some minimal conditions that are sufficient for the existence of a first-order transition. For a more quantitative analysis, it is useful to use a specific variational distribution in (27). We consider a general bosonic KCM with single spin-flip dynamics, and we define a distribution of the excitation numbers $n_i$ that is independent of the site $i$, and parameterized by a mean density $\rho$:

$$V_{\rho}(\{n_i\}) = \prod_i \sqrt{\frac{\rho^{n_i} e^{-\rho}}{n_i!}} \qquad (54)$$

From (27), we therefore have $\psi(s) \ge -N \min_{\rho} \mathcal{F}_K(\rho, s)$ with

$$\mathcal{F}_K(\rho, s) \equiv N^{-1} \frac{\langle V_{\rho} | \tilde{W}_K | V_{\rho} \rangle}{\langle V_{\rho} | V_{\rho} \rangle} \qquad (55)$$

The value of $\rho$ which minimises $\mathcal{F}_K(\rho, s)$ is denoted by $\rho_{\text{var}}(s)$. It represents a variational estimate for the order parameter $\rho_K(s)$: if the variational bound (27) is saturated then $|V\rangle$ is an eigenvector of the symmetrised operator $\tilde{W}_K$, and it follows from (A14) that $\rho_K(s) = \rho_{\text{var}}(s)/(1 - e^{-\rho N})$.

For the bosonic FA model, it is straightforward to calculate $\mathcal{F}_K(\rho, s)$. The only subtlety is that we must explicitly exclude the state with no excitations from the inner products, as discussed in Section III B. In the Doi-Peliti formalism, our choice for $V(\{n_i\})$ renders this calculation very simple: in terms of the symmetrised operator $\tilde{W}_K^{(\text{FA})}$, we have

$$\mathcal{F}_K(\rho, s) = N^{-1} \frac{\langle 0 | e^{-\sqrt{\rho} \sum_i a_i} \tilde{W}_K^{(\text{FA})} e^{-\sqrt{\rho} \sum_i a_i^\dagger} | 0 \rangle e^{-\rho N}}{1 - e^{-\rho N}} \qquad (56)$$

Hence,

$$\mathcal{F}_K(\rho, s) = 2d \frac{c + \rho - 2e^{-s} \sqrt{c\rho}}{1 - e^{-\rho N}} \rho. \qquad (57)$$

Minimising over $\rho$, we find $\lim_{N \to \infty} N^{-1} \psi(s) \ge \psi_{\text{var}}(s)$ with

$$\psi_{\text{var}}(s) = \frac{2d}{3} \rho_{\text{var}}(s) [\rho_{\text{var}}(s) - c] \qquad (58)$$

and

$$\rho_{\text{var}}(s) = \begin{cases} 0, & s > 0 \\ (c/8)(9e^{-2s} - 4 + 3e^{-s}\sqrt{9e^{-2s} - 8}), & s \le 0 \end{cases} \qquad (59)$$
<figure><img src="image_4.png" /><figcaption>FIG. 3: (Top left) Variational free-energies for the FA model (constrained dynamics, solid line), and the A-model (unconstrained dynamics, dashed line) at $s = 0$. Both models have the same thermodynamic free energy. However, the dynamical function $\mathcal{F}(\rho, s)$ reveals that the FA model has two dynamical phases while the A-model has only an active phase. (Top right) The variational estimate $\rho_{\text{var}}(s)$, for the FA model (solid line) and the A-model (dashed). For the mean-field FA model, $\rho_{\text{var}}(s)$ coincides with $\rho_K(s)$ in the limit of large system size $N$; for the A-model, $\rho_K(s) = \rho_{\text{var}}(s)$ always. (Bottom) Dependence of the variational free energy on the field $s$, in the FA model. At the phase coexistence condition, $s = 0$, the free energy has degenerate minima. For finite $s$, either the inactive or active phase is preferred.</figcaption></figure>

(Within this approach, we obtain $\rho_{\text{var}}(0) = c$ by minimising $\mathcal{F}(\rho, s)$ at fixed system size $N$, and then taking $N \to \infty$.) The bound on $\psi_K(s)$ and the corresponding estimate of $\mathcal{K}(s)$ are shown in Fig. 2. The variational estimate for $\rho_K(s)$ and the variational free energy $\mathcal{F}_K(\rho, s)$ are shown in Fig. 3.

So far, we have used Eq. (27) to obtain variational estimates for $\psi_K(s)$ and $\rho_K(s)$ for the FA model in finite dimension. For the mean-field variant of the FA model, it can be shown that these variational estimates are exact, in the limit of large system size $N$. (The factor $2d$ that appears in $\mathcal{F}(\rho, s)$ is simply an arbitrary rescaling of time in the mean-field model. Our definition of the mean-field model requires that we set $2d = 1$.) That is, the difference between the variational ansatz of (54) and the dominant eigenvector of $\tilde{W}_K$ vanishes at large $N$. Mean-field models are discussed in more detail in Section V below.

It is useful to compare these results for the FA model with the bosonic variant of the A-model, for which it can be easily verified that the large deviation function $\psi_K(s)$ coincides with the variational bound $\psi_{\text{var}}(s)$, even for finite system size $N$. In that case, we have

$$ \mathcal{F}(\rho, s) = c + \rho - 2e^{-s}\sqrt{\rho c} \quad \psi_K(s) = c(e^{-2s} - 1), \quad \rho_K(s) = ce^{-2s} \quad (60) $$

Thus, while constrained FA model and the unconstrained A-model possess the same equilibrium distribution $P_{\text{eq}}(C)$, and hence the same static free energies, their dynamical free energies show dramatic differences. For large systems, the FA model exhibits a dynamical phase transition, while the A-model does not. See Fig. 3.

The presence of the dynamical first-order transition in the FA model is intimately connected to the two minima in $\mathcal{F}_K(\rho, s)$. As shown in Fig. 3, these two minima represent an active phase, with $\rho \simeq c$ and an inactive one with $\rho \simeq 0$. The global minimum of $\mathcal{F}(\rho, s)$ lies in the active phase for $s < 0$, while it lies in the inactive phase for $s > 0$. For $s = 0$, one must consider carefully the limit of large system size $N$: we have $\mathcal{F}_K(c, 0) = 0$ while the inactive minimum occurs at $\rho_{\text{var}} = \mathcal{O}(N^{-2})$, where the value of the variational bound $\psi_{\text{var}}$ is positive. Thus, the global minimum of $\mathcal{F}_K(\rho, s = 0)$ occurs at the active state density. However, any $s > 0$ is sufficient to drive the system into the inactive phase. The effect arises because of two non-commuting limits: when minimising $\mathcal{F}(\rho, s)$, taking the limit $s \to 0$ before the limit of large $N$ results in active behaviour; on the other hand, taking the limit of large $N$ followed by a limit $s \to 0^+$ leads to the inactive phase.

We note that the dynamical phase transition in the FA model requires a limit of large system size ($N \to \infty$) as well as a limit of long trajectories ($t \to \infty$). To keep our methods well-defined, we excluded the configuration with no excitations from the initial conditions (recall Section III B). We emphasise that we have proved the existence of
a dynamical phase transition in an irreducible model, with no absorbing states (this can be compared, for example, with phase transitions in the directed percolation universality class [43]).

## C. Numerical results

### 1. Cloning method

We now present some numerical computations of the dynamical free energy $\psi_K(s)$ in KCMs. From (14), this can be obtained from the large $t$ limit of the equilibrium average $\langle e^{-sK} \rangle$. However, direct calculation of this average requires a computational effort that scales exponentially with $t$: the average is dominated by rare histories lying in the tails of the distribution of $K$. In dynamical systems [44] and in discrete time Markov processes [19], this problem can be avoided by using a cloning method similar to that used in quantum mechanical Diffusion Monte Carlo algorithms [45]. This method was generalised to continuous time Markov processes by Tailleur and Lecomte [20]. Here, we briefly summarize the algorithm for obtaining dynamical free energies.

The function $\psi_A(s)$ is obtained as the largest eigenvalue of the operator $W_A$. However, this operator does not conserve probability [that is, $W_A$ sets the time dependence of $\hat{P}_A(C, s, t)$, but the 'total probability' $\sum_C \hat{P}_A(C, s, t)$ is not a constant of the motion, except at $s = 0$]. To interpret this non-conservation, we define a new stochastic process (a 'modified dynamics') with rates $W_s(C' \to C)$, chosen so that we can decompose (16) as

$$ \partial_t \hat{P}_A(C, s, t) = \sum_{C'} [W_s(C \to C') - r_s(C)\delta_{C,C}] \hat{P}_A(C', s, t) + \delta r_s(C) \hat{P}_A(C, s, t) \quad (61) $$

with $r_s(C) = \sum_{C'} W_s(C \to C')$ and $\delta r_s(C) = r_s(C) - r(C)$. This decomposition is discussed in appendix B, and the rates $W_s(C \to C')$ are given in (B1).

For the purposes of the cloning algorithm, we note that the first term in (61) conserves probability (in the sense given above), while the second term represents the creation or destruction of copies (clones) of the system. That is, starting from a large number of copies of the system, we let each copy of the system evolve with the modified dynamics (rates $W_s$). In addition, the copies are subject to a creation/destruction process with a configuration-dependent rate $\delta r_s(C)$. In this way, the number $n_{cl}(C, t)$ of copies of the system in configuration $C$ at time $t$ has the same time evolution as $\hat{P}_A(C, s, t)$ in (61). To avoid the ensuing exponential increase or decrease of the total number of copies [which behaves as $e^{t\psi_A(s)}$], one compensates the clone creation/destruction rates of (61) with configuration-independent creation/destruction rates. The rates are adapted as the simulation proceeds, in order to keep a constant clone population [19, 20]. These adaptively determined rates can then be used to obtain the dynamical free energy $\psi_A(s)$.

### 2. Results

Using the cloning method, we investigated two classes of KCMs. In Fig. 4 we consider spin-facilitated models: the FA model in one dimension and the East model in three dimensions. We evaluated the free energy density $\frac{1}{N}\psi_K(s)$ for various system sizes. Its behavior as a function of $N$ drastically depends on the sign of $s$, as is also the case for the order parameter $\rho_K(s)$. Negative values of $s$ correspond to active histories, with non-zero mean density of particles, while for positive values of $s$, the mean number of particle in the system remains finite, leading to a zero density and activity in the infinite system size limit.

In Fig. 5, we consider two models with particle conservation: the KA and (2)-TLG models, both in two dimensions. Remarkably, the picture is very similar to the previous one, the (conserved) density being replaced with the order parameter $r_K(s)$.

In the active phase ($s < 0$), the order parameters $\rho_K(s)$ and $r_K(s)$ converge rapidly as the system size $N$ is increased, for all the KCMs that we considered. On the other hand, in the inactive phase ($s > 0$) the order parameters decrease with system size as $N^{-1}$. Comparing Fig. 4 and Fig. 2 confirms the analysis based on variational bounds on $\psi_K(s)$: in the limit of large system size, KCMs exhibit dynamical first-order transitions at $s = 0$. For models on finite lattices, the equilibrium ($s = 0$) dynamics are representative of the active phase, and the system crosses over to the inactive phase at a value of $s$ that scales as $N^{-1}$ for large $N$.
<figure><img src="image_5.png" /><figcaption>FIG. 4: Finite-size scaling in the s-ensemble associated to K, in the vicinity of the singular point $s_c = 0$ for KCMs with non-conserved number of particles: Fredrickson-Andersen model in dimension $d = 1$ (top) and East model in dimension $d = 3$ (bottom). The temperature is $T = 1/\beta = 0.91$, and linear system sizes $L$ are given on the graphs. The finite-size scaling illustrates the first order dynamical phase transition in $s_c = 0$. (Left) large deviation function $1/L^d \psi_K(s)$. (Right) density of excited sites $\rho_K(s)$.</figcaption></figure>

### D. Criticality at zero temperature and dynamical phase transition

We emphasize that while a zero temperature dynamical critical point is common to many KCMs [36], this is not a sufficient condition for dynamical phase coexistence. Rather, the relevant feature is the presence of states with subextensive escape rates, as discussed in Section IV A. In this section, we consider the AA model. The FA and AA models both have zero-temperature dynamical critical points, with the same scaling exponents and closely related correlation functions [37]. However, all states in the AA model have extensive escape rates, so we do not expect any transition at $s = 0$. In the following, we show that this is indeed the case, by discussing the AA model both within a mean field approximation and using exact results in one dimension.

#### 1. 'Mean-field' variational bound

We consider the bosonic AA model in dimension $d$. Following Section IV B, we calculate the variational Landau free energy using the Doi-Peliti representation, obtaining

$$ \mathcal{F}_K(\rho, s) = 2d[2D\rho(1 - e^{-s}) + k' + k\rho^2 - 2e^{-s}\rho\sqrt{kk'}] \quad (62) $$

and we identify $\psi_{\text{var}}(s) = -\min_{\rho} F_K[\rho, s]$ as a lower bound on $N\psi_K(s)$.

The variational estimate for the steady-state density and the variational bound are

$$ \rho_{\text{var}}(s) = e^{-s}\sqrt{\frac{k'}{k}} + (e^{-s} - 1)\frac{D}{k} \quad (63) $$

$$ \psi_{\text{var}}(s) = 2d(k\rho_{\text{var}}^2(s) - k') \quad (64) $$
<figure><img src="image_6.png" /><figcaption>FIG. 5: Finite-size scaling in the *s*-ensemble associated to *K*, in the vicinity of the singular point $s_c = 0$ for KCMs with conservation of particle number: Kob-Andersen model (**top**) and Triangular Lattice Gas model (**bottom**) Particle density is $\rho = 0.5$, and linear system sizes *L* are given on the graphs. The finite-size scaling illustrates the first order dynamical phase transition in $s_c = 0$. (**Left**) large deviation function $1/L^d \psi_K(s)$. (**Right**) order parameter $1/L^d \langle r \rangle_s$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6: Dynamical free energy density $1/N \psi_K(s)$ for the AA model in dimension one. It presents no singularity at $s = 0$. Although the AA model can be mapped to the FA model and displays the same critical properties at zero temperature, it is not subject to dynamical phase coexistence.</figcaption></figure>

The variational bound is indeed analytic for all $s$, consistent with our intuition that the AA model has no dynamical phase transition. Again, these variational estimates are exact for the mean-field AA model in the limit of large system size, if we set $2d = 1$.

## 2. AA model in one dimension

In addition to the mean-field case, we can also obtain the large deviations of the AA model in one dimension, through a mapping to a free fermion system. The evolution operator associated to *K* for AA model can be written
in a spin-$\frac{1}{2}$ representation (recall Eq. (35) and see, for example, Ref. [46]):

$$ \mathbb{W}_K = \sum_i \left\{ z[k\sigma_i^- \sigma_{i+1}^- + k'\sigma_i^+ \sigma_{i+1}^+ + D(\sigma_i^- \sigma_{i+1}^+ + \sigma_i^+ \sigma_{i+1}^-)] \right. \\ \left. - k\hat{n}_i\hat{n}_{i+1} - k'(1-\hat{n}_i)(1-\hat{n}_{i+1}) - D\hat{n}_i(1-\hat{n}_{i+1}) - D\hat{n}_{i+1}(1-\hat{n}_i) \right\} \quad (65) $$

where $z = e^{-s}$, $\sigma_i^\pm = \frac{1}{2}(\sigma_i^x \pm i\sigma_i^y)$, $\hat{n}_i = \frac{1+\sigma_i^z}{2}$ and $\sigma_i^x, \sigma_i^y, \sigma_i^z$ are the usual Pauli matrices. In the spin language, the presence (or absence) of a particle at site $i$ is coded by an up (or down) spin. We use the detailed balance property to symmetrise this operator followed by a Jordan-Wigner transformation [46, 47]

$$ \sigma_i^+ = c_i \exp \left( i\pi \sum_{j=1}^{i-1} c_j^\dagger c_j \right) \quad \sigma_i^- = \exp \left( i\pi \sum_{j=1}^{i-1} c_j^\dagger c_j \right) c_i^\dagger \quad (66) $$

which allows us to represent the spin operators in terms of fermionic creation/annihilation operators $c_j^\dagger$ and $c_j$. For values of the parameters verifying $k + k' = 2D$ [48], this puts $\tilde{W}_K$ into a quadratic form:

$$ \tilde{W}_K = -\sum_q \left[ (k-k' - z(k+k')) \cos q \right] c_q^\dagger c_q - iz\sqrt{kk'} c_q^\dagger c_{-q}^\dagger \sin q + iz\sqrt{kk'} c_{-q} c_q \sin q - k'N \quad (67) $$

where we introduced Fourier-transformed operators $c_q = \sum_j c_j e^{iqj}$ and $c_q^\dagger = \sum_j c_j^\dagger e^{-iqj}$.

We now introduce new fermionic operators $\beta_q = \cos\theta_q c_q - i\sin\theta_q c_{-q}^\dagger$, $\beta_q^\dagger = \cos\theta_q c_q^\dagger + i\sin\theta_q c_{-q}$. Taking $\pi/4 < \theta_q < \pi/4$, we write $\sin 2\theta_q = (2z\sqrt{kk'} \sin q/\Omega_q)$, with

$$ \Omega_q = \sqrt{4kk'(z^2-1) + [k+k'-z(k-k')\cos q]^2} \quad (68) $$

Then, the dynamical free energy is the largest eigenvalue of the operator

$$ \tilde{W}_K = \frac{1}{2} \sum_q [\Omega_q(1 - 2\beta_q^\dagger\beta_q) - (k + k')] \quad (69) $$

Finally, for large $N$, we convert the sum over $q$ to an integral, arriving at

$$ \psi_K(s) = \frac{N}{2} \left[ -(k + k') + \int \frac{dq}{2\pi} \Omega_q \right] \quad (70) $$

which depends on $s$ through the dependence of $\Omega_q$ on $z = e^{-s}$.

This exact result for the AA model in $d=1$ shows that the large deviation function $\psi_K(s)$ is analytic at $s=0$ (see also Fig. 6), as opposed to the one of the FA model. Despite the presence of a dynamical critical point, the AA model has no configurations with subextensive escape rate $r(C)$, and does not exhibit dynamical phase coexistence (in the vicinity of $s=0$).

# V. PROPERTIES OF TRAJECTORIES IN THE s-ENSEMBLE

We have proven the existence of a first-order dynamical phase transition in KCMs, and compared the behaviour of these models with the A and AA models. The effect of the field $s$ is to generate an ensemble of histories, biased towards small or large activity. In order to gain insight into this transition, we now discuss the histories that dominate the $s$-ensemble when $s$ is finite.

## A. Effect of temporal boundary conditions in the s-ensemble

### 1. General considerations

In steady states, the (unbiased) ensemble of histories is invariant under translation in time. Suppose that $b = b(C)$ is a configuration-dependent observable, and $B = \int_0^t dt' b(C(t'))$. Then, for trajectories of length $t$ the expectation value of the observable $b$ at time $\tau$ is

$$ \langle b(\tau) \rangle = \frac{1}{t} \langle B \rangle, \quad (71) $$
independent of the time $\tau$.

However, introducing a field $s$ biases the ensemble of histories, and, in general, time translation invariance is broken. This effect is a dynamical analog of boundary effects in classical thermodynamics: if a system is finite, the behaviour near its boundaries is different from that of the bulk. In the $s$-ensemble, we consider trajectories $C(\tau)$: the boundaries of the trajectory are $\tau = 0$ and $\tau = t$, while the analogy of the ‘bulk’ is $0 \ll \tau \ll t$. In the limit of large time, extensive quantities are dominated by the bulk: we have

$$ \langle b(\tau) \rangle_s = \frac{1}{t} \langle B \rangle_s, \qquad (72) $$

for $0 \ll \tau \ll t$. However, in general we have $\langle b(\tau) \rangle_s \neq \langle b(t) \rangle_s \neq \langle b(0) \rangle_s$. (In section V A 2 we illustrate these differences by calculating $\langle b(\tau) \rangle$ in the A-model.)

More generally, it is possible to express the average at the final time, $\langle b(t) \rangle_s$ and the time average $\frac{1}{t}\langle B \rangle_s$ by considering the eigenvalues and the eigenvectors of the operator $W_A$ (see appendix A). Using this approach, one can perform a perturbation theory around $s = 0$. In particular, when detailed balance is verified and $s$ is conjugate to an observable of type $B$, the bulk and boundary averages differ at first order in $s$: for large times

$$ \langle b(t) \rangle_s = \langle b \rangle + sb^{(1)} + \dots \qquad (73) $$

$$ \langle B \rangle_s = t [\langle b \rangle + 2sb^{(1)} + \dots], \qquad (74) $$

where an explicit expression for $b^{(1)}$ is given in Eq. (A13).

## 2. Effects of temporal boundaries in the A-model

We now illustrate the effect of temporal boundaries in the $s$-ensemble, using the (bosonic) A-model. We define the average particle density in this ensemble at a time $\tau$: that is,

$$ \rho(s; \tau, t) = \frac{1}{N} \langle n_{\text{tot}}(\tau) \rangle_s, \qquad (75) $$

with $0 \le \tau \le t$ and $n_{\text{tot}} = \sum_i n_i$. To completely specify the problem, we must set the initial conditions in (7): we take a Poisson distribution with mean density $c_0$:

$$ p_0(C_0) = \prod_i \frac{c_0^{n_i} e^{-c_0}}{n_i!} \qquad (76) $$

To proceed, we write $N_\tau[\text{hist}] = \int_0^{\min(\tau,t)} dt' n_{\text{tot}}(t')$, and we define $P(n_{\text{tot}}, N_\tau, K, t)$ to be the probability that the system contains $n_{\text{tot}}$ excitations at time $t$, having made $K$ changes of configuration, and with the observable $N_\tau[\text{hist}]$ taking a value $N_\tau$. Then, we define the generating function

$$ \hat{P}_{n_{\text{tot}}} \equiv \hat{P}(n_{\text{tot}}, h, s, t) = \sum_K \int dN_\tau e^{-hN_\tau - sK} \hat{P}(n, N_\tau, K, t) \qquad (77) $$

so that

$$ \rho(s; \tau, t) = -\frac{1}{N} \frac{\partial}{\partial \tau} \frac{\partial}{\partial h} \ln \sum_{n_{\text{tot}}} \hat{P}(n_{\text{tot}}, h, s, t) \Big|_{h=0} \qquad (78) $$

Deriving an equation of motion for $\hat{P}_{n_{\text{tot}}}$ is a straightforward generalisation of the derivation of (16): the result is

$$ \partial_t \hat{P}_{n_{\text{tot}}} = e^{-s} [cN \hat{P}_{n_{\text{tot}}-1} + (n+1) \hat{P}_{n_{\text{tot}}+1}] - [cN + n + hn\Theta(\tau - t)] \hat{P}_{n_{\text{tot}}} \qquad (79) $$

where $\Theta(t)$ is the usual Heaviside step function. With the initial condition of (76), this equation of motion is solved by a Poisson distribution with a time-dependent normalisation factor: we take $P(n_{\text{tot}}, h, s, t) = \exp[\Psi(t) - \rho_0(t)]\rho_0(t)^{n_{\text{tot}}}/(n_{\text{tot}}!)$. Then, the mean, $\rho_0(t)$, and normalisation factor, $\Psi(t)$, obey

$$ \dot{\rho}_0(t) = cN e^{-s} - [1 + h\Theta(\tau - t)]\rho_0(t) \qquad (80) $$

$$ \dot{\Psi}(t) = \dot{\rho}_0(t) + e^{-s}\rho_0(t) - cN \qquad (81) $$
<figure><img src="image_8.png" /><figcaption>FIG. 7: Mean density as a function of time in the A-model, for the s-ensemble of histories of length $t = 30$. We take $c_0 = 0.1$, $c = 0.2$, $s = -0.5$.</figcaption></figure>

with initial conditions $\rho_0(0) = c_0N$, $\Psi(0) = 0$. We identify $e^{\Psi(t)} = \sum_{n_{tot}} P(n_{tot}, h, s, t)$, so we solve for $\Psi(t)$ and use (78) to obtain

$$ \rho(s; \tau, t) = cz^2 + e^{-t}(1-z)(c_0 - cz) + e^{-\tau}z(c_0 - cz) + e^{\tau-t}cz(1-z) \quad (82) $$

where we defined $z = e^{-s}$ for ease of writing.

The $\tau$-dependent density $\rho(s; \tau, t)$ exhibits four different regimes (Fig. 7).

* Short trajectories. For $\tau < t \ll 1$, the system has a density $\rho(s; \tau, t) \simeq c_0$ close to the density at time 0.

* Long trajectories, stationary (bulk) regime. For $1 \ll \tau \ll t$, the system adopts a density $\rho(s; \tau, t) \simeq ce^{-2s}$, independent of $\tau$. For long trajectories, this average value coincides with the time averaged density $t^{-1} \int_0^t d\tau \rho(s; \tau, t)$.

* Long trajectories, initial transient regime. For early times $\tau \ll 1 \ll t$, the density depends on the value of $s$. This dependence persists even for $\tau = 0$: that is, the trajectories that dominate the s-ensemble have non-typical initial conditions as well as non-typical bulk properties. To be precise, $\rho(s; 0, t) = c_0e^{-s}$: the influence of the initial condition decays into the bulk as $\rho(s; \tau, t) \simeq (c_0e^{-s} - ce^{-2s})e^{-\tau} + ce^{-2s}$.

* Long trajectories, final transient regime: for $\tau \to t$ the density at the final time $t$ is $\rho(s; t, t) = ce^{-s}$. Moving away from this boundary, the density decays into the bulk as $\rho(s; \tau, t) \simeq (ce^{-s} - ce^{-2s})e^{\tau-t} + ce^{-2s}$.

We note that if the initial density $c_0$ is equal to the equilibrium density $c$, then the ensemble at $s = 0$ has time-reversal symmetry. Since the observable $K$ respects this symmetry, the s-ensemble is also time-reversal symmetric $\rho(s; \tau, t) = \rho(s; t - \tau, t)$, for $c = c_0$. In this case, the initial and final transient regimes are related by this symmetry.

We have used the A-model to calculate the time dependence of $\rho(\tau)$ exactly. However, we emphasise that the four regimes identified here are very general. When $t$ is large, histories in the s-ensemble are characterised by an extended intermediate (bulk) regime, with initial and final transient regimes that decay exponentially into the bulk.

## B. Landau-like theory for fluctuations within the s-ensemble

In this section, we study large deviations of observables within the s-ensemble. For example, for an s-ensemble parameterized by the observable $K$, we consider the probability of observing a history with a particular value of an observable $B$ unrelated to $K$. In particular, we connect the large deviations of the average excitation density $\rho$ to the variational free energy $F_K(\rho, s)$, defined in (55).

### 1. Variational calculation of $\psi_K(s)$ in a general mean-field model

We consider systems for which we can write the master-like equation (16) in terms of a single co-ordinate $n_{tot}$. In the mean-field FA model, this co-ordinate is the total number of excitations, but it might also represent (for example),
the total magnetisation of a mean-field Ising model [8]. To be precise, we assume the master-like operator $W_K$ has matrix elements

$$ (W_K)_{n,n'} = e^{-s}W_{n-1}^{+} \delta_{n',n-1} + e^{-s}W_{n+1}^{-} \delta_{n',n+1} - [W_{n}^{+} + W_{n}^{-}] \delta_{n',n}. \quad (83) $$

where we have abbreviated $n_{tot}$ to $n$, for compactness, and $W_n^{\pm}$ are the rates for transitions from the state $n$ to state $n \pm 1$. For an example, consider the (bosonic) mean-field variant of the FA model, for which $W_n^+ = W(n \to n+1) = cn$ and $W_n^- = W(n \to n-1) = (n-1)n/N$, as defined in Section III C. We have assumed for convenience that all processes in the system change the co-ordinate $n_{tot}$ by one: the generalisation to other cases (such as the mean-field AA model) is straightforward.

Following Section II B 6, we now symmetrise the operator $W_K$, so that the dynamical free energy $\psi_K(s)$ is given by the largest eigenvalue of the operator

$$ (\tilde{W}_K)_{n,n'} = (W_n^+ W_{n+1}^-)^{1/2} e^{-s} \delta_{n',n+1} + (W_{n'}^+ W_{n'+1}^-)^{1/2} e^{-s} \delta_{n,n'+1} - (W_n^+ + W_n^-) \delta_{n,n'} \quad (84) $$

For large systems, $(N \to \infty)$ the eigenvector associated with the largest eigenvalue takes the form $V_n = e^{-Nf(\rho)/2}$ with $\rho = n_{tot}/N$, and the function $f(\rho)$ has a unique global minimum. Then, Eq. (27) states that

$$ \psi_K(s) = \max_{f(\rho)} \frac{\sum_{n_{tot}} \{e^{-s}[W_{n_{tot}}^{+} W_{n_{tot}}^{-}]^{1/2} (e^{-f'(\rho)} + e^{f'(\rho)}) - W_{n_{tot}}^{+} - W_{n_{tot}}^{-}\} e^{-Nf(\rho)}}{\sum_{n_{tot}} e^{-Nf(\rho)}} \quad (85) $$

For any trial function $f(\rho)$, the sums over $n_{tot}$ in (85) are dominated by the occupation numbers $n_{tot}$ such that $f(\rho)$ is minimum (which implies in particular $f'(\rho) = 0$). Thus, the direct dependence on $f$ vanishes: we are left with a maximisation over the position of the minimum in $f(\rho)$. Since the form of $f(\rho)$ is irrelevant, the eigenvector $V_n$ can be written in the form given in Eq. (54). Using this choice, we arrive at

$$ \psi_K(s) = - \min_{\rho} \mathcal{F}_K(\rho, s) \quad (86) $$

where the variational free energy $\mathcal{F}_K(\rho, s)$ was originally defined in (55). For these mean-field models, it takes the form

$$ \mathcal{F}_K(\rho, s) = \frac{1}{N} \{-2e^{-s}[W_{N\rho}^{+} W_{N\rho}^{-}]^{1/2} + W_{N\rho}^{+} + W_{N\rho}^{-}\} \quad (87) $$

As discussed in Section IV B, $\mathcal{F}_K(\rho, s)$ gives a bound on $\psi_K(s)$ for all systems. However, for systems with mean-field geometry, we showed that the form of the trial distribution is irrelevant in the limit of large system size. Thus, we write Eq. (86) with an equality, and not as a bound. We now discuss the physical interpretation of this result.

## 2. Physical interpretation of the variational free energy

Consider an $s$-ensemble in which trajectories are weighted by the usual factor $e^{-sK[\text{hist}]}$, but with the further restriction that the time-integrated density be fixed. That is, we write the (unnormalized) probability, in the $s$-ensemble, to measure a time-averaged density $\rho$,

$$ \left\langle e^{-sK} \delta \left[ \rho - \frac{1}{Nt} \int_{0}^{t} d\tau \, n(\tau) \right] \right\rangle \sim e^{-Nt \mathcal{F}_{K}^{*}(\rho, s)} \quad (88) $$

where the asymptotic behaviour of the left hand side at large $t$ defines the function $\mathcal{F}_K^*(\rho, s)$. Taking a Laplace transform of (88) with respect to $\rho$, we arrive at

$$ Z_{K,\mathcal{N}}(s,h) = \left\langle \exp\left(-sK - h\mathcal{N}\right) \right\rangle = \int d\rho \, e^{-(\mathcal{F}_K^*(\rho,s)+hNt\rho)} \quad (89) $$

where we write $\mathcal{N} = \int_0^\tau d\tau \, n_{tot}(\tau)$, noting the similarities with the generating function of (78).

Now, by analogy with (14), we identify $Z_{K,\mathcal{N}}(s,h)$ as the partition function for an $(s,h)$-ensemble, in which histories are biased both by their activity $K$ and their time-integrated number of excitations $\mathcal{N}$. Repeating the analysis of section II B 3, we observe that

$$ \psi_{K,\mathcal{N}}(s,h) = \lim_{t \to \infty} \frac{1}{t} \ln Z_{K,\mathcal{N}}(s,h) \quad (90) $$
is the largest eigenvalue of an operator $W_{K,\mathcal{N}}$, whose elements are

$$ (W_{K,\mathcal{N}})_{n',n} = W_{n-1}^{+} e^{-s} \delta_{n',n-1} + W_{n+1}^{-} e^{-s} \delta_{n',n+1} - [W_{n}^{+} + W_{n}^{-} + hn] \delta_{n',n} \quad (91) $$

The largest eigenvalue of this operator can be obtained by symmetrising and repeating the variational analysis of the previous section. The result is

$$ \psi_{K,\mathcal{N}}(s,h) = -\min_{\rho}(\mathcal{F}_K(\rho,s) + h\rho) \quad (92) $$

which applies in the limit of large system size $N$ [since in that case, the maximisation over the function $f(\rho)$ can be replaced by a maximisation over the density $\rho$]. However, performing a saddle point analysis directly on (88) reveals (for large times $t$ and finite system size $N$),

$$ \psi_{K,\mathcal{N}}(s,h) = -\min_{\rho}(\mathcal{F}_{K}^{*}(\rho,s) + h\rho) \quad (93) $$

Thus, in the limit of large system size $N$, the Legendre transforms of $\mathcal{F}_K(\rho, s)$ and $\mathcal{F}_K^*(\rho, s)$ are equal. It follows that the large deviation function $\mathcal{F}_K^*(\rho, s)$ coincides with the variational free energy $\mathcal{F}_K(\rho, s)$ as long as the inverse Legendre transform can be performed. However, in Section IV B, we showed that in KCMs, $\mathcal{F}(\rho, s)$ typically has two minima, separated by a range of densities in which it is ‘non-convex’: $\partial_{\rho}^2 \mathcal{F}(\rho, s) > 0$. In this case, the inverse Legendre transformation cannot be performed. In fact, the non-convexity of $\mathcal{F}(\rho, s)$ arises because histories with some values of $\rho$ are unstable in the $s$-ensemble, as we now show.

### 3. Non-convex free energies: phase separation in time

In the thermodynamics of finite-dimensional systems, one typically has $s''(e) \le 0$ [21]. Loosely, this property arises because any energy density $e$ can be achieved by separating a system into two regions, separated by an interface whose energy cost scales subextensively with the size of the system. Thus, the total energy density is $e = (1-x)e_1 + xe_2 + \delta$ where $e_1$ and $e_2$ are the energy densities of the two regions, $x$ is the fraction of the volume of the system taken up by the second region, and $\delta$ is the energy of the interface divided by the total volume, which vanishes in the thermodynamic limit. This leads to the ‘lever rule’ $e = (1-x)e_1 + xe_2$. The total entropy density associated with these configurations is $s(e) = (1-x)s(e_1) + xs(e_2)$, and using the lever rule, it follows that $s''(e) \le 0$. However, in mean-field geometries, interfaces cannot be formed, and this argument cannot be applied.

Interestingly, in the statistics of histories, phase separation is possible even in mean-field systems. We consider the large deviation function $\mathcal{F}_K^*(\rho, s)$, at a density $\rho$ for which $\mathcal{F}_K(\rho, s)$ is non-convex. We will find that the average in (88) is dominated by histories that phase separate in time. To prove this, we use the methods of Donsker-Varadhan theory, described in appendix C. This method allows us to prove that

$$ \mathcal{F}_K^*(\rho, s) = \min_{\substack{|V| \\ \langle V | \hat{\rho} | V \rangle}} \frac{\langle V | \tilde{W}_K | V \rangle}{\langle V | V \rangle} \Big|_{\langle V | \hat{\rho} | V \rangle = \rho} \quad (94) $$

where the minimisation is over distributions $V(C)$ such that $\sum_C V(C)^2 \rho(C) = \rho$. By analogy with the thermodynamic case, we take $V(C) = (1-x)V_{\rho_1}(C) + xV_{\rho_2}(C)$, where $V_{\rho}(C)$ was defined in (54). We then minimise over the densities $\rho_1$ and $\rho_2$, choosing $x = (\rho - \rho_1)/(\rho_2 - \rho_1)$ to ensure that the mean density is $\rho$. Taking $x = 0$, we have a bound $\mathcal{F}_K^*(\rho, s) \le \mathcal{F}_K(\rho, s)$. However, if $\mathcal{F}_K(\rho, s)$ is non-convex (that is, $\partial_{\rho}^2 \mathcal{F}_K(\rho, s) < 0$) we can find a lower bound on $\mathcal{F}_K^*(\rho, s)$ that is smaller than $\mathcal{F}_K(\rho, s)$. For example, in the FA model in finite dimension (and in the limit of large system size $N$), we find that $\mathcal{F}^*(\rho, s)$ is minimised by $\rho_1 = 0$, $\rho_2 = cz^2$ and $x = \rho/\rho_2$, for $0 < \rho < ce^{-2s}$. For $0 < \rho < cz^2$, this variational approximation to $\mathcal{F}_K^*(\rho, s)$ indicates that the system separates into two phases with densities 0 and $cz^2$. We arrive at a bound

$$ \frac{\mathcal{F}_K^*(\rho, s)}{2d} \le \begin{cases} \rho c(1-z^2), & \rho \le cz^2 \\ \rho(c + \rho - 2z\sqrt{c\rho}), & \rho \ge cz^2 \end{cases} \quad (95) $$

from which we note that $\partial_{\rho}^2 \mathcal{F}^*(\rho, s) = 0$ in the two-phase regime: this is the Maxwell construction [21]. The result for the mean-field FA model is obtained by setting $2d = 1$, as in Section IV B: in that case, the bound is saturated (this follows since $\mathcal{F}^*$ is convex and its Legendre transformation is known to be equal to that of $\mathcal{F}$).

In addition to establishing the convexity of $\mathcal{F}_K^*(\rho, s)$, Donsker-Varadhan (DV) theory also provides an interpretation of the distribution $V^*(C)$ that minimises (94). (We normalise to $\sum_C V^*(C)^2 = 1$ for convenience.) At large $t$, we
consider histories in the s-ensemble with fixed average density $\rho$. The DV theorem states that this sub-ensemble is dominated by trajectories for which the fraction of time spent in configuration $C$ is $\mu^*(C) = V^*(C)^2$. If the distribution $\mu^*(C)$ is dominated by configurations $C$ with density $\rho$, we conclude that the histories in this sub-ensemble are homogeneous in time. However, if $\mu^*(C)$ is associated with a bimodal density distribution, it describes histories comprising separate periods of time, some with low excitation density $\rho_1$ and some with high excitation density $\rho_2$.

Finally, we give the interpretation of $\mathcal{F}_K(\rho, s)$. In systems described by the single co-ordinate $n_{tot}$, the DV theorem states that $e^{-Nt\mathcal{F}_K(\rho,s)}$ is the (unnormalised) probability (in the s-ensemble) of a history in which almost all configurations have density $(n_{tot}/N)$ equal to $\rho$. This can be compared with the probability $e^{-Nt\mathcal{F}_K(\rho,s)}$ of a history with a time-averaged density $\rho$. In this sense, $\mathcal{F}_K(\rho, s)$ can be interpreted as a Landau-like free energy for homogeneous trajectories, while $\mathcal{F}_K^*(\rho, s)$ is the large deviation function for the density $\rho$. In finite dimension, $\mathcal{F}_K(\rho, s)$ provides a bound on the large deviation function $\mathcal{F}_K^*(\rho, s)$, based on the assumption that histories are spatially and temporally homogeneous. That is, fluctuations in space and time are neglected. In Section V C, we discuss how these fluctuations can be taken into account.

### 4. Landau-like free energy in other s-ensembles

So far, we have considered the large deviations of the density $\rho$ in an s-ensemble that is defined in terms of the activity $K$. The variational free energy can be simply extended to s-ensembles defined as in (12). Consider again a mean-field model specified by rates $W_n^\pm$ and an observable $A$ of the form given in (9), which is incremented by $\alpha_n^\pm$ for transitions from state $n$ to $n \pm 1$. Repeating the analysis of Section V B 1, we find $\psi_A(s) = -N \min_\rho \mathcal{F}_A(\rho, s)$ with

$$ \mathcal{F}_A(\rho, s) = \frac{1}{N} \left\{ -2(W_{N\rho}^+ e^{-s(\alpha_{N\rho}^+ + \alpha_{N\rho}^-)} W_{N\rho}^-)^{1/2} + W_{N\rho}^+ + W_{N\rho}^- \right\} \quad (96) $$

For example, in the case of the complexity $Q_+$, one has $\alpha_n^\pm = \ln \frac{W_n^\pm}{W_n^+ + W_n^-}$ and

$$ \psi_+(s) = - \min_{\rho} \left\{ -2 \left[ W_{N\rho}^{+} W_{N\rho}^{-} \right]^{\frac{1-s}{2}} \left[ W_{N\rho}^{+} + W_{N\rho}^{-} \right]^s + W_{N\rho}^{+} + W_{N\rho}^{-} \right\} \quad (97) $$

where, again, this variational bound is exact because of its independence of the form of the trivial wavefunction used. By analogy with $\mathcal{F}_K(\rho, s)$, we find that $e^{-Nt\mathcal{F}_A(\rho,s)}$ give the probability of homogeneous histories with density $\rho$, in the s-ensemble.

### C. Dynamical free energy landscape (field theoretic approach)

In Section V B, we have considered large deviations of time-averaged observables, using these quantities to characterise the histories within the s-ensemble. We now discuss the calculation of dynamical correlation functions within this ensemble. We make use of a field-theoretic description of the FA model.

#### 1. Field-theory for the bosonic FA model

Using the Doi-Peliti representation of the bosonic FA model (Section III C), we use coherent states to write the partition function $Z_K(s, t)$ as a path integral over (time-dependent) functions $\{\hat{\phi}_i\}$ and $\{\phi_i\}$ [38]. Then, taking the continuum limit, we promote these functions to fields $(\phi_{x\tau}, \hat{\phi}_{x\tau})$ depending on position $x$ and time $\tau$, where $\phi_{x\tau}$ has the dimensions of a density and $\hat{\phi}_{x\tau}$ is dimensionless. Introducing sources $h$ and $\hat{h}$ for the fields $\phi$ and $\hat{\phi}$, we write

$$ \mathcal{Z}[s, t; h_{x\tau}, \hat{h}_{x\tau}] = \int \mathcal{D}[\phi_{x\tau}, \hat{\phi}_{x\tau}] \exp \left\{ -S_K[\phi, \hat{\phi}] + \int dx d\tau (h_{x\tau} \phi_{x\tau} + \hat{h}_{x\tau} \hat{\phi}_{x\tau}) \right\} \quad (98) $$

where the path integral is over histories of duration $t$, and (see, for example [37])

$$ S_K[\phi_{x\tau}, \hat{\phi}_{x\tau}] = \int dx d\tau \left\{ \hat{\phi}_{x\tau} \partial_t \phi_{x\tau} - 2dl_0^d \left[ (\hat{\phi}_{x\tau} \phi_{x\tau} + cl_0^{-d}) - e^{-s}(\phi_{x\tau} + cl_0^{-d} \hat{\phi}_{x\tau}) \right] (1 + l_0^2 \nabla^2) \hat{\phi}_{x\tau} \phi_{x\tau} \right\} \quad (99) $$

where $l_0$ is the lattice spacing, and we have taken a gradient expansion, truncating at quadratic order. We identify $Z_K(s, t) = \mathcal{Z}[s, t; 0, 0]$.
## 2. Saddle point approximation

We now show that a saddle-point analysis on the action recovers the results of the previous sections. The saddle-point equations are obtained by maximising the action with respect to $\phi$ and $\hat{\phi}$ in the absence of the sources $(h, \hat{h})$. The saddle occurs for fields that are homogeneous in space and time, with magnitudes satisfying

$$0 = 2\hat{\phi}\hat{\phi}(\phi l_0^d - e^{-s}c) + \phi(c - e^{-s}\phi l_0^d), \quad (100)$$

$$0 = 2(\hat{\phi} - e^{-s})\hat{\phi}\phi + c\hat{\phi}(1 - e^{-s}\hat{\phi}). \quad (101)$$

These reduce to a single equation if we take $cl_0^{-d}\hat{\phi} = \phi$ (this origin of this symmetry becomes clear if we use the symmetrised operator $\tilde{W}$ in the construction of the original path integral). In this single variable, the solutions are $\hat{\phi} = 0$ and

$$\hat{\phi} = \hat{\phi}_{\text{act}} \equiv \frac{3}{4}e^{-s} + \frac{1}{4}\sqrt{9e^{-2s} - 8} \quad (102)$$

To estimate the dynamical free energy, we simply identify $\psi_K(s)$ with $(-t^{-1} \min S[\phi, \hat{\phi}])$ where the minimum is over value of the action at the two saddles. The result is

$$\psi_K(s) \approx \begin{cases} 0, & s > 0 \\ Nd(c\hat{\phi}_{\text{act}})^2[\hat{\phi}_{\text{act}}e^{-s} - 1], & s < 0 \end{cases} \quad (103)$$

where the approximate equality indicates that we are working in the saddle-point approximation. We identify the time-dependent density (per site) of excitations in the s-ensemble $\langle\rho(\tau)\rangle_s = l_0^d\langle\phi^*(\tau)\phi(\tau)\rangle$. Away from temporal boundaries, we take the saddle point value for this average, obtaining

$$\rho_K(s) \approx \begin{cases} 0, & s > 0 \\ c(\hat{\phi}_{\text{act}})^2, & s < 0 \end{cases} \quad (104)$$

It is easily verified that (103) and (104) coincide with the variational estimates (58) and (59).

In principle, we can now use the tools of dynamical field theory [49] to incorporate fluctuations around the saddle points, and to calculate spatiotemporal correlation functions in the s-ensemble. For example, defining a density field $n(x, \tau)$ through a continuum limit of the original occupation variables $n_i$, we have

$$\langle n(x, \tau)n(y, \tau') \rangle_s = \frac{\delta^4}{\delta h(x, \tau)\delta h(y, \tau')\delta \hat{h}(x, \tau)\delta \hat{h}(y, \tau')} \ln \mathcal{Z}[s, t; h, \hat{h}] \Big|_{h=\hat{h}=0} \quad (105)$$

Thus, for models with a field-theoretic representation (such as the FA model), the framework described in this section provides methods for systematic calculation of correlation functions and fluctuation effects in the s-ensemble. However, these field-theoretical calculations beyond the scope of this paper. We emphasize that the analysis of Sections IV A and IV C establishes that a dynamical first-order transition does occur at $s = 0$ in finite-dimensional KCMs. Thus, while we expect fluctuations to have quantitative effects, the qualitative picture obtained through this saddle point analysis is not changed.

# VI. OUTLOOK

We have analyzed the dynamics of kinetically constrained models, using an ensemble of histories which span a long time $t$. This analysis used dynamical tools [5, 6, 7, 8, 11] constructed by analogy with the usual Boltzmann-Gibbs theory of equilibrium systems. We have established that this procedure captures physically relevant features that are not accessible from the steady state distribution of configurations in these models.

We have shown that the steady state of KCMs lies on a first-order dynamical transition line, characterised by a coexistence between active and inactive histories. This first-order line is present both in mean-field systems and in finite-dimensional models. Its existence is proven by variational bounds on the dynamical free energy, and confirmed in numerical simulations of several kinetically constrained models, including both spin-facilitated models and kinetically constrained lattice gases. We have defined dynamical Landau-like free energy, whose form is intimately connected to the existence of dynamical heterogeneities.
Earlier studies of non-equilibrium systems used a similar thermodynamic formalism for dynamics to reveal first-order transitions arising from a static phase transition [50, 51] or from an absorbing state [8]. To place our work in context, we emphasise that our dynamical phase coexistence is not related to such phenomena. However, the transitions in these models all appear as singularities in their large deviation functions, consistent with the idea [52] that phase transitions both in and out of equilibrium can be studied through the eigenvalue spectra of their master operators. Moreover, the focus of the current paper is on transitions between stationary, time-reversible dynamical states, and therefore we concentrated on large deviations of quantities that are symmetric in time: this is to be contrasted with studies that have concentrated on currents of entropy or particles [42, 53], although recent work has hinted that large deviations of time-symmetric observables may also be of importance in non-equilibrium steady states [54].

We expect our approach to be meaningful in a wider class of systems than those probed in this paper. For example, glass-forming liquids are known to be dynamically heterogeneous, and this feature can be captured in computational simulations of atomistic models. It would be interesting to establish whether this heterogeneity is linked to a dynamical phase transition similar to that present in KCMs. This could indicate a more general link between glassy properties (not necessarily related to dynamical heterogeneity) and dynamical phase transitions.

Finally, we observe that an experimental scheme for sampling the s-ensemble would be very valuable, since it would provide a direct test for the existence of a dynamical phase transition. However, the fact that the generalised master operator $W_A$ does not conserve probability makes the search for such a scheme rather challenging.

## Acknowledgments

VL would like to thank Julien Tailleur for useful and continuous discussions, and RLJ and JPG thank David Chandler for extensive discussions. This work (EP and FvW) was supported by the French Ministry of Education through Grant No ANR-05-JCJC-44482. JPG was supported by EPSRC under Grant No. GR/S54074/01. While at Berkeley, RLJ was supported initially by NSF grant CHE-0543158 and later by Office of Naval Research Grant No. N00014-07-1-068. VL was supported in part by the Swiss FNS, under MaNEP and division II.

# APPENDIX A: AVERAGES IN THE s-ENSEMBLE AND EIGENVECTORS OF $W_A$

## 1. Eigenvectors of $W_A$

In this appendix, we discuss some properties of the operator $W_A$, and their consequences for averages in the s-ensemble. We write $W_A$ in terms of its left and right eigenvectors $|L_n\rangle$ and $|R_n\rangle$: $W_A = \sum_n \lambda_n |R_n\rangle\langle L_n|$ with eigenvalues $\lambda_0 > \lambda_1 \ge \dots$. The maximal eigenvalue $\lambda_0$ is equal to $\psi_A(s)$. One can normalize eigenvectors so that

$$ \langle L_n | R_m \rangle = \delta_{nm} \quad \text{and} \quad \langle - | R_0 \rangle = 1 \qquad (A1) $$

where $\langle -| = \sum_C \langle C|$ is the projection state.

Thus, for long times, we have $e^{W_A t} = |R_0\rangle\langle L_0|e^{t\psi_A(s)} + \dots$ where the omitted terms on the right hand side are exponentially smaller than the dominant first term. Therefore, starting from an initial state $|P_0\rangle = \sum_C P_0(C)|C\rangle$, with $\langle -|P_0\rangle = 1$, one has, for large times

$$ |P(t)\rangle = e^{W_A t} |P_0\rangle \sim |R_0\rangle e^{\psi_A(s)t} \langle L_0|P_0\rangle + \dots \qquad (A2) $$

where we write the largest eigenvalue of $W_A$ as $\lambda_0 = \psi_A(s)$, and the omitted terms are exponentially smaller than the first one, for large times $t$. This allows us to identify the largest eigenvalue of $W_A$ with the dynamical free energy $\lim_{t\to\infty} t^{-1} \log Z_A(s, t)$, through Equ. (18).

## 2. Time averages

We now consider a configuration-dependent observable $b(C)$, and an s-ensemble defined as in (12), using an observable $A$ of the form given in (9). We provide a link between the eigenvectors of $W_A$ and two weighted averages: the average of $b$ at the final time $t$ in the s-ensemble

$$ \langle b(t) \rangle_s \equiv \frac{\langle b(C(t)) e^{-sA} \rangle}{\langle e^{-sA} \rangle} \qquad (A3) $$
and the time-integrated average of b in the s-ensemble

$$ \langle B \rangle_s \equiv \frac{\left\langle \int_0^t d\tau b(C(\tau))e^{-sA} \right\rangle}{\langle e^{-sA} \rangle} \qquad (A4) $$

As discussed in Section V A, $\langle B \rangle_s$ grows linearly in time, but, in general $\partial_t \langle B \rangle_s \neq \langle b(t) \rangle_s$. In operator notation, the average of b at the final time is

$$ \langle b(t) \rangle_s = \frac{\langle -|\hat{b} e^{t\hat{W}_A} |P_0 \rangle}{\langle -|e^{t\hat{W}_A} |P_0 \rangle} \qquad (A5) $$

where $\hat{b}$ denotes the diagonal operator of elements $b(C)$. Using the normalization (A1), and the large time result (A2), we arrive at

$$ \langle b(t) \rangle_s = \langle -|\hat{b}|R_0 \rangle \qquad (A6) $$

Thus the right eigenvector $|R_0\rangle$ gives the distribution over configurations $C$ at the final time $t$.

On the other hand, the integrated average $\langle B \rangle_s$ is obtained from the mean value $\langle b(\tau) \rangle_s$ in the intermediate regime $0 \ll \tau \ll t$:

$$ \frac{1}{t} \langle B \rangle_s = \langle b(\tau) \rangle_s = \frac{\langle -|e^{(t-\tau)\hat{W}_A} \hat{b} e^{\tau\hat{W}_A} |P_0 \rangle}{\langle -|e^{t\hat{W}_A} |P_0 \rangle} \qquad (A7) $$

For $0 \ll \tau \ll t$, we have

$$ e^{\tau \hat{W}_A} |P_0 \rangle = e^{\tau \psi_A(s)} |R_0 \rangle \langle L_0 | P_0 \rangle \qquad (A8) $$

$$ \langle -|e^{(t-\tau)\hat{W}_A} = \langle -|R_0 \rangle \langle L_0|e^{(t-\tau)\psi_A(s)} \qquad (A9) $$

and hence

$$ \frac{1}{t} \langle B \rangle_s = \langle b(\tau) \rangle_s = \langle L_0 | \hat{b} | R_0 \rangle \qquad (A10) $$

Thus, while the average $\langle b(t) \rangle_s$ depends only on $|R_0\rangle$, the average $\langle B \rangle_s$ depends on both $|R_0\rangle$ and $\langle L_0|$.

### 3. Dynamics with detailed balance

From (26), it follows that if a system obeys detailed balance, its master operator satisfies $\hat{W}_K^\dagger = \hat{P}_{\text{eq}}^{-1}\hat{W}_K\hat{P}_{\text{eq}}$, where $\hat{P}_{\text{eq}}$ is a diagonal operator with elements $P_{\text{eq}}(C)$. Thus, $|L_n\rangle = \hat{P}_{\text{eq}}^{-1}|R_n\rangle$. Using this property together with results from the previous section, and denoting $|R_0\rangle = \sum_C R_0(C, s)|C\rangle$ we write

$$ \langle b(t) \rangle_s = \sum_C b(C) R_0(C, s) \qquad (A11) $$

$$ \partial_t \langle B \rangle_s = \sum_C b(C) \frac{R_0(C, s)^2}{P_{\text{eq}}(C)} \qquad (A12) $$

Clearly, these averages are not the same in general. Expanding about $s = 0$ and using $R_0(C, 0) = P_{\text{eq}}(C)$, we arrive at (73), with

$$ b^{(1)} = \sum_C b(C) \frac{\partial R_0(C, s)}{\partial s} \qquad (A13) $$

Finally, we note that expectation values of the form $\partial_t \langle B \rangle_s$ take a simple form when written in terms of the eigenvectors $|V\rangle$ of the symmetric operator $\hat{W}_K$, discussed in Section II B 6. The matrix elements of this operator are $(\hat{W}_K)_{C,C'} = P_{\text{eq}}^{-1/2}(C)(\hat{W}_K)_{C,C'} P_{\text{eq}}^{1/2}(C)$, so its eigenvectors are $V_n(C, s) \propto \sqrt{L_n(C, s)R_n(C, s)} = R_n(C, s)/\sqrt{P_{\text{eq}}(C)}$. Thus, we have (for large time)

$$ \frac{1}{t} \langle B \rangle_s = \frac{\sum_C b(C) V_0(C, s)^2}{\sum_C V_0(C, s)^2} \qquad (A14) $$

which links the eigenvector $V_0(C, s)$ to physical observables such as B.
APPENDIX B: OBSERVABLES OF TYPES A AND B

Here, we discuss the connections between observables of the forms given in (9) and (22): we refer to these observables as type A and type B respectively. We begin with a result that is used in the numerical methods of [19, 20].

Consider an s-ensemble defined as in Section II B 3. That is, take a system with rates $W(C \to C')$ and modify the statistical weights of its histories by a factor $e^{-sA}$, where A is an observable of type A. In addition, we define a second stochastic process ('modified dynamics') through the transition rates

$$W_s(C \to C') = e^{-s\alpha(C,C')} W(C \to C') \quad (B1)$$

where the $\alpha(C, C')$ are obtained from the definition of the observable A, through (9). In addition we define two configuration-dependent observables, $r_s(C) = \sum_{C'} W_s(C \to C')$, and

$$\delta r_s(C) = r_s(C) - r(C) \quad (B2)$$

[Here, $r(C) = \sum_{C'} W(C \to C')$ is the escape rate for the dynamics W, as in the main text.]

Motivated by the decomposition of Equ. (61), we can establish two ways of defining the same s-ensemble. From (7), it is easily verified that

$$\mathrm{Prob}[\mathrm{hist}|W]e^{-sA} = \mathrm{Prob}[\mathrm{hist}|W_s]e^{\delta R_s} \quad (B3)$$

where the notation $\mathrm{Prob}[\mathrm{hist}|W]$ refers to the (unmodified) probability of a history in a system with dynamical rates W, and

$$\delta R_s = \int_0^t d\tau \delta r_s(C(\tau)) \quad (B4)$$

Thus, Equ. (B3) states that histories in the s-ensemble parameterized by A for the original dynamics W have the same weight as histories in an s-ensemble parameterized by $\delta R_s$, for the modified dynamics $W_s$.

Since the two ensembles are identical, it follows that all observables have the same averages: for example

$$\langle Oe^{-sA} \rangle_W = \langle Oe^{\delta R_s} \rangle_{W_s}, \quad (B5)$$

where the subscript on the average refers to the dynamical rules used for the sum over histories. Further, this result holds for histories of finite duration $t$, as long as the same initial conditions are used in both averages.

For the specific case where the observable A is the activity K then this relation takes a particularly simple form. Following Section II B 3 with $\alpha(C, C') = 1$ for all C and C', the master operator associated with this s-ensemble has matrix elements $(W_K)_{C,C'} = e^{-s} W(C \to C) - r(C)\delta_{C,C'}$. From (B1), we find $W_s(C \to C') = e^{-s}W(C \to C')$: that is, the modification to the dynamics simply involves a rescaling of time by a factor $e^{-s}$. In addition, for $B = R$, we have $\delta r_s(C) = sr(C)$, so we define an s-ensemble associated with the observable $R[\mathrm{hist}] = \int_0^t d\tau r(C(\tau))$, which is of type B. From the analysis of Section II B 5, the master operator associated with this ensemble, $W_R$, has matrix elements

$$(W_R)_{C,C'} = W(C' \to C) - (1+s)r(C)\delta_{C,C'} \quad (B6)$$

from which we can see that

$$W_K(s) = e^{-s}W_R(e^s - 1) \quad (B7)$$

This equation relates the dynamical free energies of the s-ensembles for K and R, based on the same unbiased dynamics W. The dynamical free energies $\psi_K(s)$ and $\phi_R(s)$ are given by the largest eigenvalues of $W_K$ and $W_R$: they satisfy

$$\psi_K(s) = e^{-s}\phi_R(e^s - 1) \quad (B8)$$

Hence, we can also relate the cumulants of the observables K and R. For example,

$$\langle K \rangle = \langle R \rangle \quad (B9)$$

$$\langle K^2 \rangle_c = \langle R^2 \rangle_c - \langle R \rangle \quad (B10)$$

This last equation provides an interpretation of the variance (second cumulant) of K, through

$$\frac{1}{t} \langle K^2 \rangle_c = -\langle r \rangle + \int_0^t \frac{2t'}{t} \left[ \langle r(t)r(t') \rangle - \langle r(t) \rangle \langle r(t') \rangle \right] dt' \quad (B11)$$

where the correlation function is evaluated at $s = 0$.
APPENDIX C: LINK TO DONSKER-VARADHAN THEORY

As in the main text, we consider a Markov process described by transition rates $W(C \to C')$ between configurations $\{C\}$. For a history $C(\tau)$, we define the experimental measure

$$ \bar{\mu}(C, t) = \int_{0}^{t} d\tau \, \delta_{C, C(\tau)} \qquad (C1) $$

This history-dependent observable simply counts how much time was spent in configuration $C$ between 0 and $t$. This is the central object of Donsker-Varadhan [55] theory (see also [5, 56]). For large times, the experimental measure approaches the steady state distribution, $\lim_{t \to \infty} \frac{1}{t} \langle \bar{\mu}(C, t) \rangle = P_{\text{st}}(C)$.

Donsker-Varadhan theory gives information on the large deviations of the experimental measure $\bar{\mu}(C, t)$ from the steady-state distribution, in the long time limit. Therefore, it is naturally connected to the statistics of histories and to the dynamical ensemble approach discussed in this article. For example, consider an observable $b(C)$ depending on the configuration of the system. The experimental measure $\bar{\mu}(C, t)$ determines the time-integrated value of the observable $b$ through

$$ \int_{0}^{t} d\tau \, b(C(\tau)) = \sum_{C} b(C) \bar{\mu}(C, t) \equiv t \langle b \rangle_{\bar{\mu}} \qquad (C2) $$

which defines $\langle b \rangle_{\bar{\mu}}$: the average of $b$ with respect to the experimental measure $\bar{\mu}$. In this appendix, we establish links between the s-ensemble approach and the results of Donsker and Varadhan. In particular, we develop a variational method that gives the large deviations of an observable $B$, in the s-ensemble defined for an (unrelated) observable $A$.

1. Donsker-Varadhan large deviation function

The Donsker-Varadhan (DV) theorem [55] states that in the long time limit

$$ \mathrm{Prob}[\bar{\mu}(C, t) = t\mu(C)] = e^{tJ[\mu]} \qquad (C3) $$

with (see for instance [5, 56])

$$ J[\mu] = \inf_{\rho > 0} \sum_{C, C'} \left\{ W(C \to C') \frac{\rho(C')}{\rho(C)} \mu(C) - r(C)\mu(C)\delta_{C'C} \right\} \qquad (C4) $$

where the infimum has to be taken over normalized measures $\rho(C)$, with $\sum_C \rho(C) = 1$.

If $\mathbb{W}$ obeys detailed balance with respect to an equilibrium distribution $P_{\text{eq}}(C)$, the infimum is obtained for $\rho(C) = \sqrt{\mu(C)/P_{\text{eq}}(C)}$ and the large deviation function reduces to

$$ J_{\text{eq}}[\mu] = \sum_{C, C'} \left\{ \left[ W(C \to C')W(C' \to C) \right]^{1/2} [\mu(C)\mu(C')]^{1/2} - r(C)\mu(C)\delta_{C'C} \right\} \qquad (C5) $$

for normalised measures $\sum_C \mu(C) = 1$. Writing

$$ \mu(C) = \frac{V(C)^2}{\sum_C V(C)^2}, \qquad (C6) $$

we identify

$$ J_{\text{eq}}[\mu] = \frac{\langle V | \tilde{W} | V \rangle}{\langle V | V \rangle} \qquad (C7) $$

as the function to be maximised in (27), for the case $s = 0$.
## 2. Dynamical Landau free energy at $s = 0$

We now apply the DV theorem to the large deviations of an observable $b(C)$. Integrating over a time $t$, we define the history-dependent quantity

$$B(t) = \int_{0}^{t} d\tau b(C(\tau)), \qquad (C8)$$

As discussed in Section II B 4, one expects the probability distribution of $B(t)$ to behave as

$$\Omega_{\text{dyn}}(B = bt, t) \sim e^{t\pi(b)} \qquad (C9)$$

for large times $t$.

The large-deviation function $\pi(b)$ can be obtained through the Donsker-Varadhan functional using $B(t) = \sum_C b(C) \bar{\mu}(C, t)$, so that

$$\begin{aligned} \text{Prob}[B(t) = tb] &= \langle \delta(B(t) - tb) \rangle \\ &= \int d\mu \langle \delta(B(t) - tb) \delta(t^{-1}\bar{\mu}(C, t) - \mu(C)) \rangle \\ &= \int d\mu \delta\left(\sum_C b(C)\mu(C) - b\right) \langle \delta(\bar{\mu}(C, t) - t\mu(C)) \rangle \end{aligned} \qquad (C10)$$

$$= \int d\mu \delta(\langle b \rangle_{\mu} - b) e^{tJ[\mu]} \qquad (C11)$$

where the average $\langle b \rangle_{\mu}$ was defined in (C2). Here, we have replaced an average over histories $\langle \cdot \rangle$ with an integral over possible realisations of the experimental measure $\mu$, weighted by their probabilities (which are known from the DV theorem). In the limit of large time, we maximise the argument of the exponential, subject to a constraint imposed by the $\delta$-function. Hence,

$$\pi(b) = \sup_{\substack{\mu \text{ with} \\ \langle b \rangle_{\mu} = b}} J[\mu] \qquad (C12)$$

which for systems obeying detailed balance can again be expressed in terms of the operator $\tilde{W}$, using (C6).

## 3. Dynamical Landau free energy for any $s$

We now generalise this analysis to the $s$-ensemble. We note that the values of 'type A' observables [those of the form given in (9)] cannot be obtained from the experimental measure $\bar{\mu}(C)$. To connect these observables to the DV approach, we use the results of appendix B.

The large deviations of the observable $B$ in the $s$-ensemble specified by $A$ are determined by

$$\Omega_A(s, b) = \langle \delta(B - bt)e^{-sA} \rangle_W, \qquad (C13)$$

where as in appendix B, the label on the average indicates the dynamical rules used to generate the ensemble of histories. From (B5), we can write

$$\Omega_A(s, b) = \langle \delta(B - bt)e^{-sA} \rangle_W = \langle \delta(B - bt)e^{\delta R_s} \rangle_{W_s}, \qquad (C14)$$

with an observable $\delta R_s$ and rates $W_s(C \to C')$ given in Equ. (B1) and (B4).

Now, following the analysis of the previous section, we have

$$\Omega_A(s, b) = \left\langle e^{\delta R_s} \delta(\bar{B}(t) - tB) \right\rangle_{W_s} = \int d\mu \delta(\langle b \rangle_{\mu} - B) e^{t(J[\mu|W_s] + \langle \delta r_s \rangle_{\mu})} \qquad (C15)$$

where $\langle \delta r_s \rangle_{\mu} = \sum_C \mu(C) \delta r(C)$, and $J[\mu|W_s]$ is the Donsker-Varadhan functional for stochastic process with rates $W_s$. Again, the integral over $\mu$ can be evaluated by maximising the argument of the exponential subject to the constraint on $\langle b \rangle_{\mu}$, leading to $\Omega_A(s, b) \sim \exp t\pi_A(s, b) = \text{with}$

$$\pi_A(b, s) = - \sup_{\substack{\mu \text{ with} \\ \langle b \rangle_{\mu} = b}} J_A[\mu, s] \qquad (C16)$$
with

$$J_A[\mu, s] = \sum_{C,C'} \left\{ e^{-\frac{s}{2}[\alpha(C,C') + \alpha(C',C)]} [W(C \to C')W(C' \to C)]^{1/2} [\mu(C)\mu(C')]^{1/2} - r(C)\mu(C)\delta_{C'C} \right\} \quad (C17)$$

where we emphasise that the rates $W(C \to C')$ are those of the original (unmodified) dynamics. From (C6), we identify

$$J_A[\mu, s] = \frac{\langle V | \tilde{W}_A | V \rangle}{\langle V | V \rangle} \quad (C18)$$

as the quantity to be maximised in (27) for $A = K$. Moreover, for $b$ being the occupation number $n$ and $A$ the activity $K$, one recognizes in (C18) the result (94).

We observe that these results have been derived for dynamics which obey detailed balance, but they are not restricted to that situation. For instance, (C16) holds in general, with

$$J_A[\mu, s] = \inf_{\rho>0} \sum_{C,C'} \left\{ e^{-s\alpha(C,C')} W(C \to C') \frac{\rho(C')}{\rho(C)} \mu(C) - r(C)\mu(C)\delta_{C'C} \right\} \quad (C19)$$

Finally, we note that these results amount to a generalization of the Donsker-Varadhan theorem (C3) in the $s$-ensemble: for large times,

$$\langle e^{-sA} \delta(\bar{\mu}(C, t) - t\mu(C)) \rangle = e^{t J_A[\mu, s]} \quad (C20)$$

with $J_A[\mu, s]$ given in general by (C19), which reduces to (C17) if $W_A$ can be symmetrised.

[1] M.D. Ediger, C.A. Angell and S.R. Nagel, J. Phys. Chem. **100**, 13200 (1996).

[2] C.A. Angell, Science **267**, 1924 (1995).

[3] P.G. Debenedetti and F.H. Stillinger, Nature **410**, 259 (2001).

[4] K. Binder and W. Kob, "Glassy Materials and disordered solids", (World Scientific, 2005).

[5] For a recent review see H. Touchette, arXiv:0804.0327 (to appear in Phys. Rep.).

[6] D. Ruelle, "Thermodynamic Formalism" (Addison-Wesley, Reading, 1978); J.-P. Eckmann and D. Ruelle, Rev. Mod. Phys. **57**, 617 (1985).

[7] P. Gaspard, "Chaos, scattering and statistical mechanics" (CUP, Cambridge, 1998).

[8] V. Lecomte, C. Appert-Rolland and F. van Wijland, Phys. Rev. Lett. **95**, 010601 (2005); J. Stat. Phys. **127** 51 (2007);

[9] F. Ritort and P. Sollich, Adv. Phys. **52**, 219 (2003).

[10] J.P. Garrahan, R.L. Jack, V. Lecomte, E. Pitard, K. van Duijvendijk and F. van Wijland, Phys. Rev. Lett., **98** 195702 (2007).

[11] M. Merolle, J.P. Garrahan and D. Chandler, Proc. Natl. Acad. Sci. USA **102**, 10837 (2005); R.L. Jack, J.P. Garrahan and D. Chandler, J. Chem. Phys. **125**, 184509 (2006).

[12] V. Lecomte, C. Appert-Rolland and F. van Wijland, Cr. Acad. Sc. Paris **8**, 609 (2007).

[13] M. Mézard and G. Parisi, Phys. Rev. Lett. **82**, 747 (1999).

[14] S. Franz and G. Parisi, J. Phys. C **12**, 6335 (2000).

[15] X. Xia and P.G. Wolynes, Proc. Natl. Acad. Sci. USA **97**, 2990 (2000).

[16] J.P. Bouchaud and G. Biroli, J. Chem. Phys. **121**, 7347 (2004).

[17] G. Tarjus, S.A Kivelson, Z. Nussinov and P. Viot, J. Phys. Condens. Matter **17**, R1143 (2005).

[18] W. Götze, J. Phys. Condens. Matter **11**, A1 (1999).

[19] C. Giardinà, J. Kurchan, and L. Peliti, Phys. Rev. Lett. **96**, 120603 (2006).

[20] V. Lecomte and J. Tailleur, J. Stat. Mech. (2007) P03004.

[21] R. Kubo, H. Ichimura, T. Usui and N. Hashitsume, "Thermodynamics", (North-Holland, Amsterdam, 1968).

[22] G.H. Fredrickson and H.C. Andersen, Phys. Rev. Lett. **53**, 1244 (1984).

[23] J. Jäckle and S. Eisinger, Z. Phys. B **84**, 115 (1991).

[24] W. Kob and H.C. Andersen, Phys. Rev. E **48**, 4364 (1993).

[25] J. Jäckle and A. Krönig, J. Phys.: Condens. Matter **6**, 7633 (1994); A. Krönig and J. Jäckle, J. Phys.: Condens. Matter **6**, 7655 (1994).

[26] J. Kurchan, L. Peliti, and M. Sellitto, Europhys. Lett. **39**, 365 (1997).

[27] P. Sollich and M.R. Evans, Phys. Rev. Lett. **83**, 3238 (1999).

[28] M. Einax and M. Schulz, J. Chem. Phys. **115**, 2282 (2001).

[29] J.P. Garrahan and D. Chandler, Phys. Rev. Lett. **89**, 035704 (2002).
[30] D. Aldous and P. Diaconis, J. Stat. Phys. **107**, 945 (2002).

[31] Y. Jung, J. P. Garrahan and D. Chandler, Phys. Rev. E **69**, 061205 (2004).

[32] C. Toninelli, G. Biroli and D.S. Fisher, J. Stat. Phys. **120**, 167 (2005).

[33] A. Pan, J. P. Garrahan and D. Chandler, Phys. Rev. E **72**, 041106 (2005).

[34] P.L. Geissler and D.R. Reichman, Phys. Rev. E **71**, 031206 (2005).

[35] L. Berthier and J.P. Garrahan, J. Phys. Chem. B **109**, 3578 (2005).

[36] S. Whitelam, L. Berthier and J.P. Garrahan, Phys. Rev. Lett. **92**, 185705 (2004).

[37] R. Jack, P. Mayer and P. Sollich, J. Stat. Mech. (2006) P03006.

[38] M. Doi, J. Phys. A **9**, 1465 (1976); L. Peliti, J. Physique **46**, 1469, (1985).

[39] J.L. Cardy and U.C. Täuber, Phys. Rev. Lett. **77**, 4780 (1996); J. Stat. Phys. **90**, 1 (1998).

[40] P. Gaspard and X.-J. Wang, Phys. Rep. **235** 291 (1993).

[41] P. Allegrini, J.F. Douglas, S.C. Glotzer, Phys. Rev. E **60**, 5714 (1999).

[42] J.L. Lebowitz and H. Spohn, J. Stat. Phys. **95**, 333 (1999).

[43] H. Hinrichsen, Adv. Phys. **49**, 815 (2000).

[44] J. Tailleur and J. Kurchan, Nature Physics **3**, 203 (2007).

[45] J.B. Anderson, J. Chem. Phys. **63** 1499 (1975).

[46] R. Stinchcombe, Adv. Phys. **50**, 431 (2001).

[47] P. Jordan and E. Wigner, Z. Phys. **47**, 631 (1928).

[48] M. Henkel, E. Orlandini and J. Santos, Ann. Phys. **259**, 163 (1997).

[49] J. Zinn-Justin, "Quantum Field Theory and Critical Phenomena" (Oxford University Press, 1989).

[50] F. Ritort, J. Stat. Mech. (2004) P10016.

[51] A. Imparato and L. Peliti, Phys. Rev. E **72**, 046114 (2005).

[52] B. Gaveau and L. S. Schulman, J. Math. Phys. **39**, 1517 (1998).

[53] J. Kurchan, J. Phys. A **31**, 3719 (1998).

[54] T. Bodineau and F. Lefevere, J. Stat. Phys. **133**, 1 (2008); C. Maes, K. Netočný and B. Wynants, Physica A **387**, 2675 (2008); C. Maes and K. Netočný, Europhys. Lett. **82**, 30003 (2008).

[55] M. D. Donsker, S. R. Varadhan, Comm. Pure Appl. Math. 28:1 (1975).

[56] C. Maes, K. Netočný, B. Shergelashvili, *A selection of nonequilibrium issues*, Lecture notes from the 5th Prague Summer School on Mathematical Statistical Mechanics (2006); math-ph/0701047.