Optimal sampling of dynamical large deviations in two dimensions via tensor networks

We use projected entangled-pair states (PEPS) to calculate the large deviation statistics of the dynamical activity of the two dimensional East model, and the two dimensional symmetric simple exclusion process (SSEP) with open boundaries, in lattices of up to 40 × 40 sites. We show that at long times both models have phase transitions between active and inactive dynamical phases. For the 2D East model we find that this trajectory transition is of the first-order, while for the SSEP we find indications of a second order transition. We then show how the PEPS can be used to implement a trajectory sampling scheme capable of directly accessing rare trajectories. We also discuss how the methods described here can be extended to study rare events at finite times.

Introduction.-Tensornetwork (TN) techniques, whereas most actively developed in the context of quantum many-body physics [1][2][3][4][5][6], offer powerful numerical tools for much more general problems.They are based on an efficient parametrization of the many-body state in terms of local tensors (multidimensional arrays) connected according to a graph that, in general, responds to the structure of correlations in the state.In the last few years we have seen progress in their application to compute statistical properties of dynamical trajectories in classical stochastic systems.The first application was to the long time statistics-the dynamical large deviation (LD) regime-of one-dimensional lattice systems using variational algorithms (such as density matrix renormalization group [7], or DMRG) to approximate the leading eigenvectors of tilted Markov generators by matrix product states (or MPS, see e.g.Ref. [2]) [8][9][10][11][12][13][14][15].Building on these results, we recently introduced (i) a method which exploited MPS to efficiently sample long-time rare trajectories, and (ii) an MPS time-evolution to precisely compute trajectory statistics at finite times [16].
The suitability of the TN ansatz for these problems is rooted in the fact that the targeted eigenvectors have low correlations.More concretely, for local problems, we expect them to fulfill (up to small corrections) a socalled area law [17], according to which the "entanglement" (or a mathematically analogous quantity for classical probability vectors [18]) with respect to a bipartition is upper-bounded by the size of its boundary.This scaling is captured by the MPS family in one spatial dimension.In higher dimensions, a suitable generalization with area law is provided by the projected-entangled pair states (PEPS) [19], which were recently applied to the classical asymmetric exclusion process in two dimensions in Ref. [20].A computationally cheaper alternative, without an area law but accommodating more entanglement than MPS, is that of tree tensor networks [21], used for example in Refs.[22,23] (in combination with a time-dependent variational principle [24]) to study driven problems.
Here we use PEPS to study the LDs of the dynamical activity in two paradigmatic two-dimensional models, the 2D East model (also known as North-or-East model) [25][26][27][28], and the 2D symmetric simple exclusion process (SSEP) with open boundaries where particles can be injected and removed [29].We are able to accurately estimate the leading eigenvector of the tilted generator -and thus the LDs -of these models, and construct a closeto-optimal dynamics to directly sample the corresponding rare trajectories.Such an algorithm requires efficient sampling from the PEPS, and we show how to do this in the context of trajectory sampling.We benchmark our methods, showing how the PEPS allows for a controlled accuracy of optimal dynamics.We demonstrate that both models have a phase transition between active and inactive dynamical phases, a first-order transition for the 2D East and a second-order transition for the 2D SSEP. Models.
-The models we study here live in a twodimensional square lattice of size N = L × L, with each site being occupied by a binary variable n k = 0 or 1, where k = (k x , k y ) denotes the position of the site for k x , k y = 1 • • • L. Their continuous-time dynamics is defined by a Markov generator (e.g.see Refs.[31,32]), where |x⟩ and |y⟩ are configurations on the lattice, w x→y the transition rate from x to y, and R x = y̸ =x w x→y the escape rate out of x.We can write this as W = K−R, where K contains the off-diagonal transition rates, and R the diagonal escape rates.
The first model we consider is the 2D East model [25][26][27][28], often studied in the context of the glass transition.This is a kinetically constrained model such that an excited site n k = 1, allows ("facilitates") a site to its North or East to flip stochastically, see Fig. 1 is used for the contraction part of the network, with the dimension of its virtual bonds (blue legs) controlling accuracy; see Ref. [30] for details.
transitions rates: 0 → 1 with rate c, and 1 → 0 with rate 1 − c, subject to the kinetic constraint.In addition, we choose open boundary conditions with n (1,1) = 1 fixed.This ensures the entire state space remains dynamically connected [26].The second model is the 2D SSEP.This describes particles hopping to neighbouring sites on a 2D lattice with unit rate, but only if the target site is not already occupied by a particle.We also allow particles to be injected or removed at the boundaries of the lattice with rate 1/2, see Fig. 1(b).Exact definitions of the models are given in Ref. [30].Dynamical large deviations.-Weconsider the statistics of a dynamical observable K through its probability distribution where ω t denotes a stochastic trajectory and π(ω t ) its probability.The corresponding moment generating function is Z t (s) = ωt π(ω t )e −s K(ωt) .In the t → ∞ limit, the two obey LD principles P t (K) ≍ e −tφ(K/t) and Z t (s) ≍ e tθ(s) , with the rate function φ(K/t) and scaled cumulant generating function (SCGF) θ(s) related through a Legendre transform, θ(s) = − min k [sk + φ(k)] for k = K/t (for reviews see Refs.[32][33][34][35].) We focus as an observable on the dynamical activity [36,37], which counts the number of jumps in a trajectory.The relevant operator to study is the tilted generator [32][33][34][35], which for the activity reads W s = e −s K − R, with the LD statistics encoded in the leading eigenvalue and (right and left) eigenvector(s), W s |r s ⟩ = θ(s) |r s ⟩ and ⟨l s | W s = θ(s) ⟨l s |.
Projected-entangled pair states.-TNmethods allow us to solve the problem above using a variational ansatz for |ψ s ⟩ in the PEPS family, a natural generalization of MPS for area law states and lattices in more than one spatial dimension [19].PEPS parametrize the many-body state with one rank-5 tensor per lattice site, in which the physical index has the dimension of the site variable (0, 1 in our case), and the remaining four virtual indices each have bond dimension D PEPS , determining the maximum number of parameters in the ansatz, N p = N dD 4  PEPS (see Fig. 1 for a graphical representation).Several TN algorithms exist to optimize the PEPS by maximizing the expectation of a local stochastic generator.Crucial to all of them is an efficient computation of expectation values for local operators, such as the terms in W s .We use the boundary MPS scheme [1,19] (illustrated in Fig. 1), which approximates the partial contraction of the network by a MPS, whose bond dimension χ B controls the accuracy of the contraction.A heuristic choice for local problems is χ B ∼ O(D 2 PEPS ) (see e.g.Ref. [38]).
In order to find the PEPS approximation to the leading eigenvector |ψ s ⟩ we employ time evolution, which effectively projects the ansatz onto the leading eigenvector by iteratively applying short evolution steps, decomposed in two-body terms that are applied sequentially [6,39].After each operation, the directly affected pair of tensors is updated.This requires a strategy to approximate their environment (i.e. the contraction of the remaining TN).After comparing to quasi-exact [40] results and to a more expensive strategy, we find that the computationally cheapest simple update (SU) [41], with maximal PEPS bond dimension D PEPS = 4, is enough to achieve sufficiently accurate measurements of the SCGF in our problems, and allows us reaching large sizes at low computational cost, scaling only as O(D 5 PEPS ).For further details on the numerical approach we use see Ref. [30].
Large deviations from PEPS.-The East and SSEP in 1D are known to have dynamical phase transitions in terms of the activity or other dynamical observables [11,13,36,[42][43][44][45][46][47][48][49].In two-dimensions, the SSEP has a transition in the LDs of the current [20].We now provide evidence by means of PEPS for both the 2D East and 2D SSEP having active-inactive phase transitions.Figure 2(a-c) shows the LD statistics for both the 2D East model (top) and the 2D SSEP (bottom).For the East model, we see from Fig. 2(a) that the SCGF follows a linear response, θ(s) ≈ sk(0), for small s, but at s c (L) it sharply changes to another branch.This point corresponds to a sudden drop in activity, k(s) = −θ ′ (s), which becomes discontinuous in the limit N → ∞, see Fig. 2(b).Having access to both the SCGF and the dynamical activity allows us to estimate the rate function φ(k), shown in Fig. 2(c).We see a broadening of the rate function around the mean, indicating the coexistence of active and inactive dynamics.For comparison, we also show the distribution of a simple process with the same mean activity, but which is uncorrelated in time (black dashed line).All this behaviour is characteristic of a first-order phase transition.
For the SSEP we see something different: Fig. 2(a) shows no sharp change in θ(s), and the activity in Fig. 2(b) has no discontinuity.This is indicative of a second-order transition, with the rate function showing critical broadening, see Fig. 2(c), and a divergence in the susceptibility χ(s) = θ ′′ (s), see Fig. 2(d).For both models we can extract a transition point from the drop in either first or second cumulant.The top panel of Fig. 2(d) shows how the transition point scales with L for both models (for a range of c for the 2D East).We are able to fit the data with the power laws s c (L) ∼ L −2α , as shown by the solid lines.We find the exponents α ≳ 1 for the 2D East and α ≲ 1 for the SSEP, see the inset to the top panel of Fig. 2(d).
Optimal sampling of rare trajectories from PEPS.-Sampling trajectories corresponding to the s ̸ = 0 phases is difficult as they are exponentially rare in system size and time.The optimal sampling dynam-ics at long times is given by the so-called generalized Doob transform [51][52][53][54][55], which maps the tilted generator into a true stochastic generator for the rare trajectories, where L is the leading left eigenvector of W s as a diagonal matrix.This gives a new dynamics with the transition rates wx→y = l s (y) with l s (x) = ⟨l s |x⟩.In W Doob s the counting field s appears as a physical control parameter, and running dynamics with rates (S15) gives trajectories at s ̸ = 0 on demand.While optimal, W Doob s is difficult to construct in general as one needs the exact left leading eigenvector.However, we can exploit our PEPS approximation to estimate the rates Eq.(S15), similar to Ref. [56] for 1D and MPS, see Ref. [30] for more information.
To obtain Eq. (S15) for the transitions out of a state x we calculate l s (y) from the PEPS using a boundary dimension χ B = D PEPS [57][58][59][60], thus entailing a maximum cost O(N D 6 PEPS ).If we neglect the time edges of trajectories, we can estimate an time-extensive observable by The snapshots on the right show the configurations at the marked times (black/white indicates a particle/hole).See also Refs.[30,50].
importance sampling where α t denotes a trajectory generated with Eq. (S15) (the reference dynamics), and O(α t ) is the trajectory observable.The re-weighting factor g(α t ) is where R(t ′ ) and R(t ′ ) are the escape rates of the system at time t ′ in the original dynamics and the approximate Doob dynamics, respectively.Notice that with a large enough number of trajectories, Eq. ( 3) can be used to correct on the imperfections in the reference dynamics due to an imperfect PEPS approximation.Figures 3 show results from our sampling algorithm for the 2D East with c = 0.5 and the 2D SSEP, both for system sizes N = 22 × 22.The average dynamical activity measured in trajectories (symbols) [with umbrella sampling Eqs.(3,4)] coincides with that obtained directly from the PEPS (solid line), except for D PEPS = 1 for the East model.The accuracy of our dynamics is quantified by the variance of the time integrated difference in escape rates, cf.Eq. ( 4), which vanishes for the exact Doob rates.We show this for each D in the insets of Figs.3: increasing the D PEPS consistently reduces the variance, indicating a better sampling dynamics and less need for importance sampling.
The PEPS approximation to the leading eigenvalue gives us direct access to the long time-averaged properties of the dynamics.However, the broader range of dynamical information-such as temporal correlations-can only be obtained through the simulation of the dynamics in rare dynamical regimes.The ability to exploiting the PEPS to define an efficient trajectory sampling scheme for the biased dynamics allows us to characterize the finite-time atypical behaviour beyond what is directly encoded in the PEPS approximation.Discussion.-We have shown that the dynamical LDs of two-dimensional stochastic models can be studied efficiently with PEPS, including the quasi-optimal sampling of atypical trajectories.Compared to more standard methods [61][62][63][64][65][66][67], PEPS offer the advantage of a computational cost that scales only polynomially in the system size and the tensor dimensions D PEPS .Furthermore, the algorithms produce an explicit ansatz for the leading eigenvector encoding the LDs, which allows to extract all (local) observables in the biased trajectory ensemble and, as shown above, to devise a near-to-optimal sampling dynamics.Additionally, because PEPS form a complete family, by increasing the bond dimension the quality of the solution can be systematically improved, a property that can also be utilized to estimate the error of the approximation.Our results show that the PEPS ansatz is well suited for these problems, as a moderate bond dimension suffices to explore large systems.
We showed here that both the 2D East model and the 2D SSEP have active-inactive trajectory transitions, of the first-order and second-order, respectively.Our work adds to the continuously expanding application [8-16, 20, 22, 23, 68] of tensor network methods to study the dynamical fluctuations in classical stochastic systems.
There are several interesting avenues to pursue building on this work.One is to integrate 2D trajectory sampling via TNs with a method such as transition path sampling [69] for investigating statistics of fluctuations at finite times, cf.Ref [16,70].While the current implementations with PEPS are too demanding to reasonably incorporate transition path sampling, tree tensor networks [21] are a promising alternative that could allow to reliably investigate finite time scaling.We hope to report on this is the near future.
Acknowledgements.-We acknowledge financial support from EPSRC Grant no.EP/R04421X/1 and the Leverhulme Trust Grant No. RPG-2018-181.M.C.B. acknowledges support from Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2111 -390814868.Calculations were performed using the Sulis Tier 2 HPC platform hosted by the Scientific Computing Research Technology Platform at the University of Warwick.Sulis is funded by EPSRC Grant EP/T022108/1 and the HPC Midlands+ consortium.We acknowledge access to the University of Nottingham Augusta HPC service.Example code used to generate the data can be found at https: //github.com/lcauser/2d-optimal-sampling.Data is available at Ref. [50].
We consider systems that live on a two-dimensional square lattice of size N = L × L, where each site can take the binary values n k = 0 or 1, and k = (k x , k y ) denotes the position of the lattice sites, k x , k y = 1 • • • L. The system evolves under continuous-time Markov dynamics, defined by the transition rates w x→y from configuration x to y.The average dynamics of the system can be encoded by a probability distribution P x (t), which describes the probability for the system to be in some configuration x at time t.This can be compactly written as a vector of probabilities, |P (t)⟩ = x P x (t) |x⟩, where x P x (t) = 1.The evolution of the probability distribution is given by the master equation, For the models considered here, the Markov generator can easily be written in terms of local operators.The first model we consider is the 2D East model, with the Markov generator where σ ± k are the Pauli raising/lowering operators at site k, c ∈ (0, 1/2] parameterizes the transition rates, and the kinetic constraint is P (kx,ky) = n (kx−1, ky) + n (kx, ky−1) .The second model we consider is the 2D symmetric simple exclusion process (SSEP), where ⟨k, l⟩ denotes a pair of nearest neighbours, σ x i = σ + i + σ − i , and ∂ denotes the boundary of the lattice.See the main text for a description of each.

Mapping the tilted generator onto a Hermitian operator
The long time statistics for the models are encoded in the leading eigenvalue θ(s) and eigenvectors of the tilted Markov generator, W s , retrieved by multiplying the off-diagonal components of Eqs.(S3-S4) by e −s , The probability distribution over the configuration space under these biased dynamics then behaves as for sufficiently large times.The tilted Markov generator provides an efficient way to estimate the time-averaged dynamical properties by directly targeting its leading eigenvectors.However, we can exploit the fact that each of the models considered here obeys detailed balance.We define P as the diagonal matrix whose elements are the square roots of the steady state probability of W. Then we can use the similarity transformation H s = P −1 W s P, to define the Hermitian matrix H s , with maximal eigenvalue and associated eigenvector where |ψ s ⟩ is related to the original eigenvectors by |r s ⟩ = P |ψ s ⟩ and ⟨l s | = ⟨ψ s | P −1 .The models considered here allow for a simple representation as a Hermitian matrix, This representation is convenient due to the fact its expectation with any vector |ψ⟩ is bounded by the maximal eigenvalue through the Rayleigh-Ritz variational principle, Furthermore, using this Hermitian matrix means that we only need to determine one eigenvector.Notice that in this representation, the leading eigenvector encodes the probability amplitudes of each configuration, |ψ s ⟩ = x ψ(x) |x⟩, with |ψ(x)| 2 = P (x).

Projected entangled-pair states
The long-time dynamics of deformed Markov generators can be encoded by the probability amplitude, ψ(n), with configurations n = (n k1 , n k2 , • • • , n k N ).This can be written more compactly as a vector of probability amplitudes, |ψ⟩ = n ψ(n) |n⟩.This vector has the size 2 N for N lattice sites, and quickly becomes intractable to store.However, it is often the case that the components ψ(n) are not independent of each other, such that it is possible to efficiently approximate the vector by an object with a smaller dimensionality (described by a number of parameters less than d N ).This realisation is at the heart of tensor network (TN) approximations.A TN representation of ψ(n) amounts to decomposing the N -dimensional object into a network of many smaller tensors, connected along additional virtual dimensions, which are contracted (i.e.multiplied and summed over) to retrieve the original global tensor.
For the case of a 2D square lattice with N = L × L sites, the most natural TN ansatz is the projected entangled-pair state (PEPS).In this ansatz, each system (lattice site) is assigned its own rank-5 tensor, A dj kj , where k j denotes the position of system j, and d j the state of the system.One of the dimensions corresponds to the physical dimension of the subsystem with size d, and the other four dimensions are virtual ones which connect the tensor to the tensors of the four neighbouring lattice sites.These virtual dimensions are of size D PEPS (often referred to as the bond dimension), which controls the amount of mutual information shared between the lattice sites.It follows that each tensor is parameterized by dD 4  PEPS parameters, with a maximum of N p = N dD 4 PEPS parameters for the whole PEPS.By specifying the local configuration of each site, n k , and contracting over the network, one is able to determine ψ(n), where F is a function which represents the contraction over all virtual bond dimensions.It is convenient to represent the TN pictorially, as illustrated for |ψ⟩ in Fig. S1(a).Notice that we refer to this as |ψ⟩, as none of the physical dimensions are specified.The green cubes correspond to the tensors, the grey lines represent the virtual bond dimensions to be contracted over, and the pink lines represent the physical dimensions.Similarly, we show the conjugate, ⟨ψ|, in Fig. S1(b).It is also possible to represent local operators as a tensor, such as the two-body operators in the tilted generators.Local two-body operators can be represented by a rank-4 tensor, as shown in Fig. S1(c).Figure S1(d) shows how the probability amplitude ψ(n) can be retrieved from the PEPS by specifying the local dimension for each system.This reduces the PEPS to a network of rank-4 tensors, which give ψ(n) when contracted.

Contracting PEPS
PEPS allow for an efficient representation of probability amplitude vectors for large system sizes.However, to perform any tractable calculations, such as calculating the expectation of observables, it is necessary to have a way to efficiently contract the networks.Figures S2(a, b) show the networks which must be contracted to determine the inner products ⟨ψ|ψ⟩ and ⟨ψ| Ô|ψ⟩ respectively, for some two-body operator Ô.In general, contracting exactly such a two-dimensional TN is an intractable problem [71], as any exact contraction strategy has a cost that scales at least exponentially in L. For the ansatz to be of practical use, one needs ways to approximately (but precisely) contract the TN with a tractable cost.
A popular approach, and the one we will take here, is the boundary matrix product state (MPS) scheme [1,19].This involves contracting from the edge of the TN, one row (column) at a time, and approximating the result by an MPS (in this case, a collection of rank-4 tensors), whose bond dimension χ B controls the accuracy of the approximation.We demonstrate this process in Figs.S3(a-d).The tensors at the boundary can be contracted to form the first boundary MPS, see Fig. S3(a).This boundary MPS is then contracted with the subsequent layer, and approximated by another boundary MPS, again with bond dimension χ B , see Fig. S3(a).Applying the procedure from two opposing edges of the network up to the rows (or columns) of lattice sites neighbouring the two-body operator, one is able to approximate the expectation value depicted in Fig. S2(b) by the one shown in Fig. S3(e).This network can then be contracted from the edges exactly, resulting in Fig. S3(f), which can be easily and exactly contracted to determine the approximation of ⟨ψ| Ô|ψ⟩.Contracting the networks in Fig. S3 scales as ), in which case gives the total complexity O(D 10 ). Figure S4 demonstrates the role of the bond dimension of the boundary MPS, χ B .We optimize the MPS using the simple update described below, and calculate its expectation value with respect to H s , E(χ B ) = ⟨ψ s |H s |ψ s ⟩, measured with the bond dimension χ B .We compare the measurement to that using χ = 128, which can be considered to be quasi-exact, Indeed, the results confirm that χ B ∼ O(D 2 ) is a reasonable choice for approximating the environment.We use the bond dimension χ = 50 when checking for convergence in observables during the optimization procedure, and χ = 200 when taking the final measurements to ensure convergence in the boundary approximation.

