Tensor Network Approach for Stochastic dynamics

I do this research in 2024 summer. PI: Prof. Hong-Yan Shih and Yi-Ping Huang, Academia Sinica

The tensor network (TN) framework is a robust approach for solving many-body quantum system problems. Techniques such as matrix product states (MPS) and the density matrix renormalization group (DMRG) provide efficient numerical tools for these challenges. Recently, physicists have found that these methods can also be applied to other systems, including stochastic dynamics problems.


Meeting Slides


Slides #1.
Slides #2.
Slides #3.
Slides #4.
Slides #5.
Slides #6.
Slides #7.

The SIS Model and the TN Approach

The one-dimensional SIS model is a simple but classic and useful model that can describe the non-equillibrium phase transitions. We can solve this model by sovliving differential equations or Monte Carlo method. As the following figures demonstrate:


The SIS model can be described in differential equation form
\(\displaystyle\frac{dS}{dt}=-\beta S I+\gamma I\)

\(\displaystyle\frac{dI}{dt}=\beta S I-\gamma I\)
There are two phases in the one-dimensional SIS model.

In differential equation form we can analyze the differential equations anr therefore have the theorem
$$\text{If \(\beta/\gamma\leq 1\),then }\lim_{t\xrightarrow{}\infty} (S,I)=(N,0)\;\text{(absorbing phase)},$$ $$\text{ If \(\beta/\gamma> 1\), then }\lim_{t\xrightarrow{}\infty} (S,I)=(\frac{\gamma}{\beta}N,(\frac{\gamma}{\beta}-1)N)\;\text{(endemic phase)}.$$ These two convergence behaviors correspond to the absorbing phase and the endemic phase respectively. The numerical simulations of these two phases are as follows:



We can obtain similar results by the Monte Carlo method as well, the result is as follows:




The Infinite-Size Density Matrix Renormalization Group (DMRG) for the Ising Model

The algorithm of the infinite-size DMRG is:


1. Start with two blocks containing a single site; these are the first system block and environment block.
2. Form two enlarged blocks by adding a site to the right of the system block and to the left of the environment block.
3. Form the superblock by connecting the system block to the environment block.
4. Diagonalize the superblock to find the state vector \(|\psi_i\rangle\) and the energy of the target state (usually ground state).
5. Form the reduced density matrix for the enlarged system block by first constructing the density matrix for the superblock \(\hat{\rho}_s=|\psi_i\rangle\langle\psi_i|\) and then taking the partial trace over the environment \(\hat{\rho}_s=Tr_e(\hat{\rho})\).
6. Diagonalize the reduced density matrix and order its eigenvectors by descending weight. Keep only the m₀ eigenvectors with the largest weights.
7. Form the truncation operator \(\hat{O}\) from the eigenvectors with larger eigenvalues and project all operators onto the truncated basis \(\hat{A}_{L+1}=\hat{O}^\dag\hat{A}_L\hat{O}\).
8. Reflect the system block to get the new environment block.
9. Repeat steps 2 - 8 until certain conditions are satisfied.
As an isslustration, see the following figure:



Construct two blocks (system and eviroment) with some computationable finite size, add twos sites to form a super block. Truncate the basis of the operators to form a new blocks with the same dimension.

