Quantifying Rare Events in Stochastic Reaction-Diffusion Dynamics Using Tensor Networks

The interplay between stochastic chemical reactions and diffusion can generate rich spatiotemporal patterns. While the timescale for individual reaction or diffusion events may be very fast, the timescales for organization can be much longer. That separation of timescales makes it particularly challenging to anticipate how the rapid microscopic dynamics gives rise to macroscopic rates in the non-equilibrium dynamics of many reacting and diffusing chemical species. Within the regime of stochastic fluctuations, the standard approach is to employ Monte Carlo sampling to simulate realizations of random trajectories. Here, we present an alternative numerically tractable approach to extract macroscopic rates from the full ensemble evolution of many-body reaction diffusion problems. The approach leverages the Doi-Peliti second-quantized representation of reaction-diffusion master equations along with compression and evolution algorithms from tensor networks. By focusing on a Schl\"{o}gl model with one-dimensional diffusion between $L$ otherwise well-mixed sites, we illustrate the potential of the tensor network approach to compute rates from many-body systems, here with approximately $3 \times 10^{15}$ microstates. Specifically, we compute the rate for switching between metastable macrostates, with the expense for computing those rates growing subexponentially in $L$. Because we directly work with ensemble evolutions, we crucially bypass many of the difficulties encountered by rare event sampling techniques$\unicode{x2013}$detailed balance and reaction coordinates are not needed.


I. INTRODUCTION
Connecting microscopic kinetics to emergent rates has a long history for equilibrium and nonequilibrium processes.Complex systems often settle into metastable basins with rare transitions regulating the rates of switching between the emergent states.In equilibrium, transition state theory [1,2], variational transition state theory [3], Grote-Hynes theory [4,5], and transition path sampling (TPS) [6,7] have given strategies to estimate the transition rate.Away from equilibrium, theories built upon a free energy landscape lose their applicability, yet trajectory sampling approaches can extract rates from microscopic dynamics, through noise-guided TPS approach [8,9] or forward flux sampling (FFS) [10][11][12].For all their merits, sampling approaches bring intrinsic challenges.How many samples are required to estimate the rate?Is a good approximate reaction coordinate needed, and if so, how is it found [13]?How does the expense grow with larger or more complicated microscopic systems?Motivated by those ubiquitous concerns, we set out to consider a different approach-one that effectively samples all possible trajectories as an ensemble evolution instead of propagating individual realizations.Such a proposition may appear absurd given the exploding state space presented by the many-body problem, but we show that for the class of discrete state, continuoustime, reaction-diffusion master equations (RDME), we can evolve the full distribution of a classical many-body systems in a numerically controllable manner by combining two well-developed methods, the Doi-Peliti (DP) formalism with tensor network (TN) algorithms.
While TPS and FFS have been developed for a va-riety of dynamical equations of motion (Hamiltonian, Langevin, Markov jump, etc.), we focus our attention solely on the chemical master equation (CME) dynamics, the setting where FFS was first developed [10].More specifically, we study discrete reaction-diffusion dynamics of L sites or voxels as they are often called in the RDME literature [14,15].The rationale for this focus is twofold.From a practical perspective, the RDME is important as a popular method to model chemical and biological processes [16][17][18].From a theoretical perspective, the RDME is interesting because it provides a setting where one can analyze the influence of stochasticity on the emergence of patterns [19][20][21].In many pattern-formation problems, it is possible to understand the underlying physics in terms of a mean-field partial differential equation that tracks spatiotemporal evolution of deterministically evolving fields [22,23].It is, however, well-appreciated that stochastic effects can quantitatively and even qualitatively impact kinetics [24][25][26][27][28][29].
It is therefore important that the methodology we discuss captures more than just the typical dynamics.Accurately computing transition rates between metastable states requires quantitatively resolving even rare tails of the state space distribution since these rare fluctuations can be instrumental in triggering transitions [11,30].
Historically, similar schemes for the CME ensemble evolution have been hampered by the curse of dimensionality.As a result, exact solutions are rare [31] and numerical solutions often resort to sampling methods built upon the kinetic Monte-Carlo stochastic simulation algorithm (SSA) [32,33], also known as the Gillespie algorithm.Another approach is to evolve a subset of the state space through the finite state projection (FSP) method [34].
FSP has previously been combined with TNs to approximately describe the CME [35,36], resembling our use of TNs to tame the many-body problem.Unlike those projective methods, however, the approach we describe uses the time-dependent variational principle (TDVP) to evolve the dynamics with probability conservation at all times, akin to prior work with diffusion and no reactions [37][38][39].Since we estimate rare events by measuring the distribution's small probability fluxes, conservation of probability and high-fidelity dynamical evolution are materially important.
While the formalism we describe applies quite generally to different well-mixed and reaction-diffusion models, we center our paper around a single model problem: rate calculations for transitions between basins of the Schlögl model [40] in the bistable regime, both with and without diffusion [41].In Sec.II, we first illustrate the rare event problem and extract rate constants from a propagator in the well-mixed regime.For more complicated systems (either systems with more chemical species or more voxels), the curse of dimensionality precludes the exact calculation of the propagator, but the remainder of the paper shows how to employ TNs to follow the same conceptual path.In Sec.III, we review the Doi-Peliti framework for writing reaction-diffusion dynamics in a secondquantized field theoretic form.Sec.IV then reviews the TN matrix product state (MPS) ansatz, a controllable approximation to the ensemble that can be evolved via TDVP.Those methods are applied to the many-body rate calculation in Sec.V.In Sec.VI, we show results for the RDME with diffusion between L well-mixed Schlögl voxels, demonstrating that the methodology can extend numerical transition rate calculations to regimes where the SSA becomes impractical.The approach effectively evolves the probability of occupying each microstate in an unfathomably large many-body state space, yet the calculations scale sub-exponentially with system size and do not assume detail-balanced equilibrium.Sec.VII discusses how the TDVP calculations yield not only a rate, but also a mechanism.Finally, in Sec.VIII, we discuss potential applications and challenges that should motivate future work.

II. ENSEMBLE RATE CALCULATIONS
Consider an archetypal two-state rate process between reactants and products R −− −− P, initiated entirely in the reactant state.After some reaction timescale, the concentration of products grows before eventually leveling off at its steady-state value.The reaction rate can simply be defined and measured in terms of how quickly the product concentration rises.An alternative formulation is to study the fate of an individual reactant on its path to transitioning into a product.The TPS and FFS path sampling approaches extract the reaction rate by focusing on the statistical properties of these individual trajectories.Algorithms exist to propagate those single trajecto-ries, algorithms that can be practically implemented even when the system grows large and complex.By contrast, propagating the entire ensemble of trajectories has typically been limited to simple low-dimensional problems.In this section, we illustrate an ensemble-evolution rate calculation for a simple case in order to develop notation and lay the groundwork for the many-body problem.The remainder of the paper will show that the methodology can be extended beyond such simple problems by employing TNs.

