\section{Introduction}

Connecting microscopic kinetics to emergent rates has a long history for equilibrium and nonequilibrium systems.
Complex systems often settle into metastable basins with rare transitions regulating the rates of switching between the emergent states.
In equilibrium, transition state theory~\cite{eyring1935activated, laidler1983development}, variational transition state theory~\cite{truhlar1980variational}, Grote-Hynes theory~\cite{grote1980stable,hynes1985chemical}, and transition path sampling (TPS)~\cite{dellago1998transition,bolhuis2002transition} 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 still extract rates from microscopic dynamics, either through a noise-guided TPS approach~\cite{crooks2001efficient, gingrich2015preserving} or using forward flux sampling (FFS)~\cite{allen2005sampling,allen2006simulating,allen2009forward}.
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 \cite{ma2005automatic}?
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\textemdash 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, continuous-time, 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 variety 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 \cite{allen2005sampling}.
More specifically, we study discrete reaction-diffusion dynamics of \(L\) sites or voxels, \cite{isaacson2006incorporating,erban2007practical}.
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 \cite{fange2010stochastic,smith2019spatial}. 
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 \cite{craciun2006understanding,krause2021modern,kondo2010reaction}.
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~\cite{turing1952chemical, cross1993pattern}.
It is, however, well-appreciated that stochastic effects can quantitatively and even qualitatively impact kinetics~\cite{mcquarrie1967stochastic,lee2006non,mcadams1997stochastic,thattai2001intrinsic,vilar2002mechanisms,paulsson2004summing,lee2006non}.
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~\cite{tian2006stochastic,allen2006simulating}.

Historically, similar schemes for the CME ensemble evolution have been hampered by the curse of dimensionality. 
As a result, exact solutions are rare \cite{heuett2006grand} and numerical solutions often resort to sampling methods built upon the kinetic Monte-Carlo stochastic simulation algorithm (SSA) \cite{gillespie2013perspective,cao2006efficient}, also known as the Gillespie algorithm. 
Another approach is to evolve a subset of the state space through the finite state projection (FSP) method~\cite{munsky2006finite}. 
FSP has previously been combined with TNs to approximately describe the CME \cite{kazeev2014direct,vo2017adaptive}, 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~\cite{strand2022using,strand2022computing}.
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\"ogl model~\cite{schlogl1972chemical} in the bistable regime, both with and without diffusion. 
In Sec.~\ref{sec:wellmixed}, 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.~\ref{sec:dp}, we review the Doi-Peliti framework for writing reaction-diffusion dynamics in a second-quantized field theoretic form.
Sec.~\ref{sec:tn} then reviews the TN matrix product state (MPS) ansatz that enables ensemble evolution via TDVP.
Those methods are applied to the many-body rate calculation in Sec.~\ref{sec:mbrates}.
Finally, in Sec.~\ref{sec:Results}, we show results for the RDME with diffusion between \(L\) well-mixed Schl\"ogl 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.

\section{Ensemble rate calculations}\label{sect:EstRate}
\label{sec:wellmixed}

Consider an archetypal two-state rate process between reactants and products $\ce{R <=> P}$, initiated entirely in the reactant state.
After some reaction timescale, the concentration of products grows before eventually leveling off at its equilibrium 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 trajectories, 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.

\subsection{Well-mixed Schl\"ogl Model}
\label{sec:WMSchlogl}
Consider the well-mixed Schl\"ogl model, a quintessential example of bistability in chemical reaction networks (CRNs)~\cite{vellela2009stochastic}.
The model tends to be studied under chemostatted conditions that fix the concentrations of species A and B while allowing species X's concentration to evolve according to the four reactions:
\begin{align}
\ce{$2$X + A &<=>[\tilde{c}_1][\tilde{c}_2] $3$X}\label{rxn:S1} \\ 
\ce{B &<=>[\tilde{c}_3][\tilde{c}_4] X}
\label{rxn:S2}
\end{align}
where $\tilde{c}_1,\tilde{c}_2,\tilde{c}_3$ and $\tilde{c}_4$ are the stochastic rate coefficients.
Our starting point is microscopic such that the rate \(\tilde{c}_1\) 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\) for each reaction \(r\). The relation between $k_r$ and the microscopic rate coefficients \(\tilde{c}_r\) are only equivalent for linear reactions and in the thermodynamic limit~\cite{higham2008modeling}.
By varying the number of A and B molecules along with the rate coefficients, the Schl\"{o}gl system can occupy either a unimodal or a bimodal regime.
In the bimodal regime, the number $n$ of species X fluctuates, switching between two metastable states.
A representative trajectory, shown in Fig.~\ref{fig:ER}, oscillates between a state with \(n \approx 5\) and another with \(n \approx 50\) molecules.
We first set out to compute the rate of bistable switching between those two metastable states, making use of the fact that the dynamics is one-dimensional (1D)\textemdash \(n\) executes a random walk between 0 and some large maximum number \(M\) that we impose for the sake of calculation.

\begin{figure}[htb]
\centering
\includegraphics[width=.45\textwidth]{Figures/Fig1.pdf}
\caption{\label{fig:ER}(Top) A typical trajectory in the bistable regime ($c_1 = 2.676, c_2 = 0.040, c_3 = 108.102,$ and $c_4 = 37.881$) of the well-mixed Schl\"{o}gl model of Sec.~\ref{sec:WMSchlogl}.
  The system stochastically switches between the two metastable basins, where basin $\mathcal{A}$ is the top colored region and $\mathcal{B}$ is the bottom colored region.
  (Bottom) Starting with a distribution localized in basin $\mathcal{A}$, the transition rate can be extracted from the probability growth in the linear regime of the evolution, $P(\mathcal{B}_t|\mathcal{A}_0) \propto k_{\mathcal{B}\mathcal{A}}t$.
  The inset shows that growth on a much longer timescale, revealing the eventual approach toward equilibrium.}
\end{figure}