Add two sites again. Repeat until certain conditions are satisfied.
The Hamiltonian operator of the Ising model is $$\hat{H}=-J\sum_i \sigma^x_i\otimes\sigma^x_{i+1}-h\sum_i \sigma_i^z,$$ where \(\sigma_i^\alpha\) are Pauli matrix.
Consider the system part with an additional site first. Take \(J=h=1\) and write the Hamiltonian for a single in explicit matrix form:
the first part (interaction part) is $$\hat{H}_{interaction}=-\sigma^x_i\otimes \sigma^x_{i+1} =-\begin{pmatrix}0&1\\ 1&0\end{pmatrix}\otimes \begin{pmatrix}0&1\\ 1&0\end{pmatrix}=\begin{pmatrix}0&0&0&-1\\ 0&0&-1&0\\ 0&-1&0&0\\ -1&0&0&0\\ \end{pmatrix}$$ and the second part (external part) is $$\hat{H}_{external}=-(\sigma^x_i\otimes \mathbb{I}+ \mathbb{I}\otimes\sigma^x_i) =\begin{pmatrix}-1&0\\ 0&1\end{pmatrix}\otimes \begin{pmatrix}1&0\\ 0&1\end{pmatrix}+\begin{pmatrix}1&0\\ 0&1\end{pmatrix}\otimes \begin{pmatrix}-1&0\\ 0&1\end{pmatrix}=\begin{pmatrix}-2&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&2\\ \end{pmatrix}.$$
This is the Hamiltonian for both of the system \(\hat{H}_s\) and the eviroment \(\hat{H}_e\).
$$\hat{H}_s=\hat{H}_e=\begin{pmatrix}-2&0&0&-1\\ 0&0&-1&0\\ 0&-1&0&0\\ -1&0&0&2\\ \end{pmatrix}$$
Since the additional sites of the system and the enviroment are $$\sigma^x_s=\mathbb{I}\otimes \sigma^x=\begin{pmatrix}0&1&0&0\\ 1&0&0&0\\ 0&0&0&1\\ 0&0&1&0\\ \end{pmatrix},$$ $$\sigma^x_e=\sigma^x\otimes\mathbb{I}=\begin{pmatrix}0&0&1&0\\ 0&0&0&1\\ 1&0&0&0\\ 0&1&0&0\\ \end{pmatrix}.$$ the term describing the interaction between the system and the enviroment \(\hat{H}_{s,e}\) is $$\hat{H}_{s,e}=\sigma^x_s\otimes\sigma^x_e=\begin{pmatrix}0&1&0&0\\ 1&0&0&0\\ 0&0&0&1\\ 0&0&1&0\\ \end{pmatrix}\otimes \begin{pmatrix}0&0&1&0\\ 0&0&0&1\\ 1&0&0&0\\ 0&1&0&0\\ \end{pmatrix}$$
The Hamiltonian of the super block is
\[ \begin{align*} \hat{H}_{super}&=\hat{H}_{s,e}+\hat{H}_e\otimes \mathbb{I}+\mathbb{I}\otimes\hat{H}_s \\ &\\ &=\begin{pmatrix}0&1&0&0\\ 1&0&0&0\\ 0&0&0&1\\ 0&0&1&0\\ \end{pmatrix}\otimes \begin{pmatrix}0&0&1&0\\ 0&0&0&1\\ 1&0&0&0\\ 0&1&0&0\\ \end{pmatrix}+\begin{pmatrix}0&0&0&-1\\ 0&0&-1&0\\ 0&-1&0&0\\ -1&0&0&0\\ \end{pmatrix}\otimes \mathbb{I}+\mathbb{I}\otimes\begin{pmatrix}0&0&0&-1\\ 0&0&-1&0\\ 0&-1&0&0\\ -1&0&0&0\\ \end{pmatrix}\\ &\\ &= \begin{pmatrix} -4 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & -2 & -1 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & -1 & -2 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\ -1 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & -1 & 0 & -2 & 0 & 0 & -1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & -1 & 0 & 0 & 2 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & -2 & 0 & 0 & -1 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & -1 & 0 & 0 & 2 & 0 & -1 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & -1 \\ 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 2 & -1 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & -1 & 2 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & -1 & 0 & 0 & 4 \end{pmatrix} \end{align*} \]
Obtain the ground state and its energy by diagonalize the Hamiltonian. $$|\psi_0\rangle=\begin{pmatrix} 0.902 \\ -0.000 \\ -0.000 \\ 0.225 \\ -0.000 \\ 0.109 \\ 0.235 \\ -0.000 \\ -0.000 \\ 0.058 \\ 0.109 \\ 0.000 \\ 0.225 \\ 0.000 \\ 0.000 \\ 0.058 \end{pmatrix},\;E_0=-4.7588.$$ Construct the density matrix of the ground state \(\hat{\rho}=|\psi_0\rangle\langle\psi_0|\) and the partial trace with respect to the enviroment \(\hat{\rho}_s=Tr_e(\hat{\rho})\): $$\begin{pmatrix} 0.864 & -0.000 & -0.000 & 0.216 \\ -0.000 & 0.067 & 0.032 & -0.000 \\ -0.000 & 0.032 & 0.015 & -0.000 \\ 0.216 & -0.000 & -0.000 & 0.054 \end{pmatrix}.$$ Obtain the eigenvectors of \(\hat{\rho}_s\)
\[ \begin{pmatrix} -0.242 \\ 0.000 \\ -0.000 \\ 0.970 \end{pmatrix}, \quad \begin{pmatrix} 0.000 \\ 0.429 \\ -0.903 \\ 0.000 \end{pmatrix}, \quad \begin{pmatrix} 0.000 \\ -0.903 \\ -0.429 \\ 0.000 \end{pmatrix}, \quad \begin{pmatrix} 0.970 \\ -0.000 \\ -0.000 \\ -0.242 \end{pmatrix} \] with eigenvalues \[ \lambda_1 = 0.000, \quad \lambda_2 = 0.000, \quad \lambda_3 = 0.082, \quad \lambda_4 = 0.918 \]
Select two eigenvectors with larger eigenvalues to obtain the truncation operator: \[ \hat{O} = \begin{pmatrix} 0.970 & -0.000 \\ -0.000 & -0.903 \\ -0.000 & -0.429 \\ 0.242 & -0.000 \end{pmatrix} \] Use the tuncation operator to obtain the new blocks.
The first ones are the new blocks, including both system and eviroment: \[ \sigma^z_{\text{new}} = \hat{O}^\dag \hat{H}_s \hat{O} = \hat{O}^\dag \hat{H}_e \hat{O} = \begin{pmatrix} -2.24 & -0.000 \\ -0.000 & -0.775 \end{pmatrix} \] The second one is the new interaction term of system and eviroment, which can be described by \(\sigma^x\): \[ \sigma^x_{s,\text{new}} = \hat{O}^\dag \sigma^x_s \hat{O} = \begin{pmatrix} -0.000 & -0.980 \\ -0.980 & 0.000 \end{pmatrix} \]
We can therefore construct the Hamiltonian of the new super block by these new \(\sigma^\alpha_{new}\). We can iteratively do this process until certain condition reached.