Time evolution
To find a PEPS approximation for the leading eigenvector of the tilted generator, one needs to find a suitable optimization procedure.There are many approaches to optimizing TNs, but the one taken here will be to employ time evolution (often referred to imaginary time evolution for quantum many-body systems).The main idea is to project some probability amplitude vector onto the leading eigenvector of the tilted generator H s by applying the time propagator operator, U (δ) = e δHs , to |ψ⟩ until convergence is met, |ψ s ⟩ = lim t→∞ U (t) |ψ⟩ (up to normalization).In practice, the complete time propagation operator is difficult to compute, as it requires the matrix exponentiation of the tilted generator on the complete state space.However, for small δ ≪ 1, we can approximate the time propagation by a sequence of Trotter gates, where U k,l (δ) = e δH k,l , and H k,l are the terms in H s which act only on the lattice sites k and l.This approximation is often referred to as a first-order Trotter decomposition, with each set of gates having an error of O(δ 2 ).This is the approach we take.However, it is possible to improve the accuracy by using higher order Trotter decompositions [74].Each gate can be individually applied to the PEPS, see Fig. S5(a).The goal is to approximate the application of the gate to PEPS by another PEPS with the same bond dimension, see Fig. S5(b).Naively contracting the gate moves the PEPS away from the PEPS manifold, as shown in Figs.S5(c, d).In order to restore the PEPS, we need an update scheme which restores the form of the original two tensors from the PEPS used in the contraction.Two popular approaches to achieving this task are the Simple Update (SU) [41] and the Full Update (FU) [39,72,73].In general, to find the optimal truncation, which optimizes the overlap of the new PEPS with the untruncated TN, the environment of the pair of tensors (i.e. the contraction of the remaining PEPS tensors) needs to be taken into account.This is however a costly computation, and the SU includes only a primitive, but computationally inexpensive approximation of it as a product [38], and performs a truncated singular value decomposition (SVD) to split the resulting tensor, see Figs.S5(d, e).This results in a simple algorithm, with computational cost scaling as O(D 5 PEPS ).[75].In contrast, the FU takes into account the full environment, but requires the approximate contraction of the complete TN, as illustrated in Figs.S2 and S3.This update has a greater accuracy, but with a much steeper scaling of O(χ 3 In what follows, we demonstrate the effectiveness of both methods for the system sizes N = 10 × 10.For both approaches, we use an update schedule which reduces the time step from in the range δ ∈ [10 −3 , 10 −1 ].After many iterations of time evolution, we estimate the expectation value of the state, E = ⟨ψ|H s |ψ⟩.This process is repeated until we find convergence.Figures S6(a,  SSEP respectively, and various values of s.The circles show the results of the SU, and the squares the results of the FU.While it is clear in most instances the FU can provide significant improvements on the SU, it is worth noting that even the SU provides precise results for bond dimension D = 4, with errors δE ≲ O(10 −3 ).We find these errors to be sufficiently small, and thus proceed using only the SU to allow us to reach large system sizes.

Optimal sampling of dynamics
The long time statistics of dynamical observables are encoded in the deformed Markov generator W s .While in principle this object can be used to generate the trajectories which correspond to the statistics, it is difficult to do in practice due to the unnormalized nature of W s .That is, ⟨−| W s ̸ = 0, and in general, the leading eigenvalue θ(s) ̸ = 0.In the long time limit, we are able to overcome this difficulty using the so-called Doob dynamics, which maps the tilted generator onto a proper stochastic dynamics through the transformation W Doob and escape rates Rx = R x + θ(s).
While the Doob dynamics provides an efficient way to simulate the biased dynamics at long times, it is dependent on the fact that one has access to the leading eigenvector of the tilted generator.Using the PEPS optimization methods described above, we are able to estimate the left eigenvector ⟨l s | ≈ ⟨ψ s | P −1 , where ⟨ψ s | is our approximation to the leading eigenvector of H s .Note that retrieving ⟨l s | as a PEPS is simple due to the fact that P acts locally.This allows us to implement an efficient sampling algorithm which can estimate the dynamics Eq. (S15) by sampling l s (x) directly form our PEPS.Extracting l s (x) from the PEPS is done similarly to estimating the contraction of the PEPS networks in Fig. S2.The first thing to notice is that we can reduce the PEPS to a network of rank-4 tensor by specifying the value of the local index of each tensor, which is defined by the configuration x (that is, the configuration x specifies each local n k ), see Figs.S7(a, b).The values l s (x) are then retrieved by contracting the network.As was done for the networks in Fig. S3, we can estimate the exact contraction of the network using the boundary MPS method, with some bond dimension χ B .However, this time the boundary MPS are composed of rank-3 tensors, and a heuristic choice for the bond dimension is χ B ∼ O(D PEPS ), see Fig. S7(c).By contracting from two opposing edges of the network, we can then estimate l s (x) through the exactly contractable MPS-MPS product, see Fig. S3(d).
Unlike the networks in Fig. S2-which are composed of two PEPS layers-the sampling of l s (x) only requires us to contract over a single PEPS.This leads to a significant reduction in computational cost, with each calculation of l s (x) only costing O(N D 6 PEPS ).At each Monte Carlo step in the stochastic simulation algorithm, we need to calculate l s (x) for a maximum of N configurations, and thus the cost is O(N 2 D 6 PEPS ).However, by recycling our partial contractions when calculating each l s (x), it is possible to reduce this to a cost of O(N D 6 PEPS ) for each Monte Carlo step, see Refs.[56,73].While our approach provides a way to nearly optimally sample the Doob dynamics, the PEPS used in the sampling is only approximate.These errors must be accounted for through the use of umbrella sampling, see the main text and Ref. [56] for further details.
Figures S8, S9 demonstrate trajectories sampled for the 2D East and the 2D SSEP respectively, using the approximate Doob dynamics.The vertical plot in the left shows a light yellow line each time a transition occurs, and the right panels demonstrate a configuration snapshot at some points in time, as marked in the figures.We show (a) an active trajectory, (b) a trajectory close to the transition point s c and (c) an inactive trajectory.

Comparison to other methods
A full detailed comparison to the other methods via their numerical implementation is beyond the scope of this work.However, we will provide a brief discussion explaining how this approach compares to other approaches.One popular approach to estimating large deviations is cloning methods.This approach has an exponential cost in system size, and is known to suffer from bias.This is most apparent around the transition point, meaning that while the method allows for a way to probe dynamical phase transitions, it is not reliable for an accurate finite size scaling analysis as performed here.Another popular approach is trajectory sampling methods.While this approach is unbiased, it suffers an exponential cost in both time and space.Methods such as transition path sampling [69] can be used to hinder the cost, reducing the exponent in the exponential.Nevertheless, the cost is still exponential, and can be problematic where large times are required.
Tensor network approaches estimate the long time statistics of the dynamics by directly targetting the maximal eigenvalue of a deformed Markov generator.Each iteration of the optimization methods scales only polynomially in system size for a fixed bond dimension.In the case of time evolution, the number of iterations required to reach convergence is expected to scale as the inverse gap between the two leading eigenvalues of the deformed Markov generator.Furthermore, whereas the required bond dimension for a fixed precision is not known a priori, PEPS allow for a controlled way to systematically increase the accuracy of the method by increasing the bond dimension at a cost which is only polynomial in bond dimension.As shown in these works, the PEPS can be combined with trajectory sampling algorithms.This allows us to approximate the most optimal sampling dynamics at a cost which again scales only polynomially in both space and PEPS bond dimension.In practice, the dynamics is only approximate, and errors will still exponentially propagate in time.However, the prefactor is hugely reduced from the usual setting.

FIG. 1 .
FIG. 1. Models.(a) 2D East: an occupied site (black circles) facilitates flips in neighboring sites (red cells) only in two directions.(b) 2D SSEP: particles can hop to empty neighboring sites (black arrows), symmetrically in any direction; particles can enter or leave at the boundaries (red arrows).PEPS.(c) A PEPS is parametrized by a tensor per lattice site (red boxes for the top PEPS), each local tensor with a physical index (purple line) and four virtual legs (black lines) that connect it to neighboring tensors.The expectation value of a local operator (orange box) is obtained by sandwiching it between the PEPS and its adjoint (shown in green), and contracting (i.e.multiplying and summing over) the physical (basis) indices.(d) The cost of exact contraction scales exponentially with size, so an MPS approximation (blue tensors) is used for the contraction part of the network, with the dimension of its virtual bonds (blue legs) controlling accuracy; see Ref.[30] for details.

FIG. 2 .
FIG. 2. Dynamical large deviations and active-inactive transitions from PEPS.(a) The SCGF θ(s)/L 2 for the 2D East with c = 0.3 (top) and the SSEP (bottom) for system sizes N ∈ [10 2 , 40 2 ].The black dashed line shows the linear response for small s, and the colour dotted lines show the value for s → ∞.(b) The dynamical activity k(s)/L 2 for the systems in (a).The East is on a log-log scale, and the SSEP a log-linear scale.(c) The rate function φ(k)/L 2 as a function of activity k/L 2 for the systems in (a).The dashed line shows the Poisson distribution with mean k(s = 0)/L 2 .(d) The transition points sc(L) for the 2D SSEP (black circles) and the 2D East for c ∈ [0.2, 0.5].The solid lines show the fitted power-law curves sc(L) ∼ L −2α , with the exponents shown in the inset.The black dashed line is the exponent for the SSEP, and the symbols are for the East.The symbols can be used to read the value of c in the main figure.The bottom panel shows the dynamical susceptibility χ(s) = θ ′′ (s) for the 2D SSEP.All the data was acquired using the SU except for the black markers, which show (quasi-exact) 2D DMRG data for a N = 10 lattice for comparison.

FIG. 3 .
FIG. 3. Optimal sampling of trajectories.(a) Average dynamical activity as a function of s for the 2D East model from CTMC with importance sampling (symbols), for c = 0.5 and DPEPS ∈ [1, 4].The trajectory times are chosen such that on average we expect 100 transitions per trajectory.The solid black line shows the activity measured directly from the PEPS with D = 4 for comparison.Inset: variance in the time-integrated difference of escape rates, δR 2 (see main text).Each data point is calculated from Nsp ∈ [10 3 , 10 4 ] trajectories.(b) Same but for the 2D SSEP on a 22 × 22 lattice.(c) Representative trajectories for the 2D SSEP for size N = 10 × 10 at three values of s ̸ = 0.The bars on the left of each panel show the times when particle hops occur (yellow/bright lines).The snapshots on the right show the configurations at the marked times (black/white indicates a particle/hole).See also Refs.[30,50].
Figure 3(c) illustrates this by showing sample trajectories for the SSEP at s ̸ = 0: in the active phase (s = −0.1 panel), particle hops are plentiful, as shown in the bar on the left of the panel, and the configurations visited in the trajectory have a finite density and show no appreciable clustering; in contrast in the inactive phase (s = 1.0 panel), hops are scarce and the configurations have very low density of particles; near the transition (s = 0.1 panel), configurations show more pronounced density fluctuations, related to the critical nature of the transition.
S1) where R x = y̸ =x w x→y is the escape rate out of the configuration x.It is convenient to write the master equation in terms of a Markov generator, W = x,y̸ =x w x→y |y⟩ ⟨x| − x R x |x⟩ ⟨x| , (S2) which yields ∂ t |P (t)⟩ = W |P (t)⟩.The Markov generator W conserves probability, ⟨−| W = 0, with the flat state ⟨−| = x ⟨x|, and has maximal eigenvalue zero.
FIG. S1.PEPS.(a) The vector |ψ⟩ as a PEPS.The green cubes represent tensors within the PEPS, the (open) purple lines are their physical dimensions, and the (closed) grey lines are the virtual dimensions between tensors.(b) The same for its conjugate, ⟨ψ|.The conjugate tensors are represented by red cubes and black lines.(c) The two-body operator, Ô, as a tensor.(d) The probability amplitude ψ(n) can be retrieved by specifying the local dimension of each tensor, as shown by the grey cubes.This results in a network of rank-4 tensors, which give ψ(n) when contracted.
FIG. S2.PEPS networks.(a) The inner product ⟨ψ|ψ⟩ and (b) the expectation value Ô|ψ⟩ as TNs.For visual convenience, the operator and the tensors which it is connected to are show in full colour, while the remaining tensors are semitransparent.

FIG. S4 .
FIG. S4.Accuracy of the boundary approximation.We show the relative difference in the expectation value, ∆E(χB), when a boundary MPS bond dimension χB is used, compared to χB = 128.Results are for the 2D East model with N = 10×10, c = 0.5 and s = −1.0.
FIG. S5.Time evolution.(a) The PEPS |ψ⟩ can be updated at a local neighbouring pair of sites through the application of a Trotter gate (shown by the yellow cuboid).(b) This can be approximated by a PEPS with same bond dimension.(c)The updated tensors can be retrieved using the SU.The small turqouise cubes are the "λ-matrices" retrieved from SVDs of the neighbouring tensors (see e.g.Ref.[41]).(d) First, we contract the lattice site tensors and the surrounding λ-matrices with the Trotter gate.(e) Then, through an SVD, we can restore the PEPS manifold.The λ-matrices outside of the pair of tensors are restored to their original values, while the λ-matrix between the two tensors is updated.

s=
L [W s − θ(s)I] L −1 , where L = x l s (x) |x⟩ ⟨x| is the left eigenvector as a diagonal matrix, l s (x) = ⟨l s |x⟩.It is simple to check that the flat state ⟨−| is an eigenvector of W Doob s with the maximal eigenvalue zero.The Doob dynamics has the transition rates wx→y = l s (y) l s (x) e −s w x→y , FIG. S7.Optimal sampling from PEPS.(a) The component ls(x) = ⟨ls|x⟩ can be extracted from the PEPS.(b) The physical dimensions can be contracted to give a PEPS composed of rank-4 tensors.(c) As is the case in Fig. S3, the contraction of the PEPS can be estimated through a boundary MPS, this time composed of rank-3 tensors and bond dimension χB.(d) By contracting from two opposing edges of the network, we can reduce approximate the contraction as an MPS-MPS product.
FIG. S8.Representative trajectories for the 2D East.Trajectories sampled from the approximate Doob dynamics for L = 10, c = 0.5, and (a) s = −0.1,(b) s = 0.01 and (c) s = 0.1.The bars show the times when jumps occur (yellow/bright lines).The snapshots show the configurations at the marked times (black/white indicates a occupied/unoccupied).(a) (b) (c)