Exact large deviation statistics and trajectory phase transition of a deterministic boundary driven cellular automaton

TL;DR

This paper provides an exact analysis of the large deviation statistics and phase transition behavior in a boundary-driven deterministic cellular automaton, revealing coexistence of active and inactive phases and explicit dynamics for rare events.

Contribution

It introduces an exact matrix product approach to compute large deviation functions and characterizes phase coexistence in a deterministic CA with boundary driving.

Findings

Phase coexistence between active and inactive dynamical phases

Exact finite size scaling of trajectory transitions

Explicit Doob-transformed dynamics for rare events

Abstract

We study the statistical properties of the long-time dynamics of the rule 54 reversible cellular automaton (CA), driven stochastically at its boundaries. This CA can be considered as a discrete-time and deterministic version of the Fredrickson-Andersen kinetically constrained model (KCM). By means of a matrix product ansatz, we compute the exact large deviation cumulant generating functions for a wide range of time-extensive observables of the dynamics, together with their associated rate functions and conditioned long-time distributions over configurations. We show that for all instances of boundary driving the CA dynamics occurs at the point of phase coexistence between competing active and inactive dynamical phases, similar to what happens in more standard KCMs. We also find the exact finite size scaling behaviour of these trajectory transitions, and provide the explicit "Doob-transformed" dynamics that optimally realises rare dynamical events.

Equations

$$M_{\mathrm{e}}$$

$$M_{\mathrm{o}}$$

$${\mathcal{O}}_{T}=\sum_{t=0}^{T-1}\sum_{j=1}^{N-1}\Big{[}f_{j}(n^{t}_{j},n^{t}_{j+1})+g_{j}(n^{t+1/2}_{j},n^{t+1/2}_{j+1})\Big{]}$$

$$f_{n{n}^{\prime }}^{(j-1)}{f}_{n^{\prime }{n}^{\prime }{}^{\prime }}^{(j)}{W}_n^{(j-1)}{W}_{n^{\prime }}^{(j)}{X}_{n^{\prime }{}^{\prime }}^{(j+1)}=X_n^{(j-1)}{V}_{\chi (n{n}^{\prime }{n}^{\prime }{}^{\prime })}^{(j)}{V}_{n^{\prime }{}^{\prime }}^{(j+1)},$$

$$g_{n{n}^{\prime }}^{(j-2)}{g}_{n^{\prime }{n}^{\prime }{}^{\prime }}^{(j-1)}{X}_n^{(j-2)}{V}_{n^{\prime }}^{(j-1)}{V}_{n^{\prime }{}^{\prime }}^{(j)}=W_n^{(j-2)}{W}_{\chi (n{n}^{\prime }{n}^{\prime }{}^{\prime })}^{(j-1)}{X}_{n^{\prime }{}^{\prime }}^{(j)},$$

$$f_{n{n}^{\prime }}^{(1)}{f}_{n^{\prime }{n}^{\prime }{}^{\prime }}^{(2)}⟨l_n|W_{n^{\prime }}^{(2)}{X}_{n^{\prime }{}^{\prime }}^{(3)}$$

$$\sum _{m,m^{\prime }=0,1}{R}_{n{n}^{\prime }}^{m{m}^{\prime }}{f}_{m{m}^{\prime }}^{(N-1)}|r_{m{m}^{\prime }}⟩$$

$$\sum _{m,m^{\prime }=0,1}{L}_{n{n}^{\prime }}^{m{m}^{\prime }}{g}_{m{m}^{\prime }}^{(1)}⟨l_{m{m}^{\prime }}^{\prime }|$$

$$g_{n{n}^{\prime }}^{(N-2)}{g}_{n^{\prime }{n}^{\prime }{}^{\prime }}^{(N-1)}{X}_n^{(N-2)}{V}_{n^{\prime }}^{(N-1)}|r_{n^{\prime }{}^{\prime }}^{\prime }⟩$$

$$p_{n_1,\dots ,n_N}$$

$$p_{n_1,\dots ,n_N}^{\prime }$$

$${\lambda }^4-\alpha \gamma {\mathfrak{a}}_N{\lambda }^3-\omega {\mathfrak{a}}_N^2{\lambda }^2-\beta \delta \xi {\mathfrak{a}}_N^3\lambda +\eta {\mathfrak{a}}_N^4=0,$$

$$\omega ={\mathfrak{b}}_N(1-\alpha )(1-\delta ){\beta }^{\prime }{\gamma }^{\prime }+{\mathfrak{c}}_N(1-\beta )(1-\gamma ){\alpha }^{\prime }{\delta }^{\prime }$$

$$\xi ={\mathfrak{b}}_N{\mathfrak{c}}_N(1-\alpha )(1-\beta )(1-\gamma )(1-\delta )\frac{{\alpha }^{\prime }{\beta }^{\prime }{\gamma }^{\prime }{\delta }^{\prime }}{\tilde{\alpha }\tilde{\beta }\tilde{\gamma }\tilde{\delta }}$$

$$\eta =(\alpha \beta -\tilde{\alpha }\tilde{\beta })(\gamma \delta -\tilde{\gamma }\tilde{\delta })\xi$$

$${\mathfrak{a}}_N$$

$${\mathfrak{b}}_N$$

$${\mathfrak{c}}_N$$

$$\tilde{\alpha }=\frac{f_{10}^{(1)}{g}_{11}^{(1)}}{f_{00}^{(1)}{g}_{01}^{(1)}}(1-\alpha ),\tilde{\beta }=\frac{f_{11}^{(1)}{g}_{10}^{(1)}}{f_{01}^{(1)}{g}_{00}^{(1)}}(1-\beta ),$$