The numerical result of the groud state energy is as follow:


We know that the ground state of the infinite Ising model is negative infinity, which is the same as the numericaal result.

The Finite-Size DMRG for the Ising Model

The algorithm of the finite-size DMRG is:

1. Carry out the steps of the infinite system algorithm until the superblock reaches a predetermined size \( L \). The size of the system block should now be \( l = L/2 \). Blocks should be stored for all sizes, from \( l = 1 \) to \( l = L/2 \).
2. [Start of sweep] Grow the system block to size \( l + 1 \). Retrieve a block of size \( L - l - 2 \) from storage and reflect it to get an environment block. Grow the environment block to size \( L - l - 1 \).
3. Carry out steps 3 - 7 of the infinite system algorithm.
4. Store the truncated block as the new block.
5. Set \( l = l + 1 \) and repeat steps 2 - 4 until \( l = L/2 - 2 \). At this point, set \( l = 1 \) and retrieve a block of size 1 to be the new system block.
6. [End of sweep] Repeat steps 2 - 5 until \( l = L/2 \). This concludes a sweep.
7. Repeatedly sweep over the system (steps 2 - 6) until certain conditions are satisfied.
As an isslustration, see the following figure:



Construct a finite super block, add twos sites and sweep them form the left two the right. Repeatly sweep until certain conditions are satisfied.

DMRG in Matrix Product State (MPS) and Matrix Product Operator (MPO) Language

The great success behind DMRG is MPS and MPO. The DMRG algorithm is actually the variational in the MPS space. Since the state of the quantum system is $$|\Psi\rangle=\sum_{s_1,\cdots,s_N}\Psi_{s_1s_2\cdots s_N}\bigotimes_{i=1}^N|s_i\rangle,$$ where \(\Psi_{s_1\cdots s_N}\) is a rank-\(N\) tensor and can therefore be decomposed by singular value decomposition (SVD). The quantum state is therefore $$|\Psi\rangle=\sum_{s_1,\cdots,s_N}A_1^{s_1}A_2^{s_2}\cdots A_3^{s_N}\bigotimes_{i=1}^N|s_i\rangle$$ where \(A_i^{s_i}\) are matrices. The quantum state can be described by matrix product, this state is called MPS. By the same argument, we can decompose the quantum operators by SVD as well. The quantum operators are written as $$\hat{O}=Q_1^{s_1}Q_2^{s_2}\cdots Q_3^{s_N},$$ where \(Q_i^{s_i}\) are matrices.

The Finite-Size variational Method and the DMRG