A. Well-mixed Schlögl Model
Consider the well-mixed Schlögl model, a quintessential example of bistability in chemical reaction networks (CRNs) [42,43].The model tends to be studied under chemostatted conditions that fix two species' concentrations A and B, while allowing species X's concentration to evolve according to the four reactions: where c1 , c2 , c3 and c4 are the stochastic rate coefficients.Our starting point is microscopic such that the rate c1 should be interpreted as the reciprocal of the waiting time for a particular set of 2 X molecules to react and form 3 X molecules when mediated by a fixed concentration of A molecules.In the large-system, well-mixed limit of chemical kinetics, it is common to also define macroscopic, or kinetic rate constants k r = cr (N A V ) m for each reaction r with m reactants.Here, Avogadro's number N A serves to convert between microscopic measures (counting molecules) and macroscopic kinetic measures (concentrations in molarity) for a volume V .In this paper, we focus on the microscopic kinetic picture, viewing the cr as the fundamental parameters, but these parameters can be converted into k r for linear reactions or in the thermodynamic limit [44].In that thermodynamic limit, the stochastic kinetics often tends to the solutions of the deterministic reaction rate equations.For the Schlögl model, analysis of that deterministic limit reveals that the sign of the discriminant determines the number of stable steady-state solutions [42].While holding the rate coefficients fixed, the sign of ∆ changes as a function of the concentrations of A and B, allowing one to pass from ∆ > 0 with one steadystate solution, through ∆ = 0 with a diverging variance of the number of X molecules [45], and on to ∆ < 0 with one unstable and two stable solutions.In the ∆ < 0 regime, a finite (stochastic) system sees the regions around the two stable fixed points become metastable states.The After the so-called molecular timescale τ mol , the fastest trajectories begin to reach B. Subsequently, probability leaches into region B at a constant rate with P (Bt|A0) ∝ kBA.The inset shows that growth on a much longer timescale, revealing the third scale when the system eventual approaches the steady state.number of X molecules stochastically switches between the two metastable values.A representative trajectory, shown in Fig. 1, oscillates between a state with n ≈ 5 X molecules and another with n ≈ 50.We first set out to compute the rate of bistable switching between those two metastable states in the well-mixed regime, making use of the fact that the dynamics is one-dimensional (1D)-n executes a random walk between 0 and some large maximum number M , which we impose for the sake of calculation.
Rather than track n(t) for a single trajectory as in Fig. 1, the ensemble approach tracks evolution of the distribution p t (n), measuring the probability that an individual trajectory would have n molecules at time t.Foreshadowing the analogies with quantum-mechanical methodologies, we choose to write that distribution as a ket: |p t = [p t (0), p t (1) . . ., p t (M )] T , where T indicates the transpose to give a column vector.The distribution evolves according to the master equation Here, H is a rate matrix whose off-diagonal elements H n n are the rates of transitioning from n X molecules to n and whose diagonal elements H nn are minus the rate of escaping from the state with n X molecules.The rate matrix conserves probability, a fact compactly expressed as 1| H = 0 1|, with 1| denoting a vector of ones.The distribution |p t can always be formally evolved in terms of a propagator e Ht as |p t = e Ht |p 0 .
Due to the analogy with quantum-mechanical time evolution, we equivalently call H an effective Hamiltonian.Assuming H is irreducible, the long-time limit for this distribution tends to the unique steady-state |π = lim t→∞ e Ht |p 0 .
H can be decomposed into contributions from the two different reversible mechanisms of Eqs. ( 1) and ( 2), labeled H A and H B , respectively, based on whether the transitions are mediated by the chemostatted species A or B.Then, the Schlögl model at the level of the master equations corresponds to tridiagonal rate matrices with elements and where H = H A + H B and δ i,j is the Kronecker delta.
Given the specific tridiagonal H and the restriction 0 ≤ n ≤ M , both H and the propagator e Ht are (M + 1) × (M + 1) matrices, amenable to evolving the entire distribution, effectively averaging over all the realizations.
To extract a rate from the evolution of distribution, we must initialize the distribution in one of the two metastable states and measure the rate of relaxation.In analogy with free energy basins, we will call the two metastable regions basins A and B (red and blue regions of Fig. 1 respectively) with A ∩ B = ∅.We write Pω for the non-negative idempotent ( P2 ω = Pω ) operator that projects onto the density lying in a region ω, where ω ∈ {A, B}.If |p t has support over a state space Ω, then the projection Pω |p t is an unnormalized distribution with support restricted to ω.The total probability in ω at time t can then be expressed as 1| Pω |p t .

B. Rate Constants
In general, one can compute the probability to pass from any subset A ⊂ Ω into an arbitrary non-overlapping subset B. For the current example, A and B are basins of attraction for the Schlögl model (see Fig. 1a), and we would like to measure the rate of going from A → B. Rate constants can be defined either for the unidirectional flux from A → B, the quantity we focus on in this work, or for the timescale of the relaxation to the steadystate, which also involves the B → A transitions [46].Additionally, the precise definition of the rate depends on whether one seeks trajectories that move from the boundary of one basin to another or from the interior of one basin to another.For instance, studying the time to transition from the boundary of A to the boundary of B yields the Transition Path Theory (TPT) rate associated with the committor function [47].This unidirectional rate k BA is related to the probability of having reached B from A in time t without returning to A. If we focus on trajectories short enough that those reaching B will stay (no recrossings) then k BA can be expressed in terms of the probability P (A 0 ) to start in A at time zero and the conditional probability P (B t |A 0 ) to subsequently be found in B at time t: The derivative is evaluated at some timescale exceeding τ mol , a molecular time scale of the system that is long enough for the first trajectories to reach B from A but much shorter than the typical transit time k −1 BA (see Fig. 1b).
We compute the rate of Eq. ( 8) by initializing in a distribution confined to A then propagating for time t and measuring the probability in B. A suitable initial distribution is PA |π / 1| PA |π , a normalized distribution constructed by projecting the steady-state distribution |π onto A. With this contribution, P (A 0 ) = 1, and Eq.(8) becomes To the extent that there is a separation of timescales between exploring a basin and transiting between basins, the transition event is rare and no transitions can occur faster than the molecular time scale τ mol .Beyond that time, the number of transitions accumulates over time as P (B t |A 0 ) ∝ k BA t.At times long compared to k −1 BA , P (B t |A 0 ) must eventually stop growing linearly to level off at the stationary value 1| PB |π , but our interest is in extracting k BA from the growth at much shorter times.A typical realization of the process will not require observation much beyond τ mol units of time, to ascertain the rate, meaning our ensemble evolutions must be propagated for roughly τ mol .Notice from the inset of Fig. 1(b) that at long times the recrossing events obscure the rate of growth, ultimately yielding the plateau.As discussed in Sec.V and App.D, for rare rates and long times, the rate calculations can be improved by removing these recrossing events via a modified dynamics in which region B is an absorbing state, forbidding dynamics out of that region.
Our ability to compute k BA from the ensemble evolution traces back to the fact that it is practical to numerically calculate e Ht for the well-mixed model with a single species and a modest M .While Eq. 9 generalizes naturally for harder problems (those with more than one species or with diffusion between multiple well-mixed voxels), the state space for |p t grows exponentially with additional stochastic degrees of freedom.We next show that the Doi-Peliti framework re-expresses the rate matrix in a second-quantized form that allows us to generalize beyond the case of a single well-mixed species.

III. DOI-PELITI
The Doi-Peliti (DP) formalism is a classical version of the second quantization methods from quantum field theory [48,49].Originally introduced by Doi [50] as a critique to the assumptions of the Smoluchowski equation, it was later developed by Peliti [51] and others [52] as a path integral formalism for solving chemical reaction networks [53][54][55][56].The DP framework has also found success in other areas such as studying kinetically constrained lattice models in the context of glassy dynamics [57,58], cellular signaling [59], predator prey models [60] and active matter [61].The DP representation builds operators which encode the dynamics onto the high-dimensional space that accompanies many-body problems.Crucially, writing the chemical reactions in an operator form guarantees it can be easily converted into a matrix product operator (MPO) [62,63], as we will do in Sec.IV.
To start, we go beyond Sec.II by considering wellmixed reactions with multiple dynamic species as well as heterogeneous systems where molecules can diffuse between neighboring voxels.Let X 1 , X 2 , . . .X L represent each chemical species in each site.We could be talking about a single species X with different values in L different voxels, a single well-mixed voxel with L distinct dynamic species, or a mixture with multiple species and multiple voxels.The number of molecules in a given voxel can change due to a chemical reaction or a diffusion event, with both events being cast as, where cr are again the rate coefficients for reaction r.η r j and µ r j reflect the stochiometry of the reactions and give the stochiometric vector ν r = µ r − η r which evolves the system from state n ∈ Ω to n = n + ν r .In the previous section, these reactions induced changes between the M + 1 different microstates, specified by the scalar n.Now, the microstates n are specified by n = [n 1 , n 2 , . . ., n L ], counting the occupation of each species within each voxel.
For a dilute well-mixed solution, where well-mixed means molecules locations within a voxel are unresolved, the change in probability of microstate n follows the gainloss equation [14], where the sum is taken over N different reactions of the form Eq. ( 10), each with their own stochiometric vector ν r .The transition rate between any two states is the sum of contributions from each reaction H n n = r α r n .Each propensity α r is a product of a stochastic rate coefficient cr and the combinatorial number of ways to select either η j or µ r molecules from reaction r, depending on the direction of the reaction.For a "forward" reaction the propensity is, α r (n) = cr nj ∈r nj ηj .
Encoding the reactions into a matrix form allows us to evolve the elements of the joint distribution p t (n) analogously to the previous section.While it is appealing to try and directly solve these more complicated dynamics as we did before, it is impractical.The curse of dimensionality demands a less onerous way to work with both H and p t (n).Since particles of the same species in the same voxel are indistinguishable, it is beneficial to work in an occupation number representation with The basis vectors corresponds to a possible microstate, and the set of basis vectors forms a tensor product space called a Fock space F L with |n ∈ F L and m|n = δ mn [64].Without loss of generality, we will discuss the case of a single chemical species and L physical voxels, in which case n l has the interpretation of the number of molecules in voxel l.The second-quantization procedure rewrites the rate matrix as an operator Ĥ in terms of creation and annihilation operators, x † l and x l , that respectively raise and lower the number of X molecules in voxel l.That is to say, where |0 is the vacuum state, i.e., no molecules.By x l or x † l we are implying that the operator only acts on the l th voxel: Because they act on the single-site space, the general onesite operators x and x † can be written in a matrix form as (16) That these matrices yield Eq. ( 14) can be readily confirmed.More abstractly, the existence of such creation and annihilation operators in this classical setting originates from the classical commutation relation [x i , x † j ] = δ ij and the existence of a positive semi-definite number operator nl = x † l x l , so named because nl |n l = n l |n l .Together, the set {Ω, F L , x, x † } are also sometimes called an interacting Fock space [48].As in the construction of the number operator, we will be able to decompose contributions to Ĥ in terms of the x † and x's.
In practice, the matrices of Eq. ( 16) will not be infinite in size but will rather be (M +1)×(M +1) matrices due to our enforcement of a maximum voxel occupancy of M .At the operator level, this maximum voxel occupancy supplements Eqs. ( 14) with an additional requirement that a requirement that is met by the finite truncation of Eqs.(16).Those truncated operators obey a commutation relation which is similar to that of their infinite counterpart, but the usual commutation relation is broken at the maximal occupancy state.That is to say The commutation relation for the truncated operators is therefore an expression involving the projector operator onto site i's maximal occupancy state, PMi = |M i M i |.As discussed in depth in App.A, this violation of the usual commutator relation must be accounted for when constructing an effective Hamiltonian that conserves probability.
One benefit of defining operators in terms of creation and annihilation operators is that expectation values for observables can be expressed easily in terms of inner products with operators.Consider, for example, the observable that counts all particles in all voxels, n = l nl .The operator n acts on a many-body state |n , yet each nl is a single-body operator acting on voxel l and leaving all other voxels unchanged.It is therefore understood that when local operators such as nl act on |n , it is written as a shorthand for In terms of that shorthand, nl |n = n l |n counts the number of particles in voxel l for the many-body state |n .The time-dependent expected number of particles in the whole system is thus n(t) = l 1| nl |p t , where the vector of ones 1| serves to sum over all microstates.The time-dependence of a distribution can also be extracted from the action of a second-quantized Ĥ on |p t , and this operation can be practically computed when Ĥ is expressed in terms of the single-site operators.
Mirroring the discussion in Sec.II, the generator Ĥl for the Schlögl model reactions within the l th voxel is a sum of contributions from the elementary reactions Ĥl = r Ĥr .Fig. 2a shows each reaction, expressed as a second-quantized contribution to the effective Hamiltonian, as derived in App. A. Probability conservation requires that this effective Hamiltonian include a term with a projector Pc M l ≡ I − PM l that prevents flows of probability into forbidden microstates with occupancy greater than M .This projector is complementary to the one appearing in Eq. (19) in that it projects onto states which are not maximally occupied.Combining together the forward and reverse reactions mediated by species A, we show in App.A that the reversible reaction Eq.(1) contributes to the effective Hamiltonian.For convenience we have absorbed the combinatorial terms into that rate coefficient, defining c r ≡ cr /γ!, where γ is either η r or µ r for forward and reversed reactions, respectively.Likewise, the contribution comes from Eq. ( 2).In the M → ∞ limit, Pc M l → I to give the more familiar expressions written in Fig. 2a.
The Hamiltonian operator for L independent, wellmixed Schlögl voxels is thus Ĥrxn = L l