$$\tilde{\gamma }=\frac{f_{11}^{(N-1)}{g}_{01}^{(N-1)}}{f_{10}^{(N-1)}{g}_{00}^{(N-1)}}(1-\gamma ),\tilde{\delta }=\frac{f_{01}^{(N-1)}{g}_{11}^{(N-1)}}{f_{00}^{(N-1)}{g}_{10}^{(N-1)}}(1-\delta ).$$

$$\theta_{N}(s)={\vartheta}(Ns)+{\mathcal{O}}\Big{(}\frac{1}{N}\Big{)},$$

$$0$$

$$\underset{T\to \infty }{\lim }\frac{1}{T}⟨⟨{\mathcal{O}}_T^k⟩⟩={(-{)}^k\frac{d^k}{d{s}^k}{\theta }_N|}_{s=0}\propto N^k$$

$$\lim_{T\to\infty}\frac{1}{TN}\langle{\mathcal{O}}_{T}\rangle=a+\frac{\mu b+\nu c}{2(\mu+\nu)+\alpha\gamma-\beta\delta}+\mathcal{O}\Big{(}\frac{1}{N}\Big{)}$$

$$\underset{T\to \infty }{\lim }\frac{1}{TN}\mathrm{var}{\mathcal{O}}_T=N[-\frac{2bc(1-\alpha \gamma )+\mu b^2+\nu c^2}{2(\mu +\nu )+\alpha \gamma -\beta \delta }$$

$$+\frac{3(\mu b+\nu c{)}^2}{(2(\mu +\nu )+\alpha \gamma -\beta \delta {)}^2}+\frac{2(b+c)(\mu b+\nu c)(2-\alpha \gamma )}{(2(\mu +\nu )+\alpha \gamma -\beta \delta {)}^2}$$

$$-\frac{2(\mu b+\nu c{)}^2(4+\mu +\nu -\alpha \gamma )}{(2(\mu +\nu )+\alpha \gamma -\beta \delta {)}^3}]+\mathcal{O}(1)$$