Finding the ground state in the finite-size system is $$E_0=\inf_{|\Psi\rangle}\frac{\langle\Psi|\hat{H}|\Psi\rangle}{\langle\Psi|\Psi\rangle}.$$ We can definitely obtain the ground state energy and its eigenstate (ground state) by solving the optimazation problem $$\min_{|\Psi\rangle}\frac{\langle\Psi|\hat{H}|\Psi\rangle}{\langle\Psi|\Psi\rangle}.$$ The most intuitive way to solve this is diagonalizing the Hamiltonian. However, this is impossible for huge systems. Since we know that the quantum states and the Hamiltonian can be written as MPS and MPO. We can perform the optimization tensor by tensor. Consider the effective Hamiltonian \(\hat{H}_{eff}\) with respect to \(A\), \(\langle\Psi|\hat{H}|\Psi\rangle\) without \(A^\dag\) and \(A\). The local optimization problem is $$\min_{A}(A^\dag\hat{H}_{eff}A-\lambda (A^\dag \hat{N} A-1)),$$ where \(\lambda\) is the Lagrangian multiplier to constraint that \(\langle\Psi|\Psi\rangle=1\). It is natural to consider the gradient with respect to the tensor \(A^\dag\) to perform the optimization, that is, $$\nabla_{A^\dag}(A^\dag\hat{H}_{eff}A-\lambda (A^\dag \hat{N} A-1))=0,$$ which leads to a generalized eigenvalue problem naturally $$\hat{H}_{eff}A=\lambda \hat{N} A.$$ Now we discuss why the finite-size DMRG is equivalent to performing the variational method in MPS space. This notion is similar to the conjugate gradient method, gradient descendent method, and the tangent space method. We consider the optimization in a subspace of the complete parameter space to make sure that the computational complexity will not be to large. We can approximately obtain the solution by iterations.

The finite-size DMRG is almost equivalent to the variational method explained above, the only difference is that the super blocks in the finite-size DMRG having two sites and the super block in the variational method having only one site. The algorithm of the finite-size DMRG is:

1. Generate MPSs and MPOs: Create the Matrix Product States (MPSs) and Matrix Product Operators (MPOs) for the system.
2. [Start of Right Sweep] Optimization: Optimize the matrices by diagonalizing the chosen matrices within the indexed tensor of the effective Hamiltonian.
3. Normalization: Perform Singular Value Decomposition (SVD) to normalize the matrices on the left-hand side, integrating the left matrices into the subsequent tensor.
4. [End of Right Sweep] Iteration: Use the newly obtained matrices as the starting point for the next tensor. Repeat steps 2 and 3 until the right boundary of the system is reached.
5. [Start of Left Sweep] Optimization: Optimize the matrices by diagonalizing the chosen matrices within the indexed tensor of the effective Hamiltonian.
6. Normalization: Perform SVD to normalize the matrices on the left-hand side, integrating the left matrices into the subsequent tensor.
7. [End of Left Sweep] Iteration: Use the newly obtained matrices as the starting point for the next tensor. Repeat steps 5 and 6 until the left boundary of the system is reached.
8. Iterative Sweeping: Repeat the entire sweeping process (steps 2 through 7) until convergence criteria are met.




Suppose that we have 50 sites in the system. The MPO has 50 sites also, it can be expressed as $$\hat{H}=\hat{L}\hat{M}\cdots\hat{M}\hat{R},\;\hat{M}=\begin{pmatrix} \mathbb{I}&\sigma^z&-h\sigma^x\\ 0&0&-J\sigma^z\\ 0&0&\mathbb{I} \end{pmatrix},\;\hat{L}=\begin{pmatrix} 1\\ 0\\ 0\end{pmatrix}^T\hat{M},\;\hat{R}=\hat{M}\begin{pmatrix}0\\ 0\\ 1\end{pmatrix}.$$

Take \(\chi=16\), the MPS will have the size as follows:

Size Information

Index Size Information 1 1x2x2 double 2 2x2x4 double 3 4x2x8 double 4 8x2x16 double 5 16x2x16 double 6 16x2x16 double 7 16x2x16 double 8 16x2x16 double 9 16x2x16 double 10 16x2x16 double