ĤA
l + ĤB l .To add diffusion, consider the well-mixed voxels arranged in a 1D lattice as in Fig. 2b.Each chemical species is able to move between neighboring voxels as regulated by a diffusion operator with hopping rate d = D/h 2 , consistent with macroscopic diffusion constant D for voxel width h [18].Combining reaction with diffusion gives a straightforward linear superposition at the level of the effective Hamiltonian.We emphasize that the linear superposition of diffusion Ĥdiff with reactions ĤA l and ĤB l does not imply a restriction to linear kinetics.The Schlögl model kinetics is nonlinear in that it involves elementary reactions which are not first order, and the diffusive coupling between neighbors can induce nontrivial interactions between voxels.

IV. TENSOR NETWORKS
We have now seen how to transform the matrix-form rate operator H of Sec.II into a Doi-Peliti form Ĥ, expressed in terms of local creation and annihilation operators.That change of representation is not merely cosmetic.As the number of chemical species or number of lattice sites grows, the size of H grows exponentially.For N dynamical (unchemostatted) species and L lattice sites, the state space has (M + 1) N L states.To put this scaling in perspective, the largest numerical result that FIG. 3. The ensemble evolution of the many-body state requires that one solve for the probability of every microstate at all times.That high-dimensional joint probability distribution can be well approximated by a matrix product state (MPS).
An MPS approximation for the joint distribution can be advantageous whether there are Ns species in a well-mixed reaction follows uses N = 1, M = 85 and L = 8, meaning the system contains 2.99 × 10 15 microstates.Calculations on a matrix H, even a very sparse matrix H, are untenable.By contrast, it is practical to construct a TN decomposition of the second quantized form of Ĥ in terms of an MPO, as shown in Fig. 3b, and demonstrated explicitly for the Schlögl model in App.B. Given the second quantized form, the MPO can be readily constructed in a practical form through SVD decomposition using automated packages such as AutoMPO in ITensor [65].
Much like the operator Ĥ can be compressed as an MPO, the states that operators act upon can also be expressed more compactly in terms of TNs.Consider p t (n) = p t (n 1 , n 2 , . . ., n L ), the probability of microstate n at time t.This p t (n) can be thought of as a rank L tensor.By specifying the occupancy of each of the L sites, p t (n) outputs a number-a probability-associated to that microstate.By leveraging a matrix product ansatz, p t (n) can be approximated as a partial contraction over L different low-rank tensors (see Fig. 3a).The low rank tensors of an MPS [66] are organized in a 1D chain of site dependent tensors whose upper index n l specifies the occupancy of site l and whose lower indices are so-called bond indices, nonphysical indices which are contracted over to generate correlations between nearby sites.For compactness of notation, we suppress labels Q (1) , Q (2) , etc. that would distinguish the site dependent tensors, implicitly carrying that site dependence through the labeling scheme for the physical indices.For suitably chosen Q tensors, we therefore approximate where |q t is the approximate distribution over microstates generated from the MPS.Though not explicitly written, the Q tensors are time-dependent.
There is a diagrammatic representation of the tensors, shown in Fig. 3, which aids in visualizing tensor operations.The tensors Q are given by circles, one per site, and vertical lines feeding into these circles are the physical indices that pick out a particular slice of the tensor based on the occupation number of that site.Circles also have horizontal lines emanating from them, representing the non-physical bond indices.Connection of two circles by a line represents a shared index being contracted over.Fig. 3 emphasizes that the TN decomposition is flexible enough to address varied CME's.The MPS of Fig. 3a middle, has been the focus.Each site from 1 to L reflects the number of X molecules in that site, but Fig. 3a (top) and (bottom) show that the MPS decomposition can equally well be constructed when multiple species X i , i = 1, 2, . . ., N s can occupy a single site or when N s different species diffuse between L/N s sites to yield an MPS with L sites.What is practically important is that Q's be arranged such that the physical dynamics correlates the neighbors.For example, for the RDME Schlögl model, neighboring Q's correspond to neighboring voxels in the 1D Schlögl lattice.
An attractive feature of the MPS ansatz is that the size and accuracy of the MPS are controllable through the bond dimension of the auxiliary indices s.If the dimension of those indices is allowed to grow exponentially from s 1 through s L/2 , the ansatz in Eq. ( 25) is exact, but even when the bond dimension of the indices s are capped at a maximum dimension of κ, the ansatz can be a very good approximation.Crucially, capping that maximum bond dimension makes it practical to work with |q t because the rank-L tensor |p t has been replaced by the L different Q tensors, each with no more than (M + 1)κ 2 elements.
Having approximated both |q t and Ĥ by an MPS and MPO, respectively, we can now revisit the otherwise intractable dynamics problem of Eq. ( 4), only now we wish to evolve |q t in lieu of |p t .As |q t is an MPS state, we wish for the dynamics to be constrained for all times to the manifold M of possible MPS states.Merely replacing |p t by |q t in Eq. ( 4): would not impose such a constraint because Ĥ in general evolves |q t off M. An appealing resolution is given by the Time-Dependent Variational Principle (TDVP) [67,68], which proposes to first project the right-hand side onto the tangent space of the MPS manifold: In this way, dynamics is restricted to the desired subspace of MPS states [69].It turns out that the time-discretized numerical integration of this |q t is amenable to efficient TN operations.The technical details of the time-evolution algorithm have been reported elsewhere [68,[70][71][72].We give only a brief, high-level overview of that procedure (see also App.C for more details).The central step in deriving the TDVP time-evolution algorithm is to demonstrate that the projector PT M can be decomposed in terms of a sum over L terms (associated to each site of the MPS) and L − 1 terms (associated to each bond connecting sites).In the limit of a small time step δt, the propagator e PT M Ĥδt can thus be evaluated as a composition of time-evolutions involving the 2L − 1 individual terms.That composition of time-evolutions is executed by a one-site TDVP algorithm that "sweeps" from lattice site 1 through L, advancing each term for the small time step.
The efficiency of the TDVP time-evolution algorithm relies on the flexibility to represent the same MPS state in redundant ways due to a gauge freedom.Consider, for example, neighboring tensors Q n l s l−1 s l and Q n l+1 s l s l+1 in an MPS.For an invertible R, replacing those tensors by Q n l R and R −1 Q n l+1 , respectively, does not change the MPS state; upon contracting over the neighboring sites, the identity I = RR −1 would fall out.By leveraging the gauge freedoms, it is always possible to transform Eq. ( 25) into a so-called mixed canonical form [73]: where each A is a left-orthogonal tensor satisfying A T A = I and each B is a right-orthogonal tensor satisfying BB T = I.For notational compactness, we have suppressed implied summations over the bond indices between tensors in Eq. (28).We have retained the physical indices n 1 , . . ., n L and allow their presence to remind us that although tensors at each site are named merely A and B, they are actually distinct tensors A (1) , A (2) , . . ., B (1) , B (2) , . . .that vary across sites.The mixed canonical form is said to be centered about the tensor Q n l .From that tensor's perspective, the chain of A's and chain of B's constitute the environment, and it is beneficial to contract over the links between A's and links between B's to obtain so-called environment tensors.These environment tensors inherit an orthogonality from the orthogonality of the A's and B's, significantly reducing the computational complexity to advance Q n l in time by δt provided the MPS was already centered about the Q n l we aim to evolve.
As made more explicit in App.C, the one-site TDVP algorithm initially centers the MPS around the first tensor, propagates that tensor forward in time by δt, performs a QR decomposition to change the gauge, propagates the R term of that QR decomposition backward in time by δt, then passes that term to the neighboring site to re-center the MPS around the neighboring site.The neighboring site is then propagated forward in time as the algorithm iterates.After the complete sweep through the MPS, the probability distribution |q t has been propagated by a time step δt.Whether that propagation well approximates the evolution of the full distribution |p t depends on the time-step error but also on the severity of the MPS approximation, which can be systematically controlled through the maximum bond dimension κ.
Notably, the time-evolved MPS state |q t allows us to estimate the probability of any microstate of interest, even though it would be impossible to enumerate probabilities of all microstates.Alternatively, sets of microstates can be summed over to yield tractable marginal distributions.As an example, Fig. 4 shows the results of calculating q t (n 1 , n 2 , n 3 ) for a 1D chain of L = 3 Schlögl sites.Though there are three sites, site one is traced over to leave the two-dimensional joint probability of finding n 2 X molecules in site 2 and n 3 X molecules in site 3.At time zero, the distribution is initialized within region A according to PA |π / 1| PA |π .As time progresses, a small flux of probability leaks from A to B. Because the overwhelming preference is to stay in A, Fig. 4  state approximation, while the population in B gradually amasses, the distribution over microstates between A and B is nearly stationary.The slight asymmetry between sites is due to the fact that site 2 is in the bulk while site 3 lies on the boundary.The higher-dimensional joint q t (n 1 , n 2 , n 3 ) has the anticipated symmetry under interchange of n 1 and n 3 .More generally, symmetry holds between interchanging sites l and L − l + 1, as reflected in Fig. 6.In the next section, we illustrate the fact that the MPS state's ability to capture the rare time-dependent flux is key to calculating the rate of traversing between basins.

V. MANY-BODY RATE CALCULATION
Having established the TDVP dynamics of an MPS state, we can now repeat the strategy of Eq. ( 9) to compute a transition rate between states A and B: In doing so, there are a few important differences from Eq. ( 29).Most obviously, the TDVP dynamics includes the projection onto the MPS tangent space and |π must now be an MPS approximation to the steady-state distribution.Additionally, the meaning of PA and PB must be adapted to now define many-body analogs of the projectors onto the two metastable states.The single-voxel Schlögl model had an A state defined as all microstates with occupancy n above a threshold λ * A .Similarly, B included all microstates with occupancy less than some other threshold λ * B .The L-voxel Schlögl model generalizes to involve hypercubic metastable states with basins defined by A ≡ {n l ≥ λ * A ∀l} and B ≡ {n l ≤ λ * B ∀l}.
With these generalizations, the transition rate of Eq. ( 29) would measure the flux crossing into B through its boundary, but if we are to compare with a transition path theory rate, we must exclude re-crossing events.In other words, we want to include trajectories which exit A then some time later enter B, but we want to exclude the flux that leaves B and re-enters B some time later without having returned to A. When k BA is sufficiently large, those re-crossings can be ignored, but when k BA becomes small enough, those recrossings can become significant and must be removed.We note that the TDVP dynamics can just as well be executed with a modified absorbingboundary-condition dynamics with effective Hamiltonian H = Hrxn + Hdiff , where the reaction and diffusion Hamiltonians exclude events that exit B. An explicit construction of that modified absorbing-boundary dynamics is provided in App.D, allowing us to directly remove the re-crossing events.

VI. NUMERICAL EVALUATION OF THE L-VOXEL SCHL ÖGL MODEL SWITCHING RATE
The formalism of the preceding sections has laid out a controllable approximation to extract rates between metastable states of RDMEs.Numerical experiments are necessary to evaluate whether that formalism can be practically useful.We set out to demonstrate that utility by computing the transition rates between highoccupancy and low-occupancy states of a L-voxel Schlögl model as a function of L. The model is parameterized by the four stochastic rate coefficients discussed in Sec.II.We sought rate coefficients that fall in a bistable regime, with sharp peaks of steady-state probability falling in high-and low-occupancy states.Since the dynamics The rate kBA for switching between metastable states of the L-voxel Schlögl model is calculated using the SSA (red line) and using TDVP (dashed lines) with varying bond-dimensions κ using the same rate parameters as Fig. 4. For sufficiently large bond-dimensions, TDVP agrees with a SSA rate estimate built using 400 transitions from A to B at each L. (b) As L grows a larger bond-dimension is required to capture the dynamics.Signatures of numerical instability are apparent when κ is too small, shown here for L = 3, which stably converges when κ exceeds 6.For this system, κ = 20 is sufficient for numerical stability up to L = 6, while κ = 30 gives stable solutions up to L = 8.(c) The total number of CPU hours to estimate kBA grows exponentially with SSA but polynomially with TDVP.
would need to be truncated at some maximum occupancy M , we furthermore wanted to limit how many molecules would naturally accrue in the high-occupancy state.By scanning parameters with the L = 1 Schlögl model, we found parameters (c 1 = 2.676, c2 = 0.040, c3 = 108.102,and c4 = 37.881) compatable with an M = 85 truncation and with a sharp bistability between A and B regions defined by the thresholds λ * A = 15 and λ * B = 25.By linking multiple voxels together in a 1D chain with a diffusive hopping rate d, it becomes rarer to switch between the two metastable states than it would be with a single isolated voxel.With multiple voxels, a switching event must flip the state of each voxel, and diffusion between neighboring voxels works to stabilize neighbors in identical high-or low-occupancy states.Due to that diffusive coupling between neighbors, the L-voxel Schlögl model's A to B transition is characterized by a correlated flip of all voxels, and the rarity of that correlated event depends on L, the number of correlated voxels.It is therefore straightforward for us to push the ratecalculation methodologies into challenging regimes; k BA can be decreased by growing L.
The most straightforward way to numerically detect the L-dependence of the transition rate is via sampling of trajectories.Realizations of the continuous-time RDME, subject to the M = 85 maximum voxel occupancy at fixed L, were generated using the SSA.First, 5×10 5 independent SSA realizations were initialized from a uniform distribution and propagated in parallel up to t = 5 to approach the steady state distribution.Sampled configurations that fell within A were then sampled as initial configurations for 400 independent trajectories which were evolved via the SSA until the trajectory first reached B. From the 400 trajectories, the transition rate was computed from the mean of those waiting times to reach B: k BA = τ BA −1 .Fig. 5a (red line) shows that the rate exponentially decreases as L is increased, a result consistent with the exponential dependence of rate on volume in the well-mixed Schlögl model [74].While the rate can be computed by trajectory sampling, the difficulty to do so (via brute force sampling) grows in proportion to the typical τ BA , meaning the computational cost grows exponentially with L (see Fig. 5c).
The poor scaling of brute-force SSA sampled rate calculations has motivated so-called advanced sampling methods like FFS [10], which focus computational effort on the progress along a reaction coordinate.The idea is to introduce a reaction coordinate λ that resolves progress in transitioning from a value λ * A at the boundary of A to λ * B at the boundary of B. Rather than seek long trajectories that snake from λ * A all the way to λ * B , it can be beneficial to break the transition up into smaller trajectories that make some progress in moving to states with values of λ closer to λ * B than they started.This FFS scheme is particularly clever in that it can be statistically unbiased, meaning any choice of reaction coordinate (any mapping taking a microstate to a value of λ) will yield the true rate in the infinite sampling limit.Despite that formal lack of bias, the computational cost of a FFS calculation depends strongly on the correlation between λ and the reaction's progress [11].If one already has a good idea of how the reactions rare transition events proceed, a smart λ can be constructed and the exponential scaling of the SSA can be avoided, but that identification of the reaction coordinate is seldom straightforward.For example, in the L-voxel Schlögl model, will A to B events typically be triggered by a high-to low-occupancy flip at a boundary voxel of the 1D chain or in the bulk?The boundary voxels have the benefit that they only have diffusion from one neighbor stabilizing the high-occupancy state, but the bulk voxels have the advantage of being more plentiful.In App.G we describe FFS calculations employing two reasonable choices of reaction coordinate.Those existing FFS methods can yield rates more cheaply than a brute force SSA calculation, but the speed up depends strongly on the reaction coordinate, the choice of which is nontrivial even for this simple model.
The DP/TDVP approach we have described offers a radically different strategy for circumventing the SSA scaling problem.Because the reaction-coordinate selection problem is hard, we developed the DP/TDVP strategy to bypass it entirely, computing the rate and reaction mechanism without prior information about how the transition events move through the high-dimensional state space.Fig. 5 shows that the reaction-coordinateagnostic DP/TDVP approach indeed reproduces the SSA rates at lower computational expense, reducing the scaling from exponential to polynomial in L. The SSA rate calculations can be made more accurate by collecting more trajectories, and the DP/TDVP calculations likewise offer a systematic way to improve accuracy in exchange for a greater computational expense.The two most prominent ways to adjust that accuracy are through the time step δt for the TDVP evolution and through the maximum bond dimension κ for the MPS state.In this work, a time step of δt = 10 −4 was used throughout, and the trade-off between expense and accuracy was entirely tuned through the choice of κ.The single-site TDVP algorithm we used to propagate an MPS state does not alter the bond dimension of the MPS, so κ was entirely set by the bond dimension of the initial MPS.We therefore needed to approximate the steady state |π with different values of κ.Approximations were generated by starting with an initially uniform state using a bond dimension κ set to several different values between 5 and 30.Those uniform distributions were evolved with single-site TDVP for 10,000 time-steps allowing for the natural dynamics of Ĥ to approach a stationary state with p t | ṗt = p t | Ĥ |p t ≈ 0 by time t = 1.Irreducibility of Ĥ means the unique stationary state is the steady-state distribution |π .While p t | Ĥ |p t = 0 would imply |p t = |π , p t | Ĥ |p t ≈ 0 does not guarantee that |p t ≈ |π , only that the distribution is evolving very slowly.When transitions between A and B are slow, one can imagine reaching a local steady state distribution within A and within B but having probability very slowly repartitioning between the two.For our purposes, what matters is not that t = 1 has exactly reached |π but rather that the t = 1 distribution confined to A, upon renormalization, has approached PA |π / 1| PA |π .Provided dynamics within A is comparatively fast, this local steady state within A is reached more rapidly than |π and is reasonably confirmed by checking that p t | Ĥ |p t ≈ 0.
The renormalized t = 1 distribution was then propagated under single-site TDVP with the H absorbingboundary dynamics (see App. D) for another 10,000 time steps, and projected onto B to add up the total probability to have reached B. Note that these 10,000 steps are not nearly enough to reach the steady state |π because the slow rate of transitioning from A to B is a bottleneck.As additional time steps of dynamics were taken, the growth of the transition flux from A to B was measured and plotted in Fig. 5b.That figure shows that at very short times (t < τ mol ), the flux from A to B has not reached a steady state value because there hasn't been enough time for any transitions.After that initial time τ mol , the plot plateaus at a value corresponding to the rate k BA , a plateau that is reached orders of magnitude faster than the time for a typical SSA trajectory to make a transition.We find that if the bond dimension is insufficient (e.g., κ = 5 for L = 3 in Fig. 5b), then the dynamics becomes unstable, reflected in erratic estimates for the conditional probability P (B t |A 0 ).As the bond dimension is grown (e.g., κ = 9 for L = 3 in Fig. 5b), the estimates converge, allowing a stable rate to be extracted.That convergence as a function of κ is analyzed further in App.E. We find that setting the bond dimension to roughly 10 or 20 is often sufficient, but larger values of L can only converge with a larger bond dimension.For example, Fig. 5a shows that SSA rates agree with κ = 20 calculations up to L = 6, but κ = 30 was required to push up to L = 8.Provided that larger bond dimension is used, the DP/TDVP rates reproduce the SSA rates across five orders of magnitude, up to and beyond the point that the brute-force SSA rates are practical.
Though the required bond dimension grows with L, Fig. 5c shows that the expense to estimate the rate grows sub-exponentially.The blue and black dashed lines in that plot are the total number of CPU hours time to estimate the rate using TDVP with κ = 30 and κ = 20 respectively.Both SSA and TDVP were run using codes written in the Julia programming language [75] on equivalent CPUs with details of the runtime analysis described in App.F. While the required CPU hours, of course, depends on details of the hardware, parallelization, the particular implementation, we illustrate the difference in scaling when comparing identical implementations run in identical ways across different system sizes.Fig. 5c clearly demonstrates that difference in scaling between brute-force SSA and DP/TDVP.An increase in bond dimension increases the computation time, but it does not fundamentally alter the favorable scaling in L. Comparison of computational expense with FFS is complicated by the fact that one can optimize FFS in many different ways (how many interfaces, which reaction coordinate, etc.), but App.G shows that the DP/TDVP calculations remain favorable even compared to FFS with a reasonable guess of an order parameter.

VII. DEPLOYING THE TOOL
We have described a numerical tool to extract the rate of a rare switching event from a reaction-diffusion model, but the numerical value of that rate is often not the final goal.More than the value of the rate, it is desirable to know how the rate changes in response to controllable parameters.Fig. 5a showed one such variation-the switching rate falls off exponentially with the number of voxels L. The essential physics of this exponential decay is consistent with a mechanism in which a bound-ary voxel of the 1D chain flips to create an interface between high-and low-occupation voxels.Subsequently, each voxel at the interface flips until the whole chain has flipped.Letting p end denote the small rate of flipping a boundary voxel and p int denote the small rate of flipping a voxel at the interface, it follows that at short times P (B t |A 0 ) ∝ p end p L−1 int .If ∆t is the typical time to wait between interface flips, one therefore expects the exponential decay with k BA (L) ∝ e −L ln[pint ∆t ] .
The rate calculations as a function of L thus hint at a stepwise mechanism with a flipped voxel passing from one end of the chain to the other, but that mechanistic insight can also be more directly extracted from the DP/TDVP evolution.The time-evolution of |p t contains information about the entire distribution that flows from A to B, and we can track the average number of molecules in each lattice site as a function of time.For example, Fig. 6 shows n l (t) = 1| nl |p t , the average number of molecules in voxel l with l ranging from 1 through 4 for an L = 8 voxel lattice.At time t = 0, the plot shows the normalized steady-state confined to A, the high-occupation-number state.In the steady state, the six bulk voxels, 2 through 7, have essentially indistinguishable mean occupancies, while the boundary voxels, 1 and 8, have slightly depressed means reflecting that it is marginally more likely for those two voxels to fluctuate to low-occupancy states.As the time is allowed to evolve for a timescale t < τ mol , the mean occupation of the boundary voxels drops markedly due to the influence of rare trajectories in which one or both boundary voxels can flip.Subsequent to that drop in n 1 (t) , the mean occupancies of voxels 2 then 3 also drop, though less sharply, establishing the mechanistic sequence.After the molecular timescale τ mol , the rate of decrease in n l (t) becomes linear, reflecting the steady-state flux from A to B.
In addition to varying the extent of the space, it is interesting to tune the balance between reaction and diffusion events through the diffusive hopping rate d.The sequential mechanism suggests that a limitation to the transition events emerges from p int , the probability of flipping the state of an interior voxel at an interface between low-and high-occupation states.The transition events thus have low-and high-occupation blocks separated by an interface that diffuses until the low-occupation block covers the entire system.When molecules diffuse between interface voxels, the interface itself can diffuse, and the faster the interface sweeps through the system, the faster the transition events.One might therefore anticipate that the flipping rate grows with d.In fact, we find the exact opposite.Fig. 7 shows the results from varying the diffusion over 20 evenly spaced values while holding the bond-dimension fixed at κ = 20.The rate decreases for each value of L as d is increased due to diffusion acting as an effective stabilizer of high-occupation states.
At d = 0, each voxel is an independent well-mixed system.The switching rate depends on the probability that each voxel independently experiences a collapse of popu-  lation, dislodging it from its metastable high-occupation state.As d grows, there is a diffusive mechanism to fight against that population collapse.When the population of a voxel wanders dangerously close to a tipping point, molecules could diffuse from a high-occupancy neighbor to restore the population.In the limit of very fast diffusion, the entire chain of voxels effectively behaves as a single large well-mixed system with a markedly lower switching rate because statistical fluctuations away from the metastable state get correspondingly smaller.

VIII. DISCUSSION
Using the 1D Schlögl model, we have demonstrated how TN calculations combined with the DP framework offer a new paradigm for numerical investigation of reaction-diffusion models, even when the physical bond dimension (M + 1) is not small.That demonstration opens up a large class of RDME models which represent small well-mixed reactors, connected together by diffusive coupling.As illustrated in this work, these models can exhibit interesting physics when each reactor is small enough to exhibit significant statistical fluctuations and when the diffusion between reactors induces correlations between fluctuations in different reactors.In this work, our approach was to start with a well-studied kinetic toy model and to add spatial dynamics.Spatial extensions of other toy models like the Brusselator [76] or the Oregonator [77] would be amenable to similar numerical analysis because any reaction-diffusion model which can be written as a set of elementary reactions will have a DP representation.That representation can involve creation and annihilation operators for different particle types, not only for species X, but the mapping from reactions to effective Hamiltonians is systematic.Table I shows the effective Hamiltonians for several common (unchemostatted) reactions: unimolecular, bimolecular, and auto-catalytic.Given the set of reactions, it is comparatively straightforward to derive the DP representation, and computational tools already exist to generate an MPO [62] and propagate via TDVP from that DP form of the effective Hamiltonian.

Reaction
Hamiltonian Examples of the mapping between a common chemical reaction and the Doi-Peliti form of the contribution to an effective Hamiltonian.The procedure to derive such correspondences is systematic [78] and typically discussed for a Fock space with occupation number ranging from 0 to infinity.To truncate to a finite Fock space, one must additionally add a projector to each effective Hamiltonian, as in this work, to ensure conservation of probability.
As the TN techniques become increasingly well developed, one can imagine that the computational tools will offer generic "turn the crank" analysis of microscopic models of reaction-diffusion at the ensemble level.We anticipate these numerical approaches to be beneficial to both synthetic and biophysical reaction-diffusion systems.In the case of the former, the discretization of space into reactors can be explicitly realized, as in the microfluidic fabrication of a two-dimensional array of nanolitervolume Belousov-Zhabotinsky (BZ) reactors [79].In line with our 1D Schlögl model, the strength of coupling between those BZ nanoreactors is regulated by diffusion, and there is experimental interest in pattern formation as a function of that diffusive coupling.In the case of biophysics, sub-cellular dynamics is dominated by the generation and diffusion of proteins with small copy numbers.Much work has been devoted to understanding the logic of cellular decision making through stochastic CRN models, most famously gene regulatory models, which are most frequently studied in a well-mixed limit [80].We see these TN tools as an intriguing pathway to numerically investigate the timescales and pathways for switching between metastable states in well-mixed CRNs like gene toggle switches [10] and in reaction-diffusion variants of those gene-regulatory models.Eventually, one can dream of bringing these methods to the broad array of systems biology applications where RDME models are widely simulated [81], including in large-scale GPU-accelerated simulations of systems with biologically-relevant size [82].At present those sorts of reaction-diffusion model can only be numerically accessed through sampling approaches, and we expect the ensemble perspective will some day complement those efforts.
To move from the toy models to the synthetic and biophysical systems, it will be important to expand the ability of TNs to approximate correlation between species and voxels that are not simply connected through a 1D chain.Here, because our 1D chain of Schlögl reactors are connected only to their neighbors, it was natural to approximate the microscopic distribution with an MPS.That choice allowed us to capitalize on the efficient TDVP algorithm for MPS states, but a RDME with many dynamic species will not always be easily approximated by a 1D chain of tensors.If two chemical species are strongly correlated but far apart in the MPS structure, the MPS calculation will demand a very high bond dimension.It will be important for future work to dynamically adapt the bond dimension through subspace expansion [83] or by devising systematic ways to generate tractable TN architectures with connectivity mirroring the correlations between CRN species.Those technical advances may eventually lay the groundwork to move the methodology from this Schlögl study to the complex, biologically-motivated models where SSA realizations presently reign supreme.

IX. ACKNOWLEDGMENTS
We gratefully acknowledge Hadrien Vroylandt, Nils Strand, Geyao Gu, and Luis Pedro Garcia-Pintos for many insightful discussions.We are also grateful to Miles Stoudenmire, Matthew Fishman, Steven White, and other developers of ITensor, a library for implementing tensor network calculations, upon which this work was built.This research was supported in part through the computational resources and staff contributions provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and The reverse reaction, which transitions from an occupancy of n l + 1 to n l , has (n l + 1)n l (n l − 1)/6 different ways to combine three of the X molecules together.Consequently, the change in probability for the reverse reaction is Unlike the forward reaction, the negative term does not include a factor of δ n l ,M that treats the maximal occupancy state differently than the others.The distinction is that the loss terms from reaction 2 arise when probability is lost from occupancy n l down to occupancy n l −1.
Those loss terms are possible from n l = 3 through M .If there are fewer than 3 molecules, the reaction 2 mechanism is not possible, but this restriction emerges naturally from the n l (n l − 1)(n l − 2) coefficients in Eq. (A12).Expanding in the Fock basis therefore gives Mirroring the treatment of the forward reaction, we rewrite the first |n using a lowering operator Eqs. 14 imply x †2 x3 l |n c , n l + 1 = (n l −1)n l x l |n c , n l + 1 and x †3 x 3 l |n = n l (n l − 1)(n l − 2) |n , which allow us to simplify Eq. (A14) to We recognize the second sum over n as giving |p t .The first sum ranges from n l = 1 through M +1, but the M +1 term vanishes due to the occupancy cap.Furthermore, an n l = 0 term can be added to the sum since that term also vanishes (x †2 l x 3 l |n c , 0 l = 0).Thus the first sum can be re-indexed to give |p t .Performing both of those sums yields The contribution to the effective Hamiltonian from the reverse A-mediated reaction in voxel l is thus as and the operator-valued matrix associated with voxels 2 through L − 1 as With these definitions, one can show that Ĥ = W [1] • W [2] . . .• W [L−1] • W [L] , with • representing a contraction.
This decomposition of Ĥ offers a highly efficient way to numerically encode the effective Hamiltonian as an MPO.While we have shown an explicit MPO decomposition for the nearest-neighbor situation, in practice one uses powerful automated packages such as AutoMPO in ITensor [65] to generate compressed forms of the MPO from the second-quantized expressions like Eqs. ( 21), (22), and (23).
Appendix C: TDVP algorithm As described in Sec.IV, we wish to approximate the evolution of the distribution: Here |q t is a distribution over microstates that can be expressed as an MPS using the (time-dependent) set of tensors {Q n1 , Q n2 , . . ., Q n L }.As in the main text, the distinct Q tensors are implicitly distinguished based on the labeling scheme for the physical index.Only a subset of possible distributions |p t can be written in the form of |q t , but as the bond dimension κ increases, that subset necessarily grows.We call these distributions which can be represented by MPS states the MPS manifold M. Because Ĥ generally evolves |q t off M, the Time-Dependent Variational Principle (TDVP) [67,68] seeks the evolution of |q t along the MPS manifold: where PT M is a projector mapping the infinitesimal change in probability Ĥ |q t onto the MPS manifold's nearest (with respect to the l 2 norm) tangent vector.The projector in this equation of motion ensures that |q t stays on M for all times.Remarkably, in the case of the MPS manifold, an explicit form for the projector can be written down as a sum of 2L − 1 terms, L of which have a correspondence to the L tensors of the MPS and L − 1 of which have a correspondence to the bonds.To write down that decomposition, it is beneficial to first define the environment tensors which were obliquely referenced in Sec.IV.Through gauge transformations, |q t is expressed as where the A tensors are left-orthogonal and B tensors are right-orthogonal, meaning ni,si−1 and ni,si with physical indices counting the occupancy of sites l+1 through L and a single uncontracted bond index s l .
In terms of the left and right environment tensor projectors, we can now write out the 2L − 1 terms of the projector for the MPS manifold [72,85]: The positive terms P + l are associated to the L tensors and the negative terms P − l are associated to the L − 1 bonds.
Integrating Eq. (C3) for a time step δt gives In the limit of an infinitesimal time step, we can factorize that exponential to act with one l at a time.For example, The operator e δtP + 1 Ĥ has the physical interpretation of mapping site 1 forward in time by a time step δt.To see this interpretation, we illustrate the action of e δtP + 2 Ĥ on |q t , where we assume the MPS has already be written in a mixed canonical form centered on site 2. The operator e δtP + 2 Ĥ acts globally on the MPS state, but as long as |q t is expressed in the mixed canonical form, only the tensor at site 2 will be altered.In that way, the P + 2 term is implemented as a local update to a single tensor, a local update with the structure of a time evolution with a local effective Hamiltonian.The reduction from global to local operations is seen if the global exponential operator is expanded in a Taylor series: Each term of the sum over n can be expressed graphically, where the blue triangular MPS tensors indicate their left/right orthogonality by the direction they point toward the black diamond tensor at the center of the mixed canonical form: H( 2) By contracting over all of the tensors in the shaded purple region, we obtain the local tensor we call H(2).Notice that the same tensor appears repeatedly in the higher order terms of the Taylor series.These terms can be resummed to give back a local matrix exponential acting on site 2. Hence, by contracting Q n2 (t) (the black diamond) with e δtH (2) we get Q n2 (t + δt).
Following the P + 2 Ĥ propagation, the newly updated tensor Q n2 (t + δt) is then decomposed into two blocks R s 2 s2 via QR or SVD decomposition.Like how P + 2 Ĥ acted locally on the tensor at site 2, P − 2 will act locally to propagate R s 2 s2 , but the interpretation is that the propagation is a δt time step backwards in time, owing to the negative sign in Eq. (C11).R s 2 s2 is then contracted with the next tensor in the chain, B n3 s2s3 , which is also still at time t.In doing so, we shift gauge such that the new t + δt tensor assigned to site 2 is A n2 s1s 2 and the new center site tensor at site 3 becomes s2s3 .Having shifted the center to site 3, one then iterates: acting with P + 3 Ĥ, performing another QR or SVD decomposition, acting with P − 3 Ĥ, shifting the gauge, then continuing to sweep down the line.Sweeping over the entire MPS from n 1 to n L evolves |q t by δt.

Each product over Pout
s l at multiple voxels can be simplified noting that, When calculating observables using the ensemble method, sampling error is avoided, but there are several different sources for numerical error.Each source of error is controllable in that a parameter can be increased or decreased to systematically improve accuracy at the expense of a greater numerical cost.That trade-off between accuracy and cost is not straightforward, particularly because it can switch depending on which source of error dominates.The factorization of the time propagation in Eq. (C13) relied upon a infinitesimal time step, but the computational expense grows essentially linearly with the number of timesteps [86].Consequently, one seeks a time step which is small enough but not too small.We highlight that the SSA has an effective time step that emerges naturally as the typical time between subsequent events, and that this time step shrinks as the system size grows.By contrast, the TDVP approach's time step is a parameter one can choose irrespective of the system size.In addition to the time-step error, the propagation of the local tensors discussed in App C is implemented via a Krylov expansion.We found that 30 Krylov basis vectors were sufficient, but the model system here is not very sensitive to the number used.We have not presented a thorough analysis of time-step or Krylov errors because it is clear that our is dominated by the bond dimension κ, which can severely limit the variational search space of allowable MPS distributions.
In order to get an estimate of how the bond dimension impacts the numerical errors, we repeated the TDVP rate for varying κ and L. We would ideally compare each rate against an exact rate, but getting those exact rates is not numerically feasible for all system sizes considered.We have shown (see Fig. 5c) that for L = 2 through 6, the bond dimension κ = 20 was sufficient to reproduce the SSA rates to within the size of the plot markers.For the purposes of analyzing the bond-dimension rate kBA is calculated at different values of L and different bond-dimensions κ.The relative difference in rates is plotted against the "reference" rate calculated at κ = 20.for large enough bond dimension errors shrink exponentially.
convergence, we thus take that κ = 20 calculation as the "converged" reference value and measure the approach to that value as κ increases from 2 to 20.Fig. 8 shows the absolute fractional error comparing the rate at a given bond-dimension κ and the κ = 20 reference rate, all with the same time step and Krylov dimension.We see that errors can already be very small at modest bond dimension (where the smaller tensors make the calculations substantially faster), and the errors shrink exponentially with growing κ for all values of L.
Appendix F: Runtime Analysis Fig. 5c shows the total computational expense to estimate rates from the initial A ensemble.The dominant purpose of communicating the number of CPU hours is to illustrate the scaling of computational expense with the system size.That is to say, the brute force SSA rates grow exponentially more costly as L increases while the TDVP cost grows subexponentially.It is enticing to compare the costs to each other.At a fixed L, was the SSA calculation or the TDVP calculation cheaper?Fig. 5c appears to answer that question, showing that the two methods have a cross-over with TDVP becoming cheaper for sufficiently large L. We argue that the existence of the cross-over is a generic consequence of the different scalings with L, but the exact location of the cross-over depends on specific details of the implementations of both SSA and TDVP.Considering the difference between the κ = 20 (black) and κ = 30 (blue) lines in Fig. 5c, one must be especially careful to compare one method versus the other.Each method has parameters that can be adjusted to improve accuracy, and to compare different CPU hour timings, one should have a sense that the two calculations are returning comparable accuracies.Even more importantly, one should recognize that although we have used the methods to compute the rate, the computational expense goes into two very different objects.In the case of SSA, that expense results in a collection of representative trajectories, while for TDVP it yields time-dependent joint distributions.If one were using the calculations to compute different quantities (say waiting time distributions rather than just the rate), the same TDVP joint distributions that were generated to compute the rates could be used at essentially no additional computational expense.By contrast, generating a converged waiting time distribution from SSA samples would require many more than the 400 trajectories used to estimate the rate.In a similar vein, in Fig. 6 we showed how the TDVP calculations give converged insight into the mechanism of the rare event, but the SSA approach would return only 400 representative events from which one would extrapolate.
Having laid out the caveats, we spell out explicitly how we computed the total computational expense for these two classes of rate calculations (further information about forward flux sampling calculations follows in the next appendix).The TN calculations are not trivial to parallelize since the TDVP algorithm proceeds from one tensor to the next in serial, with the output of one tensor operation feeding into the next one.We did not seek a sophisticated parallelization, but we did note some imperfect parallelization could be easily accessed from the way the ITensor software in the Julia language handles the organization and contraction of tensors.Fig. 9 shows that this imperfect parallelization gave a roughly two-fold speed-up if 8 processors were used instead of a single processor.The computational expense reported in Fig. 5(c) is the total number of CPU hours for those 8 processors to complete the single-site TDVP dynamics with time step δt = 10 −4 for 10,000 time steps with 30 Krylov dimensions starting from an already computed distribution in A. In a sense, this timing is artificially worse than it would need to be.If we were willing to run for twice the duration of real time, we could decrease the total number of CPU hours for TDVP rates in Fig. 5(c) by a factor of 4 with a serial implementation.
Whereas the TDVP timing would have room for improvement with more sophisticated parallelization schemes, the SSA calculations are already fully optimized in that they are perfectly parallelizable.Each stochastic trajectory can and should be run on a separate CPU.The CPU hours reported in Fig. 5(c) is thus the sum of the wall times for 400 independent realizations, initialized in A and run until B was first reached.The cost for SSA calculations would shift up and down the log scale if one decided to use more or less than the 400 independent trajectories.Our choice here was to seek a number of trajectories that could make the SSA errorbars on the same order as the plot markers of Fig. 5(a).The wall time to complete 100 time steps of TDVP dynamics with L = 3 and κ = 30 as a function of the number of processors.Parallelization was not explicitly designed into the TDVP sweep, so the limited parallelism emerged only from the built in capabilities of ITensor.If one could reach perfect parallelization, either through a clever parallelization scheme or by just running TDVP in serial on a single processor, the number of CPU hours needed to compute the rates in Fig. 5 would decrease by up to a factor of four.

Appendix G: Comparison with Forward Flux Sampling
It is interesting to see how the computational expense of our rate calculation compares to advanced sampling rate calculations.Here we compare against the original FFS algorithm [10].The key idea of FFS is that a (typically one-dimensional) progress coordinate is defined to decompose the single rare event from A to B down into less rare excursions that make partial progress along that order parameter.The progress coordinate is discretized to define a collection of N non-overlapping interfaces λ 1 , λ 2 , . . ., λ N .The FFS procedure decomposes the transition rate as where Φ 1,A is the flux to pass from A to the λ 1 interface and P (λ k+1 |λ k ) is the probability of reaching λ k+1 before first returning to A conditioned upon initialization at the λ k interface.To the extent that the chosen progress coordinate acts as a good reaction coordinate, the method can radically improve the efficiency of a rate calculation compared to brute force SSA, but if the order parameter is not "close enough" to the reaction coordinate, the conditional probability terms can become very expensive to compute [87].It is tempting to think that P (λ k+1 |λ k ) requires a trajectory starting at λ k and run until it hits either the next or the previous interface, but the actual procedure requires the trajectory to hit either the next interface or A. When the progress coordinate (a) The rate kBA from SSA and TDVP using κ = 30 bond-dimension compared to FFS using different order parameters.Though FFS can correctly estimate the rate even when the order parameter is not exactly the reaction coordinate, the choice of order parameter significantly impacts the computational expense.It is notable that the TDVP rate calculation can bypass the introduction of an order parameter; in fact, by analyzing the TDVP-evolved joint distribution, one can reveal the reaction mechanism as an output of the calculation.
poorly approximates the reaction coordinate, a tremendous amount of time can be wasted simulating trajectories slowly wandering back to A from one of the final interfaces.The efficiency of FFS calculations additionally depends on the number of interfaces, their distance from each other, and the number of needed to resolve the conditional distributions.Given all the tunable knobs that go into a FFS rate calculation, it only makes sense to compare the numerical expense of TDVP rate calculations against a prudently designed FFS with a reasonable choice of progress coordinate.
Without a priori knowing the mechanism of the Schlögl reaction-diffusion switching events, we designed two different order parameters and performed calculations with both.For the first order parameter, we focus on the global number of X molecules, summed over all voxels.The interfaces for this order parameter measuring global progress are λ G ≡ A, λ G 1 , . . ., λ G 8 , B with λ G i = Lλ * A +iL. Remember that λ * A = 15 and λ * B = 25, so this order parameter tracks the progress of the whole sys-tem in moving from less than 15L molecules to more than 25L molecules.The second order parameter, inspired by Fig. 6's demonstration that the voxels on the end decrease first, uses only the number of molecules in the first voxel as the sign of progress.The interfaces for this order parameter measuring voxel 1 are λ 1 ≡ A, λ 1 1 , . . ., λ 1 8 , B with each λ 1 i interface recording the addition of another X molecule in voxel 1 (λ 1 i = 15 + i), so that λ 1 1 = 16 and λ 1 8 = 23.
There are many ways to optimize both FFS calculations, for example, changing the number of interfaces, the number of trajectories to estimate Φ 1,A , and the number of trajectories to estimate the conditional probabilities.These changes impact the computational expense, and it is not generally knowable how to tune those parameters to most cheaply decrease the errorbars of the rate calculation.We found that 4,000 trajectories were sufficient to get decent estimates of the flux Φ 1,A for both order parameters, though the expense of computing that flux was much greater for λ G than for λ 1 .While the flux was more costly to compute for λ G , the conditional probabilities were easier.Each conditional probability could be well estimated with only 4,000 additional trajectories per interface.Converging those conditional probability terms for λ 1 required 10 7 trajectories per interface, a cost that was at least partially offset by the cheaper Φ 1,A calculation.For each order parameter we generated these samples to calculate ten independent estimates of the rate, allowing us to measure both the mean and standard error.From these standard errors, we determined that the FFS calculations were converged in line with the SSA errors in that plot.
A notable feature of FFS is that it will return the correct transition rates even when the order parameter is not a perfect reaction coordinate.Fig. 10a reflects that in that both λ G and λ 1 reproduce the rates given by the SSA and TDVP methods.However, the computational expense of the two are not identical.10b plots the CPU hours needed to generate the samples that yield the means and standard errors of Fig. 10a.At small L, λ 1 is faster; at large L, λ G is faster.The FFS calculations have many tunable parameters.They could likely be improved through further refinement of order parameters, changing the number of interfaces, fine-tuning the number of samples for each conditional probability, etc.Still, we find it notable that our reasonable but perhaps unoptimized FFS implementations perform with similar (and slightly greater) expense as the TDVP rate calculations.This observation is particularly notable since the TDVP calculations do not need one to guess an order parameter and since they offer converged statistical information about the mechanism in addition to just giving the rate.

FIG. 1 .
FIG. 1.(a) A typical trajectory in the bistable regime (c1 = 36.38,c2 = 0.040, c3 = 2200, and c4 = 37.6) of the wellmixed Schlögl model of Sec.II A. The system stochastically switches between the two metastable basins, where basin A is the top colored region and B is the bottom colored region.(b) Starting with a distribution localized in basin A, the evolution has three time scales.Initially, no trajectories can reach B.After the so-called molecular timescale τ mol , the fastest trajectories begin to reach B. Subsequently, probability leaches into region B at a constant rate with P (Bt|A0) ∝ kBA.The inset shows that growth on a much longer timescale, revealing the third scale when the system eventual approaches the steady state.
FIG.3.The ensemble evolution of the many-body state requires that one solve for the probability of every microstate at all times.That high-dimensional joint probability distribution can be well approximated by a matrix product state (MPS).An MPS approximation for the joint distribution can be advantageous whether there are Ns species in a well-mixed reaction [(a) top], one dynamical (unchemostatted) species with diffusion [(a) middle], or multiple species and diffusion [(a) bottom].In all three cases, each circle represents a tensor, either a rank-2 at the first and last site or rank-3 in the body.Horizontal lines represent the bond-indices connecting tensors and vertical lines are physical indices.For an RDME, the physical indices represent the possible number of molecules at a given tensor.If two lines are connected, then the index is summed over.(b) An effective Hamiltonian H can be likewise decomposed into lower-rank tensors to form a matrix product operator (MPO).Now there are two physical indices per tensor, representing transitions from state n X l to m X l .Vertical connected lines are the bond indexes which are contracted over.

75 FIG. 4 .
FIG. 4.Snapshots of the joint density in two neighboring voxels from an L = 3 site lattice of Schlögl voxels.Using c1 = 2.676, c2 = 0.040, c3 = 108.102,c4 = 37.881 and d = 8.2207, the steady-state distribution was computed with κ = 30 by evolving an initial uniform distribution according to Ĥ via TDVP.That steady-state distribution was projected into A at time 0, renormalized, then propagated according to the modified no-recrossing dynamics H using TDVP.For clarity of visualization, the density within A is not plotted as it would overwhelm the small probability that leaks into B. By t = 0.01 (left image), probability has begun to leak from A but none has reached B. Probability has just begun to reach B at at t = 0.14, roughly the τ mol timescale for this problem.By t = 0.44, probability is collecting in B, and by t = 0.75, the flux passing from A to B has attained an approximately constant steady-state value corresponding to the linear growth phase where P (Bt|A) ∝ kBAt.
FIG. 5. (a)The rate kBA for switching between metastable states of the L-voxel Schlögl model is calculated using the SSA (red line) and using TDVP (dashed lines) with varying bond-dimensions κ using the same rate parameters as Fig.4.For sufficiently large bond-dimensions, TDVP agrees with a SSA rate estimate built using 400 transitions from A to B at each L. (b) As L grows a larger bond-dimension is required to capture the dynamics.Signatures of numerical instability are apparent when κ is too small, shown here for L = 3, which stably converges when κ exceeds 6.For this system, κ = 20 is sufficient for numerical stability up to L = 6, while κ = 30 gives stable solutions up to L = 8.(c) The total number of CPU hours to estimate kBA grows exponentially with SSA but polynomially with TDVP.

46 FIG. 6 .
FIG. 6.The time evolution of the average number of molecules for different lattice sites when L = 8 and rate parameters are as in Fig.4, computed with κ = 30.Different colors correspond to different voxels.The expected number of molecules n l (t) = 1| nl |pt for voxels 1 to L/2 (shown) are mirrored by voxels L/2 + 1 to L. The initial non-linear flow of probability out of A is shown by the steep drop in average number of molecules.The system passing τ mol is shown by the transition to a linear change in particle number n l (t) .

3 FIG. 7 .
FIG.7.The switching rate kBA as a function of diffusion coefficient d for 2 ≤ L ≤ 6, computed with κ = 20 for the rates from Fig.4.As d is increased, the diffusion acts to inhibit the stochastic switching of each voxel.At large enough d the system acts similar to a large well-mixed system, resulting in a leveling off of the rate coefficients.

) 1 =
Summing the A tensors found in Eq. (C4) up to l − 1 gives the left environment tensorΦ [n1:n l−1 ] L,s l−n1,...,n l−1 s1,...,s l−2 A n1 s1 A n2 s1,s2 . . .A n l−1 s l−2 s l−1 |n 1 ,n 2 , . . .n l−1 (C7) with physical indices counting the occupancy of sites 1 through l − 1 and a single uncontracted bond index s l−1 .A similar sum over the s indices of the B tensors gives the right environment tensor Φ [n l+1 :n L ] R,s l = n l+1 ,...,n L s l+1 ,...,s L−1 D4)    where the sum over c k is now all combinations of N and Q where Q appears k times.The advantage of Eq. (D4) is that effective Hamiltonians can be projected onto only using local (single-voxel) projectors.In this work, B consists of all molecule numbers less than q * B + 1.Then, projecting onto reactions that can map out of B the effective Hamiltonians for the Schlögl model become Hrxn = (I − Pout ).(D6)Appendix E: Error in TDVP rate calculations FIG. 8.rate kBA is calculated at different values of L and different bond-dimensions κ.The relative difference in rates is plotted against the "reference" rate calculated at κ = 20.for large enough bond dimension errors shrink exponentially.
FIG.9.The wall time to complete 100 time steps of TDVP dynamics with L = 3 and κ = 30 as a function of the number of processors.Parallelization was not explicitly designed into the TDVP sweep, so the limited parallelism emerged only from the built in capabilities of ITensor.If one could reach perfect parallelization, either through a clever parallelization scheme or by just running TDVP in serial on a single processor, the number of CPU hours needed to compute the rates in Fig.5would decrease by up to a factor of four.
FIG. 10.(a)The rate kBA from SSA and TDVP using κ = 30 bond-dimension compared to FFS using different order parameters.Though FFS can correctly estimate the rate even when the order parameter is not exactly the reaction coordinate, the choice of order parameter significantly impacts the computational expense.It is notable that the TDVP rate calculation can bypass the introduction of an order parameter; in fact, by analyzing the TDVP-evolved joint distribution, one can reveal the reaction mechanism as an output of the calculation.