$$\vartheta (\sigma )=\{\begin{array}{l}-a\sigma +\ln (\alpha \gamma )+o(1),\sigma \to \infty \\ -\frac{1}{3}(b+c+3a)\sigma +\frac{1}{3}\ln (\beta \delta )+o(1),\sigma \to -\infty \end{array}$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}{\theta }_N(s)=\{\begin{array}{l}-as,s>0\\ -\frac{1}{3}(b+c+3a)s,s<0\end{array}$$

$$M_{\mathrm{Doob}}=\frac{1}{\lambda (s)}\mathcal{L}M(s){\mathcal{L}}^{-1},$$

$$M=M_{\mathrm{o}}{M}_{\mathrm{e}}$$

$$M_{\mathrm{e}}=P_{123}{P}_{345}\cdots P_{N-3,N-2,N-1}{R}_{N-1,N}$$

$$M_{\mathrm{o}}=L_{12}{P}_{234}{P}_{456}\cdots P_{N-2,N-1,N}$$

$$P_{i-1,i,i+1}={\text{1}}_2^{\otimes (i-2)}\otimes P\otimes {\text{1}}_2^{\otimes (N-i-1)},L_{12}=L\otimes {\text{1}}_2^{\otimes (N-2)},R_{N-1,N}={\text{1}}_2^{\otimes (N-2)}\otimes R$$

$$P=\left(\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right)$$

$$L=\left(\begin{array}{cccc}\alpha & 0& \alpha & 0\\ 0& \beta & 0& \beta \\ 1-\alpha & 0& 1-\alpha & 0\\ 0& 1-\beta & 0& 1-\beta \end{array}\right),R=\left(\begin{array}{cccc}\gamma & \gamma & 0& 0\\ 1-\gamma & 1-\gamma & 0& 0\\ 0& 0& \delta & \delta \\ 0& 0& 1-\delta & 1-\delta \end{array}\right).$$

$$M(s)=M_{\mathrm{o}}G(s)M_{\mathrm{e}}F(s).$$

Full text

Introduction.– Classical systems which evolve stochastically subject to constraints display complex dynamics, often beyond what can be anticipated simply from their static properties. This is what occurs in the presence of excluded volume interactions, such as in simple exclusion processes Derrida (2007); Mallick (2015), or when configuration space is restricted, such as in dimer coverings Henley (2010); Chalker (2017), or in systems where dynamical rules are subject to kinetic constraints, as for example in kinetically constrained models (KCMs) of glasses Ritort and Sollich (2003); Garrahan et al. (2011). Constrained dynamics is also proving increasingly relevant to quantum many-body systems, including problems such as slow thermalisation and non-ergodicity in the absence of disorder van Horssen et al. (2015); Smith et al. (2017); Shiraishi and Mori (2017); Lan et al. (2018); Turner et al. (2018), operator spreading and entanglement growth Nahum et al. (2017); von Keyserlingk et al. (2018); Rowlands and Lamacraft (2018); Chen and Zhou (2018); Gopalakrishnan (2018); Knap (2018); Tran et al. (2018); Gopalakrishnan et al. (2018), and in the dynamics of ensembles of Rydberg atoms Lesanovsky and Garrahan (2013); Urvoy et al. (2015); Valado et al. (2016).

Complex collective dynamics must be characterised through the statistical properties of dynamical observables, something which can be readily done by means of large deviation (LD) techniques Lecomte et al. (2007); Garrahan et al. (2007, 2009); Touchette (2009). This allows to study ensembles of trajectories of the dynamics as one would study ensembles of configurations in equilibrium statistical mechanics. Among other things, the dynamical LD approach reveals in many systems the existence of competing dynamical phases and the corresponding phase transitions between them, as for example in KCMs Garrahan et al. (2007, 2009), exclusion processes Appert-Rolland et al. (2008); Espigares et al. (2013); Jack et al. (2015); Karevski and Schütz (2017), dimer models Oakes et al. (2018), and several other classical Hedges et al. (2009); Speck et al. (2012); Weber et al. (2013); Baek et al. (2017) and quantum Garrahan and Lesanovsky (2010) systems. This rich phase structure of trajectory space is what underlies the complex dynamics of these systems.

Here we generalise the above ideas to systems whose (bulk) dynamics is deterministic and reversible. We consider specifically the “rule 54” cellular automaton (CA) of Ref.  Bobenko et al. (1993) (RCA54). The local rules that define the interactions of this CA (see below) are similar to the kinetic constraints of the simplest of KCMs, the (one-spin facilitated) Fredrickson-Andersen (FA) model Fredrickson and Andersen (1984); Ritort and Sollich (2003); Garrahan (2018). As such the RCA54 is referred to also as the “Floquet-FA” model Gopalakrishnan (2018) since it can be considered a synchronous, discrete and deterministic version of the FA model. (The RCA54 is also related to the ERCA 250R of Takesue Takesue (1987).) A remarkable property of the RCA54 is that it is integrable Prosen and Mejía-Monasterio (2016) and, in presence of stochastic driving at its boundaries, one can obtain exactly its (in general non-equilibrium) steady state distribution Prosen and Mejía-Monasterio (2016); Inoue and Takesue (2018), certain decay modes Prosen and Buča (2017), and dynamical structure factors Klobas et al. (2018) in terms of matrix product states.

In this paper we compute the exact large deviation statistics of the boundary driven RCA54 by generalising the methods of Refs. Prosen and Mejía-Monasterio (2016); Inoue and Takesue (2018); Prosen and Buča (2017); Klobas et al. (2018). Via a novel inhomogeneous matrix product ansatz we obtain the exact cumulant generating functions and rate functions of a broad class of time-extensive observables of the dynamics. We prove the existence of distinct active and inactive dynamical phases, with the dynamics of the RCA54 occurring at the phase transition point. To our knowledge, our findings here represent the only exact results for LDs in interacting models beyond those for simple exclusion processes Derrida and Lebowitz (1998); Appert-Rolland et al. (2008); Prolhac (2010); de Gier and Essler (2011); Gorissen et al. (2012); Crampé et al. (2016), and the first for bulk-deterministic systems.

Model.– We consider a system defined by binary variables $n_{i}\in\{0,1\}$ (up/down state) on sites $i\in\{1,\ldots,N\}$ of a lattice with even size $N$. A configuration at time $t$ is described by a binary string $\bm{n}^{t}=(n_{1}^{t},n_{2}^{t},\dots,n_{N}^{t})$. The dynamics in the bulk is given by the discrete, deterministic RCA54 Bobenko et al. (1993), while the dynamics on the boundary sites is stochastic Prosen and Mejía-Monasterio (2016); Prosen and Buča (2017).

The update rule is decomposed into two Floquet-like half-time steps. During the first half-time step, $\bm{n}^{t}\to\bm{n}^{t+1/2}$, only even sites are updated, so that $n_{i}^{t+1/2}=n_{i}^{t}$ for $i$ odd. For all $i$ even with $2\leq i\leq N-2$, the evolution of $n_{i}^{t}$ is deterministic through the relation $n_{i}^{t+1/2}=\chi(n_{i-1}^{t},n_{i}^{t},n_{i+1}^{t})$, where $\chi(n,n^{\prime},n^{\prime\prime})=n+n^{\prime}+n^{\prime\prime}+nn^{\prime\prime}\mod 2$ is the rule-54 function Bobenko et al. (1993), see Fig. 1(a). This local update rule is similar to the constraint of the FA model: a site can flip only if at least one of its nearest neighbours is in the up state Fredrickson and Andersen (1984). The last site is updated stochastically depending on the state of its neighbour: $n_{N}^{t+1/2}=0$ with probability $\gamma+n_{N-1}^{t}(\delta-\gamma)$ or $n_{N}^{t+1/2}=1$ otherwise. In the second half-time step, $\bm{n}^{t+1/2}\to\bm{n}^{t+1}$, for even sites we have $n_{i}^{t+1}=n_{i}^{t+1/2}$, while odd sites in the bulk evolve deterministically with rule-54, $n_{i}^{t+1}=\chi(n_{i-1}^{t+1/2},n_{i}^{t+1/2},n_{i+1}^{t+1/2})$. The first site is updated stochastically, with $n_{1}^{t+1}=0$ with probability $\alpha+n_{2}^{t+1/2}(\beta-\alpha)$ or $n_{1}^{t+1}=1$ otherwise.

The above rules define a discrete-time, irreducible and non-reversible Markov process. Physically it models a gas of solitons stochastically emitted from reservoirs at the boundaries, propagating at constant unit velocity in the bulk and interacting pairwise through a one-time step delay Prosen and Mejía-Monasterio (2016); Prosen and Buča (2017); Inoue and Takesue (2018). Depending on the boundary rates, the system is driven out-of-equilibrium by the reservoirs leading to a net flow of solitons in the stationary state.

We define $p_{\bm{n}}^{t}$ to be the probability that $\bm{n}^{t}=\bm{n}$ and $\bm{p}^{t}=\sum_{n_{1},n_{2},\dots,n_{N}\in\{0,1\}}p_{\bm{n}}^{t}\ e_{n_{1}}\otimes e_{n_{2}}\otimes\dots\otimes e_{n_{N}}$ the associated probability vector in $(\mathbb{R}^{2})^{\otimes N}$ (where $e_{0}$ and $e_{1}$ denote the elementary basis of $\mathbb{R}^{2}$). The master equation can be written as $\bm{p}^{t+1}=M\bm{p}^{t}$ where the Markov matrix $M=M_{\rm o}M_{\rm e}$ is expressed as the product of two operators associated with the even and odd half-time steps, see Fig. 1(b),

[TABLE]

The subscripts indicate on which sites of the lattice the operators are acting non-trivially. The operator $P$ is the $8\times 8$ permutation matrix (acting on three sites) that enforces the dynamical rule in the bulk, with elements $P_{nn^{\prime}n^{\prime}\!{}^{\prime}}^{mm^{\prime}m^{\prime}\!{}^{\prime}}=\delta_{n,m}\delta_{\chi(nn^{\prime}n^{\prime}\!{}^{\prime}),m^{\prime}}\delta_{n^{\prime}\!{}^{\prime},m^{\prime}\!{}^{\prime}}$ where $\delta_{n,m}$ is the Kronecker symbol. The operators $L$ and $R$ are the $4\times 4$ stochastic matrices for the boundary processes of first and last site, respectively. A typical trajectory of the RCA54 is shown in Fig. 1(c). For further details of the model see Prosen and Mejía-Monasterio (2016); Prosen and Buča (2017); Inoue and Takesue (2018) and SM .

Large deviations of time-integrated observables.– We are interested in the statistics of general (possibly inhomogeneous) space- and time-extensive observables of the form

[TABLE]

in the large time $T$ limit. These are dynamical (or trajectory) observables as they depend on the full time history $(\bm{n}^{0},\bm{n}^{1/2},\bm{n}^{1},\dots,\bm{n}^{T-1/2})$. An example is the time-integrated number of up sites (which is not conserved in the RCA54) corresponding to $f_{j}(n,n^{\prime})=(n+n^{\prime})/2$ and $g_{j}(n,n^{\prime})=0$.

For large $T$ the probability of ${\mathcal{O}}_{T}$ has a LD form, $P_{T}({\mathcal{O}})=\langle\delta({\mathcal{O}}-{\mathcal{O}}_{T})\rangle\sim_{T\to\infty}e^{-T\varphi_{N}({\mathcal{O}}_{T}/T)}$, where $\varphi_{N}(x)$ is the rate function (where the subscript indicates its size-dependence). The moment generating function also has a LD form, $Z_{T}(s)=\langle e^{-s{\mathcal{O}}_{T}}\rangle\sim e^{T\theta_{N}(s)}$, where $\theta_{N}(s)$ is called the scaled cumulant generating function (SCGF), as its derivatives at $s=0$ correspond to the cumulants of ${\mathcal{O}}_{T}$ divided by time. The LD functions play the role of free-energies for trajectories and are related by a Legendre transform, $\theta_{N}(s)=-\min_{x}\left[sx+\varphi_{N}(x)\right]$.

To obtain the SCGF we deform, or tilt, the Markov matrix Touchette (2009): we define $M(s)=M_{\rm o}\,G(s)\,M_{\rm e}\,F(s)$, where we have introduced the diagonal operators $F_{\bm{n},\bm{n^{\prime}}}(s)=\delta_{\bm{n},\bm{n^{\prime}}}\prod_{i=1}^{N-1}f_{n_{i},n_{i+1}}^{(i)}$ and $G_{\bm{n},\bm{n^{\prime}}}(s)=\delta_{\bm{n},\bm{n^{\prime}}}\prod_{i=1}^{N-1}g_{n_{i},n_{i+1}}^{(i)}$, with the shorthand notation $f_{n,n^{\prime}}^{(i)}=e^{-sf_{i}(n,n^{\prime})}$ and $g_{n,n^{\prime}}^{(i)}=e^{-sg_{i}(n,n^{\prime})}$. We then have that $\theta_{N}(s)=\ln\lambda(s)$, where $\lambda(s)$ is the largest real eigenvalue of $M(s)$.

Exact results from Matrix Ansatz.– During the last decades, a technique called Matrix Ansatz has proven to be very efficient for deriving exact results in out of equilibrium systems. It has been introduced to compute analytically the stationary state of Markov chains Derrida et al. (1993) (see also Vanicat (2017, 2018) for recent developments) and has later on been used to compute deformed ground state Gorissen et al. (2012); Finn and Vanicat (2017), eigenvectors Crampé et al. (2011); Prosen and Buča (2017) or time-evolution of particular initial states Klobas et al. (2018). We use it here to compute the ground-state of the tilted Markov operator $M(s)$. Following the approach of Prosen and Mejía-Monasterio (2016); Prosen and Buča (2017), our strategy is to look for vectors $\bm{p}$ and $\bm{p^{\prime}}$ such that $M_{\rm e}F(s)\bm{p}=\lambda_{R}(s)\bm{p^{\prime}}$ and $M_{\rm o}G(s)\bm{p^{\prime}}=\lambda_{L}(s)\bm{p}$. It then follows that $\lambda(s)=\lambda_{R}(s)\lambda_{L}(s)$ is the dominant eigenvalue of the tilted Markov operator, $M(s)\bm{p}=\lambda(s)\bm{p}$.

It turns out that one can construct four pairs of site-dependent $3\times 3$ matrices $W_{n}^{(j)}$, $V_{n}^{(j)}$, $X_{n}^{(j)}$, $Y_{n}^{(j)}$ SM which satisfy the inhomogeneous bulk relations, for $j$ even:

[TABLE]

as well as six row 3-vectors $\bra{l_{n}},\bra{l^{\prime}_{nn^{\prime}}}$ and six column 3-vectors $\ket{r_{nn^{\prime}}}$, $\ket{r^{\prime}_{n}}$, satisfying the boundary equations

[TABLE]

These equations provide a cancellation scheme implying that an eigenvector of $M$, specifically vectors $\bm{p}$ and $\bm{p}^{\prime}$, take the matrix product form

[TABLE]

An explicit expression of the matrices and boundary vectors in the 3-dimensional auxilliary space are provided in SM . The eigenvalue $\lambda=\lambda_{\rm L}\lambda_{\rm R}$ is proven to be the dominant root of a polynomial of order 4

[TABLE]

with

[TABLE]

and where

[TABLE]

and, $\alpha^{\prime}=\alpha+\widetilde{\alpha}$, $\beta^{\prime}=\beta+\widetilde{\beta}$, $\gamma^{\prime}=\gamma+\widetilde{\gamma}$, $\delta^{\prime}=\delta+\widetilde{\delta}$,

[TABLE]

When $s=0$, that is, in the non-deformed case, the polynomial factorizes as in Prosen and Buča (2017) and the largest eigenvalue becomes $\lambda=1$ as expected.

Dynamical phase transition.– From Eq. (4) we can obtain the behavior of the SCGF $\theta_{N}(s)$ in the large size limit. Since the observables we consider are extensive in system size, cf. Eq. (2), we have that $a:=-\lim_{N\to\infty}(\ln{\mathfrak{a}}_{N})/(Ns)$, $b:=-\lim_{N\to\infty}(\ln{\mathfrak{b}}_{N})/(Ns)$ and $c:=-\lim_{N\to\infty}(\ln{\mathfrak{c}}_{N})/(Ns)$ exist and are finite. The SCGF then takes the scaling form

[TABLE]

where the function ${\vartheta}(\sigma)$ is defined such that $\exp[{\vartheta}(\sigma)]$ is the largest real root of the polynomial

[TABLE]

The scaling form (6) provides us immediately with system size behavior of the long-time cumulants of ${\mathcal{O}}_{T}$

[TABLE]

where $\langle\!\langle\cdot\rangle\!\rangle$ indicates the cumulant. The supra-linear dependence on size for $k\geq 2$ indicates the presence of a singularity at $s=0$ in the large size limit.

We can extract explicitly from (Exact large deviation statistics and trajectory phase transition of a deterministic boundary driven cellular automaton) the exact asymptotic of the first few cumulants. From $k=1$ we get the average observable per unit time

[TABLE]

while from $k=2$ the corresponding susceptibility

[TABLE]

with $\mu=\gamma(1-\alpha)+\beta(1-\gamma)$ and $\nu=\delta(1-\alpha)+\alpha(1-\gamma)$.

The scaling function ${\vartheta}(\sigma)$ has the following properties: (i) at $\sigma=0$ it vanishes as the polynomial (Exact large deviation statistics and trajectory phase transition of a deterministic boundary driven cellular automaton) trivially factorizes; (ii) it is a convex function and ${\vartheta}^{\prime\prime}(\sigma)$ admits a global maximum $\sigma^{*}$; (iii) if $(b+c)>0$ it has the asymptotic behavior

[TABLE]

[if $(b+c)<0$ the asymptotic behavior is obtained by $\sigma\to-\sigma$]. We can thus deduce that the SCGF converges to the limit shape

[TABLE]

when $(b+c)>0$ [for $(b+c)<0$ the shape is obtained by $s\to-s$]. The singularity at $s=0$ corresponds to a first-order phase transition.

Figure 2 shows the SCGF for one choice of the observable (the time-integrated number of up sites). As $N$ grows, $\theta_{N}(s)$ approaches the piecewise linear form (11). The order parameter, $-\theta^{\prime}(s)/N$, changes from a large value for $s$ negative to one close to zero for $s$ positive, the change becoming discontinuous for $N\to\infty$. The increasing sharpness of the crossover is manifested in the behaviour of the susceptibility $\theta^{\prime\prime}(s)/N$, whose peak grows as $N$. The value of $s$ at its peak indicates the location of the finite size crossover, which goes as $s_{c}\propto N^{-1}$. The finite size scaling of the transition point is similar to that expected in the FA model Bodineau et al. (2012); Bodineau and Toninelli (2012); Nemoto et al. (2017), while the scaling of the susceptibility is different.

The Inset to Fig. 2 shows the rate function $\varphi(x)$ where $x={\mathcal{O}}_{T}/TN$ for various sizes, as obtained from the SCGF via the Legendre transform. As the size increases $\varphi(x)$ progressively broadens. For finite $N$ the broadening is indicative of large fluctuations, and a precursor of the phase transition. In the limit $N\to\infty$, it takes the shape of a flat square well, corresponding to the Maxwell construction due to the first-order coexistence of the two dynamical phases, the inactive one with $x_{\rm min}=a$ and the active one with $x_{\rm max}=\frac{1}{3}(b+c+3a)$, cf. Eq. (11) ($\varphi=\infty$ elsewhere). Due to the integrable nature of the RCA54 – and in contrast to facilitated models Garrahan et al. (2007, 2009); Garrahan (2018) – there are no fluctuations within each dynamical phase, which means that for finite $N$ the rate function should have the shape of a “tilted ellipse” Jordan and Sukhorukov (2004); Lambert et al. (2015); Brandner et al. (2017).

Doob transformation and optimal dynamics.– The dynamical phase transition above corresponds to a singular change at the level of fluctuations: if the ensemble of trajectories is reweighed by $e^{-s{\mathcal{O}}_{T}}$ (the so-called $s$-ensemble Hedges et al. (2009); Garrahan (2018)), there is a singular change in the nature of atypically active trajectories ($s<0$) to atypically inactive ones ($s>0$). These reweighed ensembles can be sampled from the original dynamics by post-processing, but this is exponentially costly in $T$. However, they can be optimally accessed in terms of an “auxiliary” Jack and Sollich (2010) or “driven” Chetrite and Touchette (2015a) Markov process, by means of a so-called generalised Doob transformation; see also Borkar et al. (2003); Chetrite and Touchette (2015b); Jack and Sollich (2015); Garrahan (2016); Carollo et al. (2018); Derrida and Sadhu (2018).

From the matrix product construction of the leading eigenvector of $M(s)$ we can obtain the exact long-time Doob operator Jack and Sollich (2010); Chetrite and Touchette (2015a); Garrahan (2016):

[TABLE]

where ${\cal L}$ is a diagonal operator formed out of components of the leading left eigenvector $\bm{q}$ of $M(s)$, i.e. ${\cal L}_{\bf{n},\bf{n^{\prime}}}:=\delta_{\bf{n},\bf{n^{\prime}}}q_{\bf{n}}$. The exact matrix product construction of $\bm{q}$ in given in SM . The operator (12) is a Markov matrix for a stochastic dynamics whose trajectories are guaranteed to coincide - for long-times - with those of the $s$-ensemble of the original $M$ Jack and Sollich (2010); Chetrite and Touchette (2015a, b); Jack and Sollich (2015); Garrahan (2016). Figure 3 shows how $M_{\rm Doob}$ allows to sample the fluctuations of the RCA54 parameterised by $s$ in an optimal way.

Note that due to the deterministic nature of the bulk dynamics of the RCA54, the Doob operator amounts to a non-trivial modification only of the boundary probabilities, which in $M_{\rm Doob}$ depend on the configuration of the whole lattice. For details see SM .

Conclusion.– We studied the statistics of a general class of dynamical observables in a cellular automaton with stochastic boundary driving. We provided an exact expression of their scaled cumulant generating functions by means of an inhomogeneous matrix product expression for the leading eigenvector of the corresponding tilted Markov operator. Our results give a precise analytical description of the phase transition between active and inactive dynamical phases observed in a wide range of other models.

We foresee extensions of our work here in several directions, including computing the large deviation statistics of currents, and even the complete “level 2.5” statistics for the empirical measure and fluxes Maes and Netocny (2008); Bertini et al. (2015); Barato and Chetrite (2015); Hoppenau et al. (2016); Bertini et al. (2018). The analytic inhomogeneous Matrix Ansatz introduced here could also be used to address similar questions in more complicated models, for instance for cellular automata with asymmetric constraints, and for systems with stochastic dynamics in the bulk.

Acknowledgements.– The work here has been supported by Advanced Grant 694544 – OMNES of the European Research Council (ERC), Grants No. N1-0025, N1–0055 of the Slovenian Research Agency (ARRS), EPSRC programme grant EP/P009565/1, EPSRC Grant No. EP/R04421X/1, Leverhulme Grant RPG-2018-181.

S1 Markovian dynamics

S1.1 Markov operator

We provide here details about the Markov matrix associated to the model Prosen and Mejía-Monasterio (2016); Prosen and Buča (2017). The matrix acts on the vector space $(\mathbb{R}^{2})^{\otimes N}$ constructed as the tensor product of $N$ copies of the elementary space $\mathbb{R}^{2}$ (each associated to one site of the lattice). It can be written as the product of two operators, corresponding to the two half-time steps of the stochastic dynamics

[TABLE]

with

[TABLE]

We have used the notation

[TABLE]

where $\text{1}_{2}$ is the identity matrix acting on the elementary space $\mathbb{R}^{2}$. $P$ is a $8\times 8$ permutation matrix encoding the deterministic dynamics in the bulk and acting on the space $(\mathbb{R}^{2})^{\otimes 3}$ (i.e. it acts on three sites of the lattice) $L$ and $R$ are $4\times 4$ matrices encoding the stochastic dynamics at the boundaries and acting on the space $(\mathbb{R}^{2})^{\otimes 2}$ (i.e. they act on two sites of the lattice). These local matrices read explicitly Prosen and Mejía-Monasterio (2016); Prosen and Buča (2017)

[TABLE]

and

[TABLE]

S1.2 Tilted operator

In order to study the statistics of the time-integrated observables we introduce a deformation of the Markov operator

[TABLE]

It will be convenient for the following, when deriving the cancellation scheme for the matrix product expression of the leading eigenvector, to introduce the following diagonal matrices

[TABLE]

Then the extensive tilt operators $F(s)$ and $G(s)$ can be expressed as

[TABLE]

where we have again used subscript index to denote the sites of the lattice on which the operators are acting non-trivially. A direct computation shows that we can distribute those local operators $F^{(j)}$ and $G^{(j)}$ on the local matrices encoding the dynamics of the model

[TABLE]

where

[TABLE]

and, for $i$ even,

[TABLE]

The structure of the tilted Markov operator given in (S9) as the product of local operators will be very convenient to prove the matrix ansatz construction of the leading eigenvector.

S2 Matrix product expression of the leading eigenvector

This section is devoted to the construction of the leading eigenvector of the tilted operator in a matrix product form. Recall that our strategy is to look for a pair of vectors $\bm{p}$ and $\bm{p^{\prime}}$ such that

[TABLE]

and such that $\lambda(s)=\lambda_{R}\lambda_{L}$ is the dominant eigenvalue of the tilted Markov operator $M(s)$. This would indeed imply $M(s)\bm{p}=\lambda(s)\bm{p}$.

We are going to explain the exact construction of these vectors in a matrix product form

[TABLE]

where $W_{n}^{(k)},V_{n}^{(k)}$, $n\in\{0,1\}$, are $3\times 3$ matrices acting in an auxiliary space, and where $\bra{l_{n}},\bra{l^{\prime}_{n,n^{\prime}}}$, $n,n^{\prime}\in\{0,1\}$, are row vectors in the auxiliary space (with $3$ entries) and $\ket{r^{\prime}_{n}},\ket{r_{n,n^{\prime}}}$, $n,n^{\prime}\in\{0,1\}$, are column vectors in the auxiliary space (with $3$ entries). The proof of the matrix product structure of this pair of vectors will be done in two steps. First we will present a cancellation scheme, which shows that the pair of relations (S12) can be obtained if the auxiliary matrices and vectors arising in the ansatz (S13) satisfy very simple algebraic relations. Second we will construct explicit solutions of these algebraic relations using $3\times 3$ matrices and $3$-components vectors.

S2.1 Cancellation scheme

In order to prove efficiently our exact construction and to present the cancellation scheme, we introduce the vectors with matrix components

[TABLE]

with $X_{n}^{(j)}$, $n\in\{0,1\}$, $i\in\{2,\ldots,N-1\}$ being $3\times 3$ matrices acting in the auxiliary space, which will play an important role in the proof of the matrix ansatz. We also introduce

[TABLE]

Using those new objects we can rewrite concisely the pair of vectors

[TABLE]

We used once again tensor product notations: the subscript indices stand for the tensor product component (i.e. the site of the lattice) the vectors are belonging to. As an illustration of our convention, we have for instance

[TABLE]

The tensor product notation (S16) will be very convenient to present the cancellation mechanism.

The cancellation scheme consists of two key relations. First, the inhomogeneous bulk relations, written for even site index $j=2i$:

[TABLE]

Second, the boundary equations:

[TABLE]

Using tensor product notations, the inhomogeneous bulk equation can be concisely rewritten, for even $j$:

[TABLE]

Similarly, the boundary equations can be rewritten, as

[TABLE]

for the left boundary, and

[TABLE]

for the right boundary. A direct computation shows that the algebraic relations (S20), (S21) and (S22) imply the eigenvalue equations (S12). We illustrate the mechanism in the case of a lattice of length $N=6$

[TABLE]

The relation $M_{\rm o}G(s)\bm{p^{\prime}}=\lambda_{L}\bm{p}$ is established in a very similar way. Note that the value of the parameters $\lambda_{L}$ and $\lambda_{R}$ are still free at this stage. They will be fixed when solving explicitly the algebraic relations (S20), (S21) and (S22) with finite dimensional matrices and vectors, and will provide the eigenvalue $\lambda(s)=\lambda_{L}\lambda_{R}$.

S2.2 Solution to the inhomogeneous bulk algebra

We proceed to show how to explicitly solve the bulk relations (S20). We start by the matrix ansatz

[TABLE]

[TABLE]

[TABLE]

We stress here that these matrices are acting on a 3-dimensional auxiliary space which is different from the physical space of local configuration $\mathbb{R}^{2}$. This 3-dimensional vector space is traced out in the matrix product expression (S13) thanks to the action of the boundary vectors.

Plugging the matrix ansatz into (S20), we end up with the following two-parameter $(\rho,\kappa)$ solution

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

where

[TABLE]

with the convention that $f_{nn^{\prime}}^{(-1)}=f_{nn^{\prime}}^{(0)}=1$ and $g_{nn^{\prime}}^{(-1)}=g_{nn^{\prime}}^{(0)}=1$.

Note that the 3-dimensional auxiliary space, arising in the explicit representation of the inhomogeneous bulk algebra, plays essentially the role of detecting the left and right moving solitons in a given configuration. More precisely, in order to detect if a soliton is located at a given site of the lattice, it is necessary to know the content of the two neighboring sites (the number of left and right moving solitons is indeed a locally conserved charge with support of size three). When scanning a configuration from left to right, the auxiliary space allows to store the information about the content of the two last sites and thus permits the detection of the solitons. The $3\times 3$ matrices associate a weight (that could be thought as a momentum) to each soliton, which depends on its position on the lattice and on its direction of propagation. We will investigate further in future work the connection between the conserved charges of the model and the matrix product construction of its eigenvectors.

S2.3 Solution to the boundary relations

We first recall the convenient notations $\alpha^{\prime}=\alpha+\widetilde{\alpha}$, $\beta^{\prime}=\beta+\widetilde{\beta}$, $\gamma^{\prime}=\gamma+\widetilde{\gamma}$, $\delta^{\prime}=\delta+\widetilde{\delta}$,

[TABLE]

It will be also useful to introduce the notations

[TABLE]

with the convention that $f_{nn^{\prime}}^{(N)}=f_{nn^{\prime}}^{(N+1)}=1$ and $g_{nn^{\prime}}^{(N)}=g_{nn^{\prime}}^{(N+1)}=1$.

The resolution of the left boundary relation (S21) fixes the value of the parameters $\rho$ and $\kappa$ appearing in the matrices $W^{(i)}_{n}$ and $V^{(i)}_{n}$

[TABLE]

and imposes the following explicit expression for the left boundary vectors

[TABLE]

where we have introduced the parameter

[TABLE]

The resolution of the right boundary relation fixes

[TABLE]

and imposes the following explicit expression for the right boundary vectors

[TABLE]

S2.4 Leading eigenvalue

Imposing the equality between the two constraints (S34) and (S37), we obtain that the eigenvalue $\lambda(s)=\lambda_{R}\lambda_{L}$ has to be the dominant root of a polynomial of order $4$

[TABLE]

with

[TABLE]

Note that without counting fields, i.e. for $s=0$, we have the factorization

[TABLE]

and we recover that the dominant eigenvalue is $\lambda=1$.

S3 Doob transformation

S3.1 Explicit form of the Doob transform

The long-time Doob transform can be defined as Carollo et al. (2017),

[TABLE]

where $\lambda(s)$ is the leading eigenvalue of the tilted propagator $M(s)$ and ${\cal L}$ is a diagonal operator formed out of components of the leading left eigenvector $\bm{q}$ of $M(s)$,

[TABLE]

where ${\bf{n}}:=(n_{1},n_{2},\dots,n_{N})$, $n_{k}\in\{0,1\}$. As mentioned in the main text, the final result of this transform and the deformation is only a change to the stochastic boundary matrices $R$ and $L$ that makes the transition rates dependent on configuration of the ingoing and outgoing state in the probability vector. The bulk dynamics remain deterministic. The explicit expression can be written as,

[TABLE]

[TABLE]

where we defined,

[TABLE]

where we used that $[L_{12},P_{234}]=0$

S3.2 Leading left eigenvector

To calculate the Doob transformation we need the left leading eigenvector $\bm{q}M(s)=\lambda(s)\bm{q}$. In order to reduce the computation of the left eigenvector to a similar problem as the computation of the right eigenvector, and thus use partially our previous result, we define,

[TABLE]

It is clear that the leading right eigenvector of $N(s)$, $N(s)\bm{o}=\lambda(s)\bm{o}$ is related to $\bm{q}$ via a similarity transformation, $\bm{q}=(F(s)\bm{o})^{\intercal}$. We also have,

[TABLE]

Note that $N(s)$ has formally a very similar structure than $M(s)$. The differences are that the odd and even half-time steps are switched, $F(s)$ and $G(s)$ are also interchanged and $R_{N-1,N}$ (respectively $L_{12}$) is replaced by $\left(R_{N-1,N}\right)^{\intercal}$ (respectively $\left(L_{12}\right)^{\intercal}$). As previously, we will construct a pair of vector $\bm{o}$ and $\bm{o^{\prime}}$ satisfying the relations $M_{\rm o}^{\intercal}F(s)\bm{o}=\lambda_{R}\bm{o^{\prime}}$ and $M_{\rm e}^{\intercal}G(s)\bm{o^{\prime}}=\lambda_{L}\bm{o}$. Taking $\bm{o}$ (respectively $\bm{o^{\prime}}$) to be of the same matrix product form as $\bm{p^{\prime}}$ (respectively $\bm{p}$),

[TABLE]

we arrive to a similar set of equations for the cancellation scheme as we did for $\bm{p^{\prime}}$ and $\bm{p}$.

The inhomogeneous bulk relations now read (for even site index $j=2i$):

[TABLE]

They are similar to (S18) but with $f^{(k)}_{nn^{\prime}}$ and $g^{(k)}_{nn^{\prime}}$ interchanged. An explicit representation of these inhomogeneous bulk relations is thus obtained using our previous results

[TABLE]

The boundary equations read

[TABLE]

They are similar to (S19) but with $f^{(k)}_{nn^{\prime}}$ and $g^{(k)}_{nn^{\prime}}$ interchanged and with the matrix $R$ (respectively $L$) replaced with $R^{\intercal}$ (respectively $L^{\intercal}$). They can be solved explicitly. The solution for the right boundary equations gives,

[TABLE]

Likewise, we have for the right boundary vectors in the MPS,

[TABLE]

The resolution of the left boundary relation fixes the value of the parameters $\rho$ and $\kappa$

[TABLE]

and imposes the following explicit expression for the left boundary vectors

[TABLE]

where we have introduced the parameter

[TABLE]

Curiously, the solution for the left eigenvectors $\bm{o}$ is quite similar to the one for the right eigenvectors. We note that we could have solved for the left eigenvectors alternatively. By first multiplying the first bulk exchange relation (S20) by $\left(f^{(2k-1)}_{nn^{\prime}}f^{(2k)}_{n^{\prime}n^{\prime\prime}}\right)^{-1}P_{nn^{\prime}n^{\prime\prime}}$ and the second one in (S20) by $\left(g^{(2k-2)}_{nn^{\prime}}g^{(2k-1)}_{n^{\prime}n^{\prime\prime}}\right)^{-1}P_{nn^{\prime}n^{\prime\prime}}$ and then replacing $g^{(k)}_{nn^{\prime}}\to 1/f^{(k)}_{nn^{\prime}}$ and $f^{(k)}_{nn^{\prime}}\to 1/g^{(k)}_{nn^{\prime}}$, we arrive at new bulk exchange relations. These new exchange relations are suitable for taking the same ansatz (S50) and solving $\bm{q}M_{\rm o}G(s)=\lambda_{R}\bm{q^{\prime}}$ and $\bm{q^{\prime}}M_{\rm e}F(s)=\lambda_{L}\bm{q}$ directly. We also would arrive to the same set of boundary equations that may be solved. We choose the above because we find that the similarity of the form of the explicit solution to the one for the right eigenvectors is more pronounced.