Index Size Information 11 16x2x16 double 12 16x2x16 double 13 16x2x16 double 14 16x2x16 double 15 16x2x16 double 16 16x2x16 double 17 16x2x16 double 18 16x2x16 double 19 16x2x16 double 20 16x2x16 double

Index Size Information 21 16x2x16 double 22 16x2x16 double 23 16x2x16 double 24 16x2x16 double 25 16x2x16 double 26 16x2x16 double 27 16x2x16 double 28 16x2x16 double 29 16x2x16 double 30 16x2x16 double

Index Size Information 31 16x2x16 double 32 16x2x16 double 33 16x2x16 double 34 16x2x16 double 35 16x2x16 double 36 16x2x16 double 37 16x2x16 double 38 16x2x16 double 39 16x2x16 double 40 16x2x16 double

Index Size Information 41 16x2x16 double 42 16x2x16 double 43 16x2x16 double 44 16x2x16 double 45 16x2x16 double 46 16x2x16 double 47 16x2x8 double 48 8x2x4 double 49 4x2x2 double 50 2x2x1 double



\(\hat{L}\) and \(\hat{R}\) are rank-3 tensors with \(\hat{M}\) is rank-4 tensor. The matrices in the MPS are rank-3 tensors except those on the boundary are rank-2 tensors.
The diagram of \(\langle\Psi|\hat{H}|\Psi\rangle\) can be expressed as:



With 4 sweeps (4\(\times\)2\(\times\)50 iterations), we can have the result as follows:



The value of \(\chi\) does not change the convergence rate of the algorithm but is highly related to the upper bound of its effectiveness. Moreover, the smaller the number of sites is, the more exact the result is. Suppose that we have a twenty-site system, the reslut is as follows:



We can see that the error can get to about \(10^{-8}\) with 4 sweeps.

Notice that we take \(J=H=1\), which is the critical point, here. The divergence of the correlated length makes we need larger bond dimension to get the accurate result. If we chage the parameters to \(J=50\) and \(h=1\), which is far form the critical point, then the result is as follows:



We can see the the results obtained from larger and smaller bond dimension are not very differenta and the resluts are more accurate.

The TN Approach fot the SIS Model



The Hamiltonian of the SIS model is $$\hat{W}=\lambda\sum_{i=1}^{N-1}(\hat{n}_i\hat{w}_{i+1}^{0\xrightarrow{}1} +\hat{w}_i^{0\xrightarrow{}1}\hat{n}_{i+1})+\sum_{i=1}^{N-1}\hat{w}_i^{1\xrightarrow{}0} +\hat{W}_{driv}(\alpha),$$ where \(\hat{n}=|1\rangle\langle1|\), \(\hat{w}^{0\xrightarrow{}1}=|1\rangle\langle 0|-|0\rangle\langle 0|\) , \(\hat{w}^{1\xrightarrow{}0}=|0\rangle\langle 1|-|1\rangle\langle 1|.\)
The MPO is $$\hat{H}=\hat{L}\hat{M}\cdots\hat{M}\hat{R},\;\hat{M}=\begin{pmatrix} \mathbb{I}&0&0&0\\ \lambda\hat{w}^{0\xrightarrow{}1}&0&0&0\\ \lambda\hat{n}&0&0&0\\ \hat{w}^{1\xrightarrow{}0}&\hat{n}&\hat{w}^{0\xrightarrow{}1}&\mathbb{I} \end{pmatrix},\;\hat{L}=\begin{pmatrix} 0\\ 0\\ 0\\ 1\end{pmatrix}^T\hat{M},\;\hat{R}=\hat{M}\begin{pmatrix}1\\ 0\\ 0\\ 0\end{pmatrix}.$$







As an accurate enough \(\lambda=0\) is enough, we also get the eigenstate, that is, the NESS. However, we can not employ it directly since it is in MPS form. To utilize it, we need to construct some operators to project them to certain quantities.
If we want to compute the expected value of number of gaps with length \(k\), we can construct an operator $$\hat{L}_k\hat{M}_k\cdots\hat{M}_k\hat{R}_k,$$ where \(\hat{M}_k=\begin{pmatrix}\mathbb{I}&\hat{n}&0&\cdots&&0\\ 0&0&\hat{v}&&&0\\ &&&\ddots&&\vdots\\ \vdots&&&&\hat{v}&0\\ &&&&&\hat{n}\\ 0&&&\cdots&0&\mathbb{I} \end{pmatrix},\;\hat{L}_k=\begin{pmatrix}1&0&\cdots&0\end{pmatrix}\hat{M}_k,\;\hat{R}_k=\hat{M}_k\begin{pmatrix}0\\ \vdots\\0\\1\end{pmatrix}\).