Rather than track \(n(t)\) for a single trajectory as in Fig.~\ref{fig:ER}, 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: \(\ket{p_t} = \left[p_t(0), p_t(1) \hdots, p_t(M)\right]^T\), where \(T\) indicates the transpose to give a column vector.
The distribution evolves according to the master equation
\begin{equation}
\frac{d\ket{p_t}}{dt} = H\ket{p_t}.
\label{eq:ClSch}
\end{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 $\bra{\mathbf{1}}H = 0\bra{\mathbf{1}}$, with $\bra{\mathbf{1}}$ denoting a vector of ones.
The distribution \(\ket{p_t}\) can always be formally evolved in terms of a propagator \(e^{Ht}\) as 
\begin{equation}
\ket{p_t} = e^{Ht}\ket{p_0}.
\label{eq:Propagator}
\end{equation}
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 $\ket{\boldsymbol{\pi}} = \lim_{t\to \infty} e^{Ht}\ket{p_0}$. 

\(H\) can be decomposed into contributions from the two different reversible mechanisms of Eqs.~\eqref{rxn:S1} and~\eqref{rxn:S2}, labeled \(H^{A}\) and \(H^{B}\), respectively, based on whether the transitions are mediated by the chemostatted species A or B.
Then, the Schl\"ogl mechanism at the level of the master equations corresponds to tridiagonal rate matrices with elements
\begin{align}
  \label{eq:tridiagonal}
H^A_{n'n} = &\frac{\tilde{c}_1 n (n-1)}{2}\delta_{n'-1,n} \frac{\tilde{c}_2 n (n-1)(n-2)}{6}\delta_{n'+1,n} - \nonumber \\
&\frac{n (n-1)}{2}\left(\tilde{c}_1 + \frac{\tilde{c}_2(n-2)}{3}\right)\delta_{n',n}
\end{align}
and
\begin{equation}
H^B_{n',n} = \tilde{c}_3\,\delta_{n'-1,n} + \tilde{c}_4n \delta_{n'+1,n} - \left(\tilde{c}_3 + \tilde{c}_4n \right)\delta_{n',n},
\end{equation}
where $H = H^A + H^B$.
Given the specific tridiagonal \(H\) and the restriction \(0 \leq n \leq M\), both \(H\) and the propagator \(e^{Ht}\) are \((M+1) \times (M+1)\) matrices, amenable to evolving the entire distribution of realizations, 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 \(\mathcal{A}\) and \(\mathcal{B}\) (red and blue regions of Fig.~\ref{fig:ER} respectively) with \(\mathcal{A} \cap \mathcal{B} = \emptyset\).
We write \(\hat{\mathcal{P}}_{\omega}\) for the non-negative idempotent ($\hat{\mathcal{P}}_{\omega}^2 = \hat{\mathcal{P}}_\omega$) operator that projects onto the density lying in a region \(\omega\), where \(\omega \in \{\mathcal{A},\mathcal{B}\}\).
If \(\ket{p_t}\) has support over a state space \(\Omega\), then the projection \(\hat{\mathcal{P}}_\omega \ket{p_t}\) is an unnormalized distribution with support restricted to \(\omega\).
The total probability in \(\omega\) at time \(t\) can then be expressed as \(\bra{\mathbf{1}}\hat{\mathcal{P}}_\omega\ket{p_t}\).

\subsection{Rate Constants}
\label{sec:rateconstants}
In general, basins can be arbitrary subsets of $\Omega$.
For the current example, $\mathcal{A}$ and $\mathcal{B}$ are basins of attraction for the Schl\"{o}gl model (see Fig.~\ref{fig:ER}a), and we would like to measure the rate of going from $\mathcal{A} \rightarrow \mathcal{B}$.
Rate constants can be defined either for the unidirectional flux from $\mathcal{A} \rightarrow \mathcal{B}$, the quantity we focus on in this work, or for the timescale of the relaxation to equilibrium, which also involves the $\mathcal{B} \rightarrow \mathcal{A}$ transitions~\cite{roux2022transition}.
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 \(\mathcal{A}\) to the boundary of \(\mathcal{B}\) yields the Transition Path Theory (TPT) rate related to the committor function~\cite{e2010transition}.
Here, we will have access to the distribution over all states in the interior of $\mathcal{A}$ and $\mathcal{B}$, allowing for the calculation of the unidirectional rate \(k_{\mathcal{B} \mathcal{A}}\).
That $k_{\mathcal{B}\mathcal{A}}$ is the rate to start in $\mathcal{A}$ at time zero with probability $P(\mathcal{A}_0)$, and transition to any state in $\mathcal{B}$ at time $t$ with conditional probability $P(\mathcal{B}_t|\mathcal{A}_0)$.
We initialize with a distribution over microstates $\hat{\mathcal{P}}_\mathcal{A}\ket{\boldsymbol{\pi}}/\bra{\mathbf{1}}\hat{\mathcal{P}}_{\mathcal{A}}\ket{\boldsymbol{\pi}}$, formed by projecting the equilibrium distribution onto $\mathcal{A}$ and normalizing, making $P(\mathcal{A}_0)=1$ by construction.
With this initialization, the rate for transitioning from $\mathcal{A}$ to $\mathcal{B}$ can then be written in terms of the time derivative of the conditional probability of occupying $\mathcal{B}$ a time $t$ later \cite{allen2005sampling}:
\begin{align}
k_{\mathcal{B} \mathcal{A}} &\equiv \frac{d}{dt} P(\mathcal{B}_t | \mathcal{A}_0)P(\mathcal{A}_0)\biggr\rvert_{t >\tau_{\rm mol}},\nonumber \\
 &=  \frac{d}{dt} \frac{\bra{\mathbf{1}}\hat{\mathcal{P}}_{\mathcal{B}}\,e^{Ht}\, \hat{\mathcal{P}}_{\mathcal{A}} \ket{\boldsymbol{\pi}}}{\bra{\mathbf{1}}\hat{\mathcal{P}}_\mathcal{A}\ket{\boldsymbol{\pi}}}\biggr\rvert_{t >\tau_{\rm mol}}.
    \label{eq:rate}
\end{align}
That derivative is evaluated at some timescale exceeding $\tau_{\rm mol}$, a molecular time scale of the system that is long enough for the first trajectories to reach \(\mathcal{B}\) from \(\mathcal{A}\) but much shorter than the typical transit time $k_{\mathcal{B}\mathcal{A}}^{-1}$.
To the extent that there is a separation of timescales between the time to explore the basin and the time to transit between basins, distinctions between this rate and a TPT rate are insignificant.
For rare events, no transitions can occur faster than the molecular time scale \(\tau_{\rm mol}\), beyond that point, the number of transitions accumulates over time as \(P(\mathcal{B}_t|\mathcal{A}_{0}) \propto k_{\mathcal{B}\mathcal{A}}  t\).
At times long compared to \(k_{\mathcal{B} \mathcal{A}}^{-1}\), \(P(\mathcal{B}_t|\mathcal{A}_{0})\) must eventually stop growing linearly to level off at the stationary value \(\bra{\mathbf{1}}\hat{\mathcal{P}}_\mathcal{B} \ket{\boldsymbol{\pi}}\), but our interest is in extracting \(k_{\mathcal{B} \mathcal{A}}\) from the growth at much shorter times.
A typical realization of the process will not require observation much beyond \(\tau_{\rm mol}\) units of time, to ascertain the rate, meaning our ensemble evolutions must be propagated for roughly $\tau_{\rm mol}$.

Our ability to compute \(k_{\mathcal{B} \mathcal{A}}\) 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.~\ref{eq:rate} generalizes naturally for harder problems (those with more than one species or with diffusion between multiple well-mixed voxels), the state space for \(\ket{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.
 
\begin{figure*}[!ht]
\centering
\includegraphics[width=0.95\textwidth]{Figures/Fig2.pdf}
\caption{\label{fig:one} The reaction-diffusion master equation (RDME) consists of voxels in which well-mixed chemical reactions occur (a) and between which species can diffuse (b).
  This work centers around a one-dimensional chain of $L$ such voxels, each with a fluctuating number of molecules of species X in a Schl\"{o}gl chemical reaction.
  The Doi-Peliti (DP) framework turns the reaction-diffusion evolution into a second-quantized effective Hamiltonian that governs the evolution of the distribution over classical reaction-diffusion microstates.
  (c) That effective Hamiltonian is represented by a high-dimensional tensor, with each vertical leg taking one of $(M+1)$ possible values corresponding to the number of molecules of species X which could occupy each voxel.
  The DP form of the effective Hamiltonian lends itself to a compact representation of the effective Hamiltonian tensor of rates for transitioning from one microstate to another.}
\end{figure*}

\section{Doi-Peliti}

\label{sec:dp}
The Doi-Peliti (DP) formalism is a classical version of the second quantization methods from quantum field theory \cite{ohkubo2013algebraic}.
Originally introduced by Doi \cite{doi1976stochastic} as a critique to the assumption of the Smoluchowski equation, it was later developed by Peliti \cite{peliti1985path} and others \cite{rey1999asymptotic} as  a path integral formalism for solving chemical reaction networks.
The DP framework has also found success in other areas such as studying kinetically constrained lattice models in the context of glassy dynamics~\cite{garrahan2007dynamical,garrahan2009first}.
The DP representation builds operators which encode the dynamics onto the high-dimensional space that accompanies many-body problems.
Crucially, writing the chemical mechanism in an operator form guarantees it can be easily converted into a matrix product operator (MPO)~\cite{schollwock2011density,causer2022finite}, as we will do in Sec.~\ref{sec:tn}.\\

To start, we go beyond Sec.~\ref{sec:wellmixed} by considering well-mixed reactions with multiple dynamic species  as well as heterogeneous systems where molecules can diffuse between neighboring voxels.
Let $X_1, X_2,\ldots X_L$ represent each chemical species in each voxel.
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,
\begin{equation}
\ce{ $\sum_j \eta_j^r$ X$_j$ ->[\tilde{c}_r]  $\sum_j \mu_j^r$ X$_j$}
\label{eq:rxn}
\end{equation}
where $\tilde{c}_r$ are again the rate coefficients for reaction $r$.
$\eta_j^r$ and $\mu_j^r$ reflect the stochiometry of the reactions and give the stochiometric vector $\nu^r = \mu^r - \eta^r$ which evolves the system from state $\mathbf{n} \in \Omega$ in state space $\Omega$, to $\mathbf{n}' = \mathbf{n} + \nu^r$.
In the previous section, these reactions induced changes between the \(M+1\) different microstates, specified by the scalar \(n\).
Now, the microstates $\mathbf{n}$ are specified by \(\mathbf{n} = [n_1, n_2, \hdots, n_L]\), counting the occupation of each species within each voxel. 

For a dilute well-mixed solution, where well-mixed means molecules are evenly distributed around a given sub-volume, the change in probability  of microstate $\mathbf{n}$ follows the gain-loss equation \cite{isaacson2006incorporating},
\begin{equation}
\frac{dp_t(\mathbf{n})}{dt} = \sum_{r}^N\left[ \alpha_{\mathbf{n} - \nu^r}^r p_t({\mathbf{n} - \nu^r}) - \alpha_{\mathbf{n}}^r p_t({\mathbf{n}})\right],
\label{eq:chgPr}
\end{equation}
where the sum is taken over \(N\) different reactions of the form Eq.~\eqref{eq:rxn}, each with their own stochiometric vector \(\nu^r\).
The transition rate between any two states is the sum of contributions from each reaction \(H_{\mathbf{n}'\mathbf{n}} = \sum_r \alpha^r_{\mathbf{n}}\).
Each propensity $\alpha^r$ is a product of a stochastic rate coefficient $\tilde{c}_r$ and the combinatorial number of ways to select either $\eta_j$ or $\mu^r$ molecules from reaction $r$, depending on the direction of the reaction.
For a ``forward" reaction the propensity is, $\alpha_r(\mathbf{n}) = \tilde{c}_r\prod_{n_j\in r}{n_j\choose\eta_j}$.

Encoding the reactions into a matrix form allows us to evolve the elements of the joint distribution $p_t(\mathbf{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(\mathbf{n})\).

Since particles of the same species in the same voxel are indistinguishable, it is beneficial to work in an occupation number representation with
\begin{equation}
    \ket{p_t} = \sum_{\mathbf{n}} p_t(\mathbf{n}) \ket{\mathbf{n}}.
\end{equation}
The basis vectors
\begin{align}
    \nonumber \ket{\mathbf{n}} &= \ket{n_1, n_2, \hdots n_L}\\
    &= \ket{n_1} \otimes \ket{n_2} \otimes \hdots \otimes \ket{n_L}
\end{align}
corresponds to a possible microstate, and the set of basis vectors forms a tensor product space called a Fock space $\mathcal{F}^L$ with $\ket{\mathbf{n}} \in \mathcal{F}^L$ and \(\bra{\mathbf{m}}\ket{\mathbf{n}} = \delta_{\mathbf{m}\mathbf{n}}\)~\cite{del2022probabilistic}.
Without loss of generality, we will discuss the case of a single chemical species and \(L\) physical sites, in which case \(n_l\) has the interpretation of the number of molecules at site \(l\).
The second-quantization procedure rewrites \(H\) in terms of creation and annihilation operators, \(x_l^\dagger\) and \(x_l\), that respectively raise and lower the number of X molecules at site \(l\).
That is to say,
\begin{align}
x_l^\dagger\ket{n_l} &= \ket{n_l+1},\nonumber \\
x_l\ket{n_l} &= n_l\ket{n_l-1},\nonumber \\
x_l\ket{0} &= 0.
\label{eq:raiselower}
\end{align}
where $\ket{0}$ is the vacuum state, i.e., no molecules.
By $x_l$ or $x_l^\dagger$ we are implying that the operator only acts on the $l^{\rm th}$ site:
\begin{equation}
  x_l = \mathbb{I}_1 \otimes \mathbb{I}_2\otimes\ldots\otimes \mathbb{I}_{l-1}\otimes  x \otimes \mathbb{I}_{l+1} \otimes \ldots \otimes \mathbb{I}_L
\end{equation}
Because they act on the single-site space, the general one-site operators $x$ and $x^\dagger$ can be written as the \(M+1 \times M+1\) matrices
\begin{equation}
x^\dagger = \begin{bmatrix} 
0 && 0 && 0 && \ldots \\
1 && 0 && 0 && \ldots \\
0 && 1 && 0 && \\
\vdots && && \vdots && \ddots
\end{bmatrix} \text{ and }
\qquad
x = \begin{bmatrix} 
0 && 1 && 0 && \ldots \\
0 && 0 && 2 && \ldots \\
0 && 0 && 0 && \\
\vdots && && \vdots && \ddots
\end{bmatrix}.
\end{equation}
That these matrices yield Eq.~\eqref{eq:raiselower} can be readily confirmed.
More abstractly, the existence of such raising and lowering operators in this classical setting originates from the classical commutation relation \([x_i x_j^\dagger] = \delta_{ij}\) and the existence of a positive semi-definite number operator $\hat{n}_l = x_l^\dagger x_l$, so named because \(\hat{n}_l\ket{n_l} = n_l \ket{n_l}\).
Together, the set $\{\Omega,\mathcal{F}^L,x,x^\dagger\}$ are also sometimes called an interacting Fock space \cite{ohkubo2013algebraic}.
Our use of the raising and lowering operators is less formal and more pragmatic.
As in the construction of the number operator, we will be able to decompose contributions to \(\hat{H}\) in terms of the \(x^\dagger\) and \(x\)'s.

\begin{figure*}[!ht]
\centering
\includegraphics[width=.9\textwidth]{Figures/Fig3.pdf}
\caption{\label{fig:mps} 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 $N_s$ 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_l^X$ to $m_l^X$.
  Vertical connected lines are the bond indexes which are contracted over.}
\end{figure*}

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, \(\hat n = \sum_l \hat n_l\).
The operator \(\hat{n}\) acts on a many-body state \(\ket{\mathbf{n}}\), yet each \(\hat{n_l}\) is a single-body operator acting on site $l$ and leaving all other sites unchanged.
It is therefore understood that when local operators such as $\hat n_l$ act on $\ket{\mathbf{n}}$, it is written as a shorthand for
\begin{equation}
\hat n_l = \mathbb{I} \otimes \hdots \otimes x^\dagger_lx_l \otimes \hdots \otimes \mathbb{I}.
\end{equation}
In terms of that shorthand, $\hat n_l \ket{\mathbf{n}} = n_l \ket{\mathbf{n}}$ counts the number of particles in site $l$ for the many-body state $\ket{\mathbf{n}}$.
The time-dependent expected number of particles in the whole system is thus \(\langle n(t) \rangle = \sum_l\bra{\mathbf{1}}\hat n_l \ket{p_t}\), where the vector of ones \(\bra{\mathbf{1}}\) serves to sum over all microstates.
The time-dependence of a distribution can also be extracted from the action of a second-quantized \(\hat{H}\) on \(\ket{p_t}\), and this operation can be practically computed when \(\hat{H}\) is expressed in terms of the single-site operators.

Mirroring the discussion in Sec.~\ref{sec:wellmixed}, the generator $\hat H_l$ for the Schl\"ogl model reactions within the  $l^{\rm th}$ voxel is a sum of contributions from the elementary reactions $\hat H_l = \sum_r \hat H^r$.
Fig.~\ref{fig:one}a shows each reaction, expressed as a second-quantized contribution to the effective Hamiltonian, as derived in App.~\ref{app:DPEnc}.
Combining together the forward and reverse reactions mediated by species A, we see that the reversible reaction Eq.~\eqref{rxn:S1} contributes 
\begin{equation}
   \hat{H}^A_l = x_l^{\dagger 2}\left[\left(c_1x^\dagger_l + c_2\left(\mathbb{I} - x_l^\dagger\right)x_l\right)x_l^2 - c_1x_l^2\hat{\mathcal{P}}^c_{M_l}\right].
  \label{eq:HA}
\end{equation}
to the effective Hamiltonian. For convenience we have absorbed the combinatorial terms into that rate coefficient, defining $c_r \equiv \tilde{c}_r / \gamma!$, where $\gamma$ is either $\eta^r$ or $\mu^r$ for forward and reversed reactions, respectively. Likewise, the contribution
\begin{equation}
  \hat{H}_l^B = x_l^\dagger(c_3-c_4)x_l + c_4x_l - c_3\hat{\mathcal{P}}^c_{M_l}.
  \label{eq:HB}
\end{equation}
comes from Eq.~\eqref{rxn:S2}.
We show in App.~\ref{app:DPEnc} that the truncation of the state space at finite \(M\) requires the addition of a projector, $\hat{\mathcal{P}}^c_{M_l}$ to prevent \(\hat{H}_l\) from promoting any probability from the \(M^{\rm th}\) state into an \((M+1)^{\rm th}\) that has been truncated away.
In the $M \rightarrow \infty$ limit, the projector can be discarded, as reflected by the expressions in Fig.~\ref{fig:one}a.

The Hamiltonian operator for \(L\) independent, well-mixed Schl\"ogl voxels is thus \(\hat{H}^{\rm rxn} = \sum_{l}^L \hat{H}^A_l + \hat H^B_l\). To add diffusion, consider the well-mixed voxels arranged in a 1D lattice as in Fig.~\ref{fig:one}b.
Each chemical species is able to move between neighboring voxels as regulated by a diffusion operator
\begin{equation}
  \hat{H}^{\rm diff} = \sum_{l = 1}^{L-1} d \left[\left(x^\dagger_{l+1} - x^\dagger_l\right)\left(x_l - x_{l+1}\right)\right],
  \label{eq:Hdiff}
\end{equation}
with hopping rate \(d = D/h^2\), consistent with macroscopic diffusion constant $D$ for voxel width $h$  \cite{smith2019spatial}.
Combining reaction with diffusion:
\begin{equation}
  \hat{H} = \hat{H}^{\rm rxn} + \hat{H}^{\rm diff},
\end{equation}
is straightforward at the level of the effective Hamiltonian, but the diffusive coupling between neighbors has the consequence of introducing nontrivial interactions between voxels.

\section{Tensor Networks}\label{sec:tn}

We have now seen how to transform the matrix-form rate operator $H$ of Sec.~\ref{sect:EstRate} into a Doi-Peliti form $\hat{H}$, 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)^{NL}$ states.
To put this scaling in perspective, the largest numerical result that follows uses $N=1$, $M=85$ and $L=8$, meaning the system contains $2.99\times 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 $\hat H$ in terms of an MPO, as shown in Fig.~\ref{fig:mps}b, and demonstrated explicitly for the Schl\"ogl mechanism in Appendix \ref{app:MPO}.
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~\cite{itensor}.

Much like the operator \(\hat{H}\) can be compressed as an MPO, the states that operators act upon can also be expressed more compactly in terms of TNs.
Consider \(\ket{p_t} = p_t(\mathbf{n}) = p_t(n_1, n_2, \hdots, n_L)\), the probability of microstate \(\mathbf{n}\) at time \(t\).
This \(\ket{p_t}\) can be thought of as a rank \(L\) tensor.
By specifying the occupancy of each of the \(L\) sites, \(\ket{p_t}\) outputs a number\textemdash a probability\textemdash associated to that microstate.
By leveraging a matrix product ansatz, $\ket{p_t}$ can be approximated as a partial contraction over $L$ different low-rank tensors (see Fig.~\ref{fig:mps}a).
 The low rank tensors of an MPS~\cite{orus2014practical} are organized in a 1D chain of tensors \(Q_{s_1}^{n_1}, Q_{s_1 s_2}^{n_2}, \hdots, Q_{s_{L-2} s_{L-1}}^{n_{L-1}},\) and \(Q_{s_{L-1}}^{n_L}\) 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 suitably chosen \(Q\) tensors, we therefore approximate
\begin{align}
\ket{p_t} &= \sum_{n_1, n_2, \hdots n_L} p_t(n_1,n_2,\ldots,n_L) \ket{n_1, n_2, \hdots n_L} \nonumber\\
&\approx \hspace{-3mm}\sum_{s_1,\ldots,s_{L-1}}\hspace{-3mm} Q^{n_1}_{s_1}Q^{n_2}_{s_1s_2}\ldots Q^{n_L}_{s_{L-1}}\ket{n_1, n_2, \hdots n_L} \nonumber\\
&\equiv \ket{q_t},
\label{eq:Pr_MPS}
\end{align}
where $\ket{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.~\ref{fig:mps}, 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.~\ref{fig:mps} emphasizes that the TN decomposition is flexible enough to address varied CME's.
The MPS of Fig.~\ref{fig:mps}a middle, has been the focus.
Each site from 1 to \(L\) reflects the number of X molecules in that voxel, but Fig.~\ref{fig:mps}a (top) and (bottom) show that the MPS decomposition can equally well be constructed when multiple species X$_i$, $i = 1,2,\ldots,N_{\rm s}$  can occupy a single voxel or when $N_{\rm s}$ different species diffuse between $L/N_{\rm s}$ voxels 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\"{o}gl model, neighboring \(Q\)'s correspond to neighboring sites in the 1D Schl\"{o}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 $s$ auxiliary indices.
If the dimension of those indices is allowed to grow exponentially from $s_1$ through $s_{L/2}$, the ansatz in Eq.~\eqref{eq:Pr_MPS} can be exact, but even when the bond dimension of the indices $s$ are capped at a maximum dimension of $\kappa$, the ansatz can be a very good approximation.
Crucially, capping that maximum bond dimension makes it practical to work with $\ket{q_t}$ because the rank-$L$ tensor $\ket{p_t}$ has been replaced by the $L$ different $Q$ tensors, each with no more than $(M+1)\kappa^2$ elements.

Having approximated both $\ket{q_t}$ and $\hat{H}$ by an MPS and MPO, respectively, we can now revisit the otherwise intractable dynamics problem of Eq.~\eqref{eq:ClSch}, only now we wish to evolve $\ket{q_t}$ in lieu of $\ket{p_t}$.
 As $\ket{q_t}$ is an MPS state, we wish for the dynamics to be constrained for all times to the manifold $\mathcal{M}$ of possible MPS states.
  Merely replacing $\ket{p_t}$ by $\ket{q_t}$ in Eq.~\eqref{eq:ClSch}:
\begin{equation}
  \frac{d \ket{q_t}}{dt} = \hat{H} \ket{q_t}
  \label{eq:qeom}
\end{equation}
would not impose such a constraint because $\hat{H}$ in general evolves $\ket{q_t}$ off $\mathcal{M}$.
An appealing resolution is given by the Time-Dependent Variational Principle (TDVP)~\cite{dirac1930note,haegeman2011time}, which proposes to first project the right-hand side onto the tangent space of the MPS manifold:
\begin{equation}
  \frac{d \ket{q_t}}{dt} = \mathcal{P}_{T_{\mathcal{M}}} \hat{H} \ket{q_t}.
  \label{eq:tdvpeom}
\end{equation}
In this way, dynamics is restricted to the desired sub-space of MPS states~\cite{vanderstraeten2019tangent}.


It turns out that the time-discretized numerical integration of this $\ket{q_t}$ is amenable to efficient TN operations.
The technical details of the time-evolution algorithm have been reported elsewhere \cite{paeckel2019time,haegeman2011time,haegeman2016unifying,lubich2015time}.
We give only a brief, high-level overview of that procedure (see also App.~\ref{app:tdvp} for more details).
The central step in deriving the TDVP time-evolution algorithm is to demonstrate that the projector $\mathcal{P}_{T_{\mathcal{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 $\delta t$, the propagator $e^{\mathcal{P}_{T_{\mathcal{M}}} \hat{H} \delta 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_ls_{l+1}}$ in an MPS.
  For an invertible $R$, replacing those tensors by $Q^{n_l}_{s_{l-1}s_l} R$ and $R^{-1} Q^{n_{l+1}}_{s_ls_{l+1}}$, respectively, does not change the MPS state; upon contracting over the neighboring sites, the identity $R R^{-1}$ would fall out.
  By leveraging the gauge freedoms, it is always possible to transform Eq.~\eqref{eq:Pr_MPS} into a so-called mixed canonical form~\cite{vidal2003efficient}:
\begin{equation}
  \ket{q_t} = A^{n_1} \hdots A^{n_{l-1}} Q^{n_l} B^{n_{l+1}} \hdots B^{n_l} \ket{n_1, \hdots, n_L},
  \label{eq:mixedcanonical}
\end{equation}
where each $A$ is a left-orthogonal tensor satisfying $A^T A = \mathbb{I}$ and each $B$ is a right-orthogonal tensor satisfying $B B^T = \mathbb{I}$, where here, $T$ represents the transpose. For notational compactness, we have suppressed implied summations over the bond indices between tensors in Eq.~\eqref{eq:mixedcanonical}.
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 $\delta t$ provided the MPS was already centered about the $Q^{n_l}$ we aim to evolve.

\begin{figure*}[hbt]
\centering
\includegraphics[width=1\textwidth]{Figures/Fig4.pdf}
\caption{\label{fig:Snaps} Snapshots of the joint density in two neighboring voxels from an $L=3$ site lattice of Schl\"{o}gl voxels.
  Using $c_1 = 2.676, c_2 = 0.040, c_3 = 108.102, c_4 = 37.881$ and $d = 8.2207$, the steady-state distribution was computed with $\kappa = 30$ by evolving an initial uniform distribution according to $\hat{H}$ via TDVP.
  That steady-state distribution was projected into $\mathcal{A}$ at time 0, renormalized, then propagated according to the modified no-recrossing dynamics $\tilde{H}$ using TDVP.
  For clarity of visualization, the density within $\mathcal{A}$ is not plotted as it would overwhelm the small probability that leaks into $\mathcal{B}$.
  By $t = 0.01$ (left image), probability has begun to leak from $\mathcal{A}$ but none has reached $\mathcal{B}$.
  Probability has just begun to reach $\mathcal{B}$ at at $t = 0.14$, roughly the $\tau_{\rm mol}$ timescale for this problem.
  By $t = 0.44$, probability is collecting in $\mathcal{B}$, and by $t = 0.75$, the flux passing from $\mathcal{A}$ to $\mathcal{B}$ has attained an approximately constant steady-state value corresponding to the linear growth phase where $P(\mathcal{B}_t|\mathcal{A}) \propto k_{\mathcal{B}\mathcal{A}}t$.}
\end{figure*}

As made more explicit in App.~\ref{app:tdvp}, the one-site TDVP algorithm initially centers the MPS around the first tensor, propagates that tensor forward in time by $\delta t$, performs a QR decomposition to change the gauge, propagates backward in time by $\delta t$ the R term of that QR decomposition 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 $\ket{q_t}$ has been propagated by a time step $\delta t$.
Whether that propagation well approximates the evolution of the full distribution $\ket{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 $\kappa$.

Notably, the time-evolved MPS state $\ket{q_t}$ allows us to estimate the probability of \emph{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.~\ref{fig:Snaps} shows the results of calculations on a 1D chain of $L = 3$ Schl\"{o}gl voxels.
Though there are three voxels, site one is traced over to leave the two-dimensional joint probability of finding $n_2^X$ X molecules in site 2 and $n_3^X$ X molecules in site 3.
At time zero, the distribution is initialized within region $\mathcal{A}$ according to $\hat{\mathcal{P}}_\mathcal{A}\ket{\boldsymbol{\pi}}/\bra{\mathbf{1}}\hat{\mathcal{P}}_{\mathcal{A}}\ket{\boldsymbol{\pi}}$.
As time progresses, a small flux of probability leaks from $\mathcal{A}$ to $\mathcal{B}$.
Because the overwhelming preference is to stay in $\mathcal{A}$, Fig.~\ref{fig:Snaps} highlights the rare transitions by only plotting the joint distribution outside of $\mathcal{A}$.
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.

\section{Many-body rate calculation}\label{sec:mbrates}
Having established the TDVP dynamics of an MPS state, we can now repeat the strategy of Eq.~\eqref{eq:rate} to compute a transition rate between states $\mathcal{A}$ and $\mathcal{B}$:
\begin{equation}
  \label{eq:tdvpdynamicsrate}
k_{\mathcal{B} \mathcal{A}} =  \frac{d}{dt} \frac{\bra{\mathbf{1}}\hat{\mathcal{P}}_{\mathcal{B}}\,e^{\mathcal{P}_{T_{\mathcal{M}}} \hat{H}t}\, \hat{\mathcal{P}}_{\mathcal{A}} \ket{\boldsymbol{\pi}}}{\bra{\mathbf{1}}\hat{\mathcal{P}}_\mathcal{A}\ket{\boldsymbol{\pi}}}\biggr\rvert_{t >\tau_{\rm mol}}.
\end{equation}
In doing so, there are a few important differences from Eq.~\eqref{eq:tdvpdynamicsrate}.
Most obviously, the TDVP dynamics includes the projection onto the MPS tangent space and $\ket{\boldsymbol{\pi}}$ must now be an MPS approximation to the steady-state distribution.
Additionally, the meaning of $\hat{\mathcal{P}}_{\mathcal{A}}$ and $\hat{\mathcal{P}}_{\mathcal{B}}$ must be adapted to now define many-body analogs of the projectors onto the two metastable states.
The single-voxel Schl\"{o}gl model had an $\mathcal{A}$ state defined as all microstates with occupancy $n$ above a threshold  $q_A^*$. Similarly, $\mathcal{B}$ included all microstates with occupancy less than some other threshold $q_B^*$.
The $L$-voxel Schl\"{o}gl model generalizes to involve hypercubic metastable states; $\mathcal{A}$ requires that each voxel satisfies $n \geq q_A^*$ and $\mathcal{B}$ requires that each voxel satisfies $n \leq q_B^*$.

With these generalizations, the transition rate of Eq.~\eqref{eq:tdvpdynamicsrate} would measure the flux crossing into $\mathcal{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 $\mathcal{A}$ then some time later enter $\mathcal{B}$ without having re-entered $\mathcal{A}$, but we want to exclude the flux that leaves $\mathcal{B}$ and re-enters $\mathcal{B}$ some time later without having returned to $\mathcal{A}$.
In Sec.~\ref{sec:rateconstants}, $k_{\mathcal{BA}}$ was sufficiently large that those re-crossings could be ignored, but as we generalize to handle smaller $k_{\mathcal{BA}}$, those recrossings can become significant.
We note that the TDVP dynamics can just as well be executed with a modified absorbing-boundary-condition dynamics with effective Hamiltonian
\begin{equation}
  \label{eq:modified}
\tilde H = \tilde H^{\rm rxn} +\tilde H^{\rm diff},
\end{equation}
where the reaction and diffusion Hamiltonians exclude events that exit $\mathcal{B}$.
An explicit construction of that modified dynamics is provided in App.~\ref{app:modifieddynamics}, allowing us to directly remove the re-crossing events.

\begin{figure*}[!ht]
\centering
\includegraphics[width=.9\textwidth]{Figures/Fig5.pdf}
\caption{\label{fig:RateL} (a) The rate $k_{\mathcal{B}\mathcal{A}}$ for switching between metastable states of the $L$-voxel Schl\"{o}gl model is calculated using the SSA (Red line) and using TDVP (dashed lines) with varying bond-dimensions $\kappa$.
  For sufficiently large bond-dimensions, TDVP agrees with a SSA rate estimate built using $400$ transitions from $\mathcal{A}$ to $\mathcal{B}$ at each $L$.
  (b) As $L$ grows a larger bond-dimension is required to capture the dynamics. For $L=3$ we see that for $\kappa=5$ and $\kappa=6$ the flow of probability from $\mathcal{A}$ to $\mathcal{B}$ is unstable.
  As $\kappa$ grows, the solutions converge.
  For this system, $\kappa=20$ is stable up to $L=6$.
  While $\kappa=30$ gives stable solutions up to $L=8$.
(c) The total computation time in hours to estimate $k_{\mathcal{B}\mathcal{A}}$ from both the SSA and TDVP.
Despite using different bond-dimensions to estimate rates on larger lattices, we see that the total time to evolve the TDVP ensemble evolution scales sub-exponentially while the SSA scales exponentially.}
\end{figure*}

\section{Numerical Evaluation of the $L$-Voxel Schl\"{o}gl Model}\label{sec:Results}
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 high-occupancy and low-occupancy states of a $L$-voxel Schl\"{o}gl model as a function of $L$.
The model is parameterized by the four stochastic rate coefficients discussed in Sec.~\ref{sect:EstRate}: $c_1 = 2.676, c_2 = 0.040, c_3 = 108.102,$ and $c_4 = 37.881$, chosen so the $L=1$ Schl\"{o}gl model will exhibit the bistablity reflected in Fig.~\ref{fig:ER}.
With those chosen rates, $\mathcal{A}$ and $\mathcal{B}$ are identified by choosing $n_a = 15$ and $n_b = 25$, and the dynamics lends itself to a finite maximum-occupancy truncation of $M = 85$.
By linking multiple voxels together in a 1D chain with a diffusive hopping rate $d = 8.2207$, 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 from a neighbor can stabilize the original state of each voxel.
Increasing $L$ thus gives a way to systematically probe rarer transition rates.

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, were generated using the SSA.
First, $5\times 10^5$ independent SSA realizations were initialized from a uniform distribution and propagated up to $t=5$ to ensure the steady state distribution was reached.
Sampled configurations that fell within $\mathcal{A}$ were then uniformly drawn as an initial configuration for 400 independent trajectories which were evolved via the SSA until the trajectory first reached $\mathcal{B}$.
From each of those $400$ trajectories, the time between the final exit of $\mathcal{A}$ and the first entry of $\mathcal{B}$, $\tau_{\mathcal{BA}}$ was recorded and the transition rate was computed from the mean of those waiting times: $k_{\mathcal{BA}} = \langle \tau_{\mathcal{BA}}\rangle^{-1}$.
Fig.~\ref{fig:RateL}a (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\"{o}gl model~\cite{vlysidis2018differences}. While the rate can be computed by trajectory sampling, the difficulty to do so (via brute force sampling) grows in proportion to the typical $\tau_{\mathcal{BA}}$, meaning the computational cost grows exponentially with $L$ (see Fig.~\ref{fig:RateL}c).
The poor scaling of brute-force SSA sampled rate calculations has motivated so-called advanced sampling methods like FFS, which focus computational effort on the progress along a reaction coordinate. 
Fig.~\ref{fig:RateL} shows that the DP/TDVP approach provides a radically different strategy for circumventing the scaling problem, reproducing the SSA rates at lower computational expense.

  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 by increasing the computational expense.
The two most prominent ways to adjust that accuracy are through the time step $\delta t$ for the TDVP evolution and through the maximum bond dimension $\kappa$ for the MPS state.
In this work, a time step of $\delta t = 10^{-4}$ was used throughout, and the trade-off between expense and accuracy was entirely tuned through the choice of $\kappa$.
The single-site TDVP algorithm we used to propagate an MPS state does not alter the bond dimension of the MPS, so $\kappa$ was entirely set by the bond dimension of the initial MPS.
We therefore needed to approximate the steady state $\ket{\boldsymbol{\pi}}$ with different values of $\kappa$.
Approximations were generated by starting with an initially uniform state using a bond dimension $\kappa$ 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 $\hat{H}$ to approach the limit $\ket{\boldsymbol{\pi}}$.
We found that evolving that $10,000$ time steps was sufficient to approximate the steady-state distribution, as reflected by the average rate of change of the distribution.
That is to say, $\bra{p_t}\ket{\dot{p}_t} = \bra{p_t}\hat{H}\ket{p_t} \approx 0$ by time $t = 1$, indicating proximity to the actual steady state, which satisfies $\ket{\boldsymbol{\dot \pi}} = \ket{\mathbf{0}}$.


To get rates, the $\kappa$-dependent MPS approximation to the steady state was projected onto $\mathcal{A}$ and renormalized.
That renormalized distribution was then propagated under single-site TDVP with the modified dynamics for another $10,000$ time steps, and projected onto $\mathcal{B}$ to add up the total probability to have reached $\mathcal{B}$.
Note that these 10,000 steps are not nearly enough to reach the steady state $\ket{\boldsymbol{\pi}}$ because the slow rate of transitioning from $\mathcal{A}$ to $\mathcal{B}$ is a bottleneck, one that did not need to be surmounted to relax from the uniform distribution to $\ket{\boldsymbol{\pi}}$ in $10,000$ steps.
As additional time steps of dynamics were taken, the growth of the transition flux from $\mathcal{A}$ to $\mathcal{B}$ was measured and plotted in Fig.~\ref{fig:RateL}b.
That figure shows that at very short times ($t < \tau_{\rm mol}$), the flux from $\mathcal{A}$ to $\mathcal{B}$ has not reached a steady state value because there hasn't been enough time for any transitions.
After that initial time $t_{\rm mol}$, the plot plateaus at a value corresponding to the rate $k_{\mathcal{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 (i.e., $\kappa = 5$ for $L = 3$), then the dynamics becomes unstable, reflected in erratic estimates for the conditional probability $P(\mathcal{B}_t|\mathcal{A}_0)$.
As the bond dimension is grown to $\kappa = 9$, the estimates converge, allowing a stable rate to be extracted.
Larger values of $L$ demand an even greater value of bond dimension.
Fig.~\ref{fig:RateL}a shows that SSA rates agree with $\kappa = 20$ calculations up to $L = 6$, but $\kappa = 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.~\ref{fig:RateL}c shows that the total time to estimate the rate grows sub-exponentially.
The blue and black dashed lines in that plot are the total computational time to estimate the rate using TDVP with $\kappa=30$ and $\kappa=20$ respectively.
Both SSA and TDVP were run using codes written in the Julia programming language \cite{bezanson2017julia} on equivalent CPUs.
While the exact wall time, of course, depends on details of the hardware and the implementation, the plot clearly demonstrates the 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$.

In addition to the favorable scaling at large $L$, we highlight that the DP/TDVP formalism allows us to extract many observables with similar ease.
For example, Fig.~\ref{fig:ExpMol} shows $\langle n_l(t)\rangle = \bra{\mathbf{1}}\hat n_l\ket{p_t}$, the average number of molecules in site $l$ with $l$ ranging from 1 through 4 for an $L = 8$ site lattice.
The state $\ket{p_t}$ was initialized in the normalized steady-state confined to $\mathcal{A}$, the high-occupation-number state.
For $t < \tau_{\rm mol}$, a rapid loss of particles near the boundary of $\mathcal{A}$ leads to a sharp drop.
After that $\tau_{\rm mol}$ timescale, the rate of decrease in $\langle n_l(t)\rangle$ becomes linear, reflecting the steady-state flux from $\mathcal{A}$ to $\mathcal{B}$.

\begin{figure}[th]
\centering
\includegraphics[width=.4\textwidth]{Figures/Fig6.pdf}
\caption{\label{fig:ExpMol}Average number of molecules for different lattice sites when $L=8$. Colors correspond to a different lattice sites. The expected number of molecules $\langle n_l(t)\rangle = \bra{\mathbf{1}}\hat n_l\ket{p_t}$ for sites $1$ to $\lfloor L/2\rfloor$ (shown) are mirrored by sites $\lfloor L/2+1\rfloor$ to $L$. The initial non-linear flow of probability out of $\mathcal{A}$ is shown by the steep drop in average number of molecules. The system passing $\tau_{mol}$ is shown by the transition to a linear change in particle number $\langle n_l(t)\rangle$.}
\end{figure}

\section{Discussion}

In this work we have shown how macroscopic rates can be efficiently calculated from the microscopic dynamics of a RDME.
We have illustrated that capability in a minimal toy model, the $L$-voxel Schl\"{o}gl model, with each voxel chemostatted such that it contains a single dynamical species (X).
In actuality, the formalism and numerical methodology extend far beyond this initial demonstration.
Any reaction-diffusion mechanism which can be written as a set of elementary reactions will have a DP representation, one which can involve raising and lowering operators for different particle types, not only for X.
For example, Table~\ref{tbl:Rxn} 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~\cite{schollwock2011density} and propagate via TDVP from that DP form of the effective Hamiltonian.
\begin{table}[ht]
\centering
\begin{tabular}{c|c}
Reaction                 & Hamiltonian                                                       \\ \hline
\ce{X ->[$\tilde{c}$] Y}        & $c[(y^\dagger - x^\dagger) x]$                           \\[1pt]
\ce{X +Y ->[$\tilde{c}$] Z}     & $c[(z^\dagger - x^\dagger y^\dagger)x y]$                \\[1pt]\ce{X + Y ->[$\tilde{c}$] 2X} & $c[x^{\dagger }(x^\dagger - y^\dagger)x y]$ \\[1pt] \hline
\end{tabular}
\caption{\label{tbl:Rxn}
  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~\cite{schulz2005exact} 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.
  }
\end{table}
As those computational tools become increasingly well developed, one can imagine that the TN tools will offer generic ``turn the crank'' analysis of microscopic models of reaction-diffusion at the ensemble level.
That vision motivates more future work at the interface of reaction-diffusion kinetics and tensor network methods.
Here, we have approximated the microscopic distribution with an MPS, allowing us to capitalize on the efficient TDVP algorithm for MPS states.
The downside is that an 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 either to avoid the need for high bond dimension, for example by leveraging TN topologies which mirror that of an underlying CRN or to adapt the bond dimension through subspace expansion \cite{hubig2015strictly}.
Those technical advances may eventually lay the groundwork to move the methodology from this Schl\"{o}gl study to the complex, biologically-motivated models where SSA realizations presently reign supreme.
 
\section{Acknowledgments}
We gratefully acknowledge Hadrien Vroylandt, Nils Strand, 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.
  The material presented in this manuscript is based upon work supported by the National Science Foundation under Grant No.\ 2141385.
  TRG is thankful that Phill Geissler could see this work in a preliminary form, but is deeply saddened that Phill's much-too-soon passing prevented him from sharing more of his wisdom.

 