Application in Ecology

Consider a one-dimensional lattice with \(L\) sites, where each site can be in one of four states: empty, occupied by a prey, occupied by a predator, or occupied by both a prey and a predator. The dynamics of the model are governed by the following algorithm:
1. Initialize the lattice with a uniform distribution of preys and predators, ensuring that no site is initially occupied by both a prey and a predator.
2. Prey Reproduction: Implement this step with probability \(r\). At each time step, for every site occupied by a prey, the prey can reproduce by adding an offspring to a neighboring site that is either empty or occupied by a predator (i.e., a site without a prey).
3. Predator Predation and Reproduction: Implement these step with probability \(1-r\). At each time step, if a site is occupied by both a prey and a predator, the predator consumes the prey. Subsequently, the predator may reproduce by adding an offspring to a neighboring site that is either empty or occupied by a prey (i.e., a site without a predator).
4. Predator Death: Implement these step with probability \(1-r\). At each time step, if a site is occupied by a predator only (with no prey present), the predator dies due to starvation.
The parameter of this model is \(r\), the critical point of this model is \(r\approx 0.177\). When \(r < 0.177 \), the system is in the absorbing phase and the predators will ultimately extinct. When \(r > 0.177 \), the system is in the active phase and the predators will keep competing.




The TN Approach

$$\begin{align*} \hat{W}=&\sum_{i=1}^{N-1}\bigg((\hat{n}_i^1\hat{w}^{01}_{i+1}+\hat{w}^{01}_i\hat{n}^1_{i+1}) +(\hat{n}_i^2\hat{w}^{13}_{i+1} +\hat{w}^{13}_i\hat{n}^2_{i+1})+(\hat{n}_i^1\hat{w}^{23}_{i+1} +\hat{w}^{23}_i\hat{n}^1_{i+1}) +(\hat{n}_i^1\hat{w}^{01}_{i+1}+\hat{w}^{01}_i\hat{n}^1_{i+1})\\ &+(\hat{n}_i^3\hat{w}^{01}_{i+1}+\hat{w}^{01}_i\hat{n}^3_{i+1}) +(\hat{n}_i^3\hat{w}^{13}_{i+1} +\hat{w}^{13}_i\hat{n}^3_{i+1})+(\hat{n}_i^3\hat{w}^{23}_{i+1} +\hat{w}^{23}_i\hat{n}^3_{i+1}) +(\hat{n}_i^3\hat{w}^{01}_{i+1}+\hat{w}^{01}_i\hat{n}^3_{i+1})\bigg)\\ &+\sum_{i=2}^{N-1}\bigg((\hat{n}_{i-1}^1\hat{w}^{03}_i\hat{n}_{i+1}^2 +\hat{n}^2_{i-1}\hat{w}^{03}_i\hat{n}^1_{i+1})+(\hat{n}_{i-1}^3\hat{w}^{03}_i\hat{n}_{i+1}^2 +\hat{n}^3_{i-1}\hat{w}^{03}_i\hat{n}^1_{i+1}) +(\hat{n}_{i-1}^1\hat{w}^{03}_i\hat{n}_{i+1}^3 +\hat{n}^2_{i-1}\hat{w}^{03}_i\hat{n}^3_{i+1})\\ &+(\hat{n}_{i-1}^3\hat{w}^{03}_i\hat{n}_{i+1}^3 +\hat{n}^3_{i-1}\hat{w}^{03}_i\hat{n}^3_{i+1})\bigg)+(1-r)\sum_{i=1}^{N-1}\bigg(\hat{n}_i^3\hat{w}^{03}_{i+1}+\hat{w}^{03}_i\hat{n}^3_{i+1}\bigg) +\sum_{i=1}^{N}\bigg(\hat{w}^{20}_{i}+\hat{w}^{32}_i\bigg) \end{align*}$$

References


[1] Efficient simulations of epidemic models with tensor networks: Application to the one-dimensional susceptible-infected-susceptible model.
[2] Density-matrix renormalization group: a pedagogical introduction.
[3] Non-equilibrium phase transitions.