Numerical Simulation of Critical Dissipative Non-Equilibrium Quantum Systems with an Absorbing State

TL;DR

This paper explores the use of tensor network methods, specifically matrix product states and time-evolving block decimation, to simulate the critical dynamics of a non-equilibrium quantum system with an absorbing state, revealing insights into its universality class.

Contribution

It demonstrates the effectiveness of the Heisenberg picture and quantum trajectories in improving simulation accuracy for the quantum contact process at criticality.

Findings

Heisenberg picture improves simulation accuracy over Schrödinger approach

Quantum trajectories can simulate longer timescales

Quantum contact process does not belong to directed percolation universality class

Abstract

The simulation of out-of-equilibrium dissipative quantum many body systems is a problem of fundamental interest to a number of fields in physics, ranging from condensed matter to cosmology. For unitary systems, tensor network methods have proved successful and extending these to open systems is a natural avenue for study. In particular, an important question concerns the possibility of approximating the critical dynamics of non-equilibrium systems with tensor networks. Here, we investigate this by performing numerical simulations of a paradigmatic quantum non-equilibrium system with an absorbing state: the quantum contact process. We consider the application of matrix product states and the time-evolving block decimation algorithm to simulate the time-evolution of the quantum contact process at criticality. In the Lindblad formalism, we find that the Heisenberg picture can be used to improve the accuracy of simulations over the Schr\"{o}dinger approach, which can be understood by considering the evolution of operator-space entanglement. Furthermore, we also consider a quantum trajectories approach, which we find can reproduce the expected universal behaviour of key observables for a significantly longer time than direct simulation of the average state. These improved results provide further evidence that the universality class of the quantum contact process is not directed percolation, which is the class of the classical contact process.

Tables

Table 1. Table 1: Relevant exponents of the QCP and CCP. The exponents δ , z 𝛿 𝑧 \delta,z and Θ Θ \Theta can be associated to non-equilibrium observables when starting from an initial seed-state, see ( 9 ) - ( 11 ). The exponent α 𝛼 \alpha can be associated to the decay of excitation density when starting from a homogenous fully-occupied state, Henkel  et al. ( 2008 ) . The exponents of the CCP are given by the DP universality class Hinrichsen ( 2000 ) . The QCP estimates of δ , z 𝛿 𝑧 \delta,z and Θ Θ \Theta are given in this work (see Section VI ) with two standard errors, while the value of α 𝛼 \alpha indicated by an asterisk is estimated in Carollo  et al. ( 2019 ) . The two final columns give exponents for the QCP with the inclusion of classical branching at the tricritical point where quantum and classical branching compete, estimated using functional renormalisation group Buchhold  et al. ( 2017 ) . Note that for DP α = δ 𝛼 𝛿 \alpha=\delta , a consequence of the rapidity reversal symmetry that is characteristic of the class Henkel  et al. ( 2008 ) .

	QCP	$1d$ DP	$2d$ DP	$1d$ Ref. Buchhold et al. (2017)	$2d$ Ref. Buchhold et al. (2017)
$\delta$	$0.26\pm 0.04$	0.16	0.45	$-$	$-$
$z$	$1.61\pm 0.16$	1.58	1.77	1.93	1.97
$\Theta$	$0.26\pm 0.05$	0.31	0.23	$-$	$-$
$\alpha$	$0.36$ $\pm$ ${0.08}^{\ast }$	0.16	0.45	0.21	0.35

Equations

$$\frac{d}{dt}\rho (t)=\mathcal{L}[\rho (t)],$$

$$\mathcal{L}[X]=-i[H,X]+\mathcal{D}[X].$$

$$\mathcal{D}[X]=\sum _{\mu }\left(J_{\mu }X{J}_{\mu }^{†}-\frac{1}{2}\left\{X,J_{\mu }^{†}{J}_{\mu }\right\}\right),$$

$${\mathcal{D}}_{\text{cl}}[\rho ]=\mathcal{D}[\rho ]+{\mathcal{D}}_{\mathrm{br}}[\rho ].$$

$$\mathcal{D}[\rho ]=\gamma \sum _{k=1}^L\left({\sigma }_{-}^{(k)}\rho {\sigma }_{+}^{(k)}-\frac{1}{2}\left\{n^{(k)},\rho \right\}\right),$$

$$\begin{array}{rl}{\mathcal{D}}_{\mathrm{br}}[X]& =\Gamma \sum _{k=1}^{L-1}{J}_{\mathrm{L}}^{(k)}X(J_{\mathrm{L}}^{(k)}{)}^{†}-\frac{1}{2}\{(J_{\mathrm{L}}^{(k)}{)}^{†}{J}_{\mathrm{L}}^{(k)},X\}+\\ & +\Gamma \sum _{k=1}^{L-1}{J}_{\mathrm{R}}^{(k)}X(J_{\mathrm{R}}^{(k)}{)}^{†}-\frac{1}{2}\{(J_{\mathrm{R}}^{(k)}{)}^{†}{J}_{\mathrm{R}}^{(k)},X\},\end{array}$$

$$H=\Omega \sum _{k=1}^{L-1}\left({\sigma }_1^{(k)}{n}^{(k+1)}+n^{(k)}{\sigma }_1^{(k+1)}\right).$$

$$P_{\text{sur}}(t)$$

$$N_{\text{a}}(t)$$

$$n_{\text{seed}}(t)$$

$$P_{\text{sur}}(t)$$

$$N_{\text{a}}(t)$$

$$n_{\text{seed}}(t)$$

$$|\psi ⟩=\sum _{s_1,s_2,...,s_L=1}^K{\mathbf{M}}_{s_1}^{[1]}{\mathbf{M}}_{s_2}^{[2]}...{\mathbf{M}}_{s_L}^{[L]}|s_1{s}_2...s_L⟩,$$

$$S({\varrho }_A)=-\text{tr}\left[{\varrho }_A\log {\varrho }_A\right]=-\text{tr}\left[{\varrho }_B\log {\varrho }_B\right],$$

$${\varrho }_{A/B}={\text{tr}}_{B/A}|\psi ⟩⟨\psi |,$$

$$|{\varphi }_{t^{\prime }}⟩={\text{argmin}}_{|\varphi ⟩}||{\psi }_{t^{\prime }}⟩-|\varphi ⟩|,$$

$${ϵ}_{\chi }^{[k]}=\sum _{j\ge \chi +1}{\left({\lambda }_j^{[k]}\right)}^2,$$

$$\rho (t)=\sum _{mn}{\rho }_{mn}(t)|m⟩⟨n|\to |\rho (t)⟩=\sum _{mn}{\rho }_{mn}(t)|m⟩\otimes |n⟩.$$

$$\frac{d}{dt}|\rho (t)⟩=\mathbb{L}|\rho (t)⟩,$$

$$\mathbb{L}=-i{H}_D+{\mathcal{D}}_D,$$

$$H_D=\left(H\otimes \text{mathds}1-\text{mathds}1\otimes H^T\right),$$

$${\mathcal{D}}_D=\sum _{\mu }\left(J_{\mu }\otimes J_{\mu }^{\ast }-\frac{1}{2}{J}_{\mu }^{†}{J}_{\mu }\otimes \text{mathds}1-\frac{1}{2}\text{mathds}1\otimes J_{\mu }^T{J}_{\mu }^{\ast }\right),$$

$$O(t)$$

$$\tilde{S}=-\text{tr}\left[{\tilde{\varrho }}_A\log {\tilde{\varrho }}_A\right]=-\text{tr}\left[{\tilde{\varrho }}_A\log {\tilde{\varrho }}_A\right],$$

$${\tilde{\varrho }}_{A/B}={\text{tr}}_{B/A}|\rho ⟩⟨\rho |.$$

$$|O(t)⟩=e^{t{\mathbb{L}}^{†}}|O⟩.$$

Full text

I Introduction

Despite the recent experimental progress in probing the emergent behaviour of out-of-equilibrium ensembles of cold atoms or trapped ions, Syassen et al. (2008); Kim et al. (2010); Barreiro et al. (2011); Bohnet et al. (2016); Lienhard et al. (2018); Wade et al. (2018), a clear understanding of these quantum non-equilibrium systems remains a major challenge. While in classical settings the study of such systems – of their phases and critical phenomena – are well developed, options for going beyond semi-classical treatments or the physics of exactly solvable quantum models are rather limited. This is especially true for open (dissipative) systems, which are of interest as they have potential to display a rich set of novel non-equilibrium physics – e.g. critical dynamical behaviour or dynamical phase transitions – not possible in closed (unitary) settings.

Here, we explore the simulation of critical dynamics in dissipative quantum many body systems, in the case of the quantum contact process (QCP) Marcuzzi et al. (2016); Buchhold et al. (2017); Roscher et al. (2018); Jo et al. (2019); Carollo et al. (2019). The QCP is an attractive model to study for a number of reasons: Firstly, the QCP is the coherent version of the well-understood classical contact process (CCP) Hinrichsen (2000), which exhibits a non-equilibrium phase transition (NEPT) in the directed percolation (DP) universality class, even in $1d$. Formulating both the CCP and QCP in the Lindblad formalism then allows for a direct comparison between the performance of a given numerical approach in the classical and quantum cases.

Secondly, the QCP contains an absorbing state, which has been suggested to make simulations of quantum systems more challenging Carollo et al. (2019). This idea can be explored by comparing the simulation of dynamics performed in the Schödinger picture with that of the Heisenberg picture: In the contact process, as in other dissipative systems, there is an asymmetry between dynamics in the Schrödinger picture and Heisenberg picture. In particular, in the Heisenberg picture the absorbing state is absent. It is then of interest to investigate any difference in performance between simulations in the two. Finally, the critical physics of the $1d$ QCP is similar to that of the CCP, in the sense that key observables display power-law behaviour but with different exponents. This means that, at criticality, different numerical methods or approaches to the dynamics can be compared by their ability to reproduce the expected power-laws. Furthermore, since the universality class of the QCP is currently debated, Marcuzzi et al. (2016); Roscher et al. (2018); Carollo et al. (2019), it is of considerable interest in its own right to make estimates of critical exponents, comparing these with previous estimates and known cases.

To simulate the non-equilibrium dynamics of the QCP, we apply matrix product states (MPSs) and the time-evolving block-decimation (TEBD) algorithm Vidal (2004); Schollwöck (2011); Paeckel et al. (2019). This algorithm belongs to a more general class of tensor network (TN) methods, well established for the simulation of closed quantum systems in $1d$, which have also been applied to dissipative quantum systems previously in a number of cases Bonnes and Läuchli (2014); Jaschke et al. (2019); Cui et al. (2015); Mascarenhas et al. (2015); Gangat et al. (2017); Kshetrimayum et al. (2017); De las Cuevas et al. (2013, 2016); Werner et al. (2016). In the context of studying dissipative quantum dynamics, a key question for TN methods is whether different approaches, such as quantum trajectories (QTs) as opposed to the Lindblad master equation, can lead to substantially different accuracies. This question has been explored previously in Bonnes and Läuchli (2014); Jaschke et al. (2019) and it has been suggested that in high-entanglement scenarios QTs might prove more accurate.

In the case of simulating the critical QCP with TEBD, we find the following key results: First, in the Lindblad formalism, the Heisenberg picture can be used to improve the accuracy of simulations beyond that of the Schrödinger picture and this can be explained by considering the evolution of operator-space entanglement entropy. Second, we find that a QTs approach leads to a significant improvement in the approximation of key universal observables, i.e., the reproduction of power-law behaviour for longer times. Finally, using the results from QTs, we show that the estimated exponents of the QCP lie far from those of DP, providing further evidence that the QCP belongs to a universality class different to DP.

While we focus on the application of TEBD to the critical dynamics of the QCP, we note there are also a number of other methods available with which it would be interesting to compare results. For example, cluster mean-field Jin et al. (2016), variational minimisation Weimer (2015) and algorithms based on time-dependent variational Monte Carlo with neural networks have all shown promise for the simulation of open quantum systems Carleo and Troyer (2017); Nagy and Savona (2019); Vicentini et al. (2019); Hartmann and Carleo (2019); Yoshioka and Hamazaki (2019).

The layout of the paper is as follows: In Section II we discuss the classical and quantum contact processes. In Section III we provide a brief overview of the TEBD algorithm and MPSs for the unfamiliar reader. Section IV then shows the results for the simulation of the QCP in the Lindblad formalism, comparing with the CCP. Section V examines the Heisenberg picture for the QCP, while Section VI examines a QTs approach. Conclusions and outlook are contained in Section VII.

II The Classical and Quantum Contact Processes

II.1 Lindblad Formalism for Open Quantum Dynamics

In many physically relevant scenarios, it is not possible to characterise the time-evolution of quantum system by means of unitary closed dynamics. In fact, the simple presence of thermal surroundings or of stochastic processes, e.g. measurements of the system or of its output, makes the study of the time-evolution of the system much more involved. These open quantum systems (OQSs) are not described by usual Schrödinger equation. Instead, their average dynamics is implemented, in the Markovian regime, by means of Lindblad-type master equations. These describe the evolution of the average quantum system state $\rho(t)$, where average might mean expectation over a stochastic process or over the degrees of freedom of an environment, according to the physical scenario one has in mind.

Concretely, the quantum state of the system $\rho(t)$ obeys the Lindblad differential equation Lindblad (1976),

[TABLE]

where the map $\mathcal{L}$, also known as the Lindblad generator, is made of two different contributions,

[TABLE]

The first piece represents the usual commutator with the system Hamiltonian, which implements the coherent part of the dynamics. On the other hand, the term $\mathcal{D}$, often called the dissipator, contains information about the effects of the environment, or of stochasticity in general, on the system. This must have the following form,

[TABLE]

with $\{A,B\}=AB+BA$, in order to guarantee a physical meaning for the evolved quantum state $\rho(t)$. Note that the dynamics implemented by (1) does not preserve the purity of states and, therefore, in general $\rho(t)$ will represent a mixed state.

II.2 Model Definitions

As illustrated in Fig. 1, in a contact process, lattice sites can be either occupied ($\ket{\bullet}$) or empty ($\ket{\circ}$) and evolve under processes of spontaneous decay and branching. The CCP, i.e. a contact process with classical branching [Fig. 1(a)], can be represented in the Lindblad formalism, thus representing the CCP as an OQS. The evolution of the CCP is then given by a Lindblad master equation with a zero Hamiltonian and dissipative contribution given by two terms,

[TABLE]

For a system of size $L$, the dissipative term $\mathcal{D}$ is defined as,

[TABLE]

where $\sigma_{-}\ket{\bullet}=\ket{\circ}$, $\sigma_{-}\ket{\circ}=0$, $\sigma_{+}=(\sigma_{-})^{\dagger}$, $n$ is the excitation density operator, $n\ket{\bullet}=\ket{\bullet}$, $n\ket{\circ}=0$, and the parameter $\gamma$ sets the strength of the spontaneous decay. The term $\mathcal{D}_{\rm br}$ instead describes the classical branching/coagulation process and can be written as,

[TABLE]

where $J^{(k)}_{\rm L}=\sigma_{1}^{(k)}n^{(k+1)},J^{(k)}_{\rm R}=\sigma_{1}^{(k+1)}n^{(k)}$ with $\sigma_{1}\ket{\bullet/\circ}=\ket{\circ/\bullet}$.

In the limit $L\to\infty$, the CCP exhibits a NEPT between a unique steady state devoid of particles and a degenerate steady state with finite particle density Hinrichsen (2000). The former steady state is known as the absorbing state, as once it is reached during the dynamics it cannot be escaped, and the corresponding phase is called the absorbing or inactive phase.

Absorbing phases occur when the strength of the decay process, characterised by the parameter $\gamma$, is sufficient to overcome the process of branching, characterised by the parameter $\Gamma$. Conversely, when the process of branching is strong enough relative to decay, i.e. the dimensionless branching is large, $\Gamma/\gamma\gg 1$, a steady state with finite particle density can been found, in what is known as the active phase.

In the QCP, branching is coherent and represented by the Hamiltonian,

[TABLE]

The parameter $\Omega$ thus sets the strength of coherent evolution and, as a consequence, of the branching/coagulation process. The relevant dimensionless parameter is then $\Omega/\gamma$, which determines the relative strength of coherent evolution, tending to create excitations, and dissipation, represented by $\mathcal{D}$, tending to remove them.

As with the CCP, in the limit $L\to\infty$, the QCP displays a non-equilibrium absorbing state phase transition Marcuzzi et al. (2016); Buchhold et al. (2017); Roscher et al. (2018). In fact, below the critical points, $\Omega<\Omega_{c}$ and $\Gamma<\Gamma_{c}$, in both models the stationary state is the same unique absorbing state, $\rho_{\text{a}}=\bigotimes_{k}\ket{\circ}\bra{\circ}_{k}$. Above the critical point, instead, the stationary state is degenerate and has a finite density of particles. Thus, both models have an active phase, though the stationary state in this phase will differ between the CCP and the QCP.

II.3 Non-equilibrium Setting

The specific setting we investigate is that of an initial single seed state evolving under the contact process dynamics Henkel et al. (2008), see Fig. 1. In this scenario, the initial state of the system is a product state. This state is unoccupied at all sites except the central one, $x_{\text{seed}}=\text{floor}(L/2)+1$. Therefore, the initial state can be written as, $\rho(0)=\sigma_{+}^{(x_{\text{seed}})}\rho_{\text{a}}\sigma_{-}^{(x_{\text{seed}})}$.

Starting from an initial seed state, the contact process can be characterised by the survival probability, $P_{\text{sur}}$, total density, $N_{\text{a}}$, and seed-site density, $n_{\text{seed}}$, defined as:

[TABLE]

where $n(t,k)=\operatorname{Tr}\left[\rho(t)n^{\left(k\right)}\right]$ is the density profile.

As illustrated in Figs. 1(c) and 1(d), in the absorbing phase, all clusters generated from a single seed die out so that $P_{\text{sur}}(t),N_{a}(t)$ and $n_{\text{seed}}(t)$ all tend to zero as $t\to\infty$. In contrast, within the active phase in the limit $L\to\infty$, these observables tend to non-zero values as $t\to\infty$. At criticality, the observables are characterised by universal power-law behaviour, which defines the exponents $\delta$, $\Theta$ and $z$ as:

[TABLE]

II.4 Universality Classes of the Contact Processes

While the CCP has been established to belong to the $1d$ DP universality class Hinrichsen (2000), the universal properties of the QCP – in particular the values of exponents and to which class the QCP belongs – is under debate. However, there have been a number of studies performed from different perspectives. For instance, the QCP has been examined from the perspective of mean-field theory Marcuzzi et al. (2016); Jo et al. (2019), functional renormalisation group Buchhold et al. (2017); Roscher et al. (2018) and tensor networks Carollo et al. (2019). The transition has been argued to be continuous Roscher et al. (2018), with the critical point being estimated as $\Omega_{c}\approx 6\gamma$, Carollo et al. (2019). A selection of critical exponents for the QCP have been estimated in Carollo et al. (2019) as well as in Buchhold et al. (2017), where the latter includes the effects of classical branching. Table 1 collects a number of these estimated exponents for the QCP, as well as those of the DP universality classes for comparison. We also emphasise that the QCP can in principle be experimentally realised in Rydberg quantum simulators Marcuzzi et al. (2016); Bloch et al. (2012); Bernien et al. (2017); Kim et al. (2018); Barredo et al. (2018), providing a possible check for theoretical results in the future.

III Time-Evolving Block Decimation in Closed Quantum Systems

The numerical simulations we perform are based on MPSs and the Time-Evolving Block Decimation (TEBD) algorithm, Vidal (2004), which we briefly summarise in this section. The TEBD algorithm has been applied extensively in closed quantum systems, on which we focus.

We consider many-body quantum systems (in $1d$) made of $L$ equal components with single-site Hilbert space having dimension $K$. In TEBD, the state $\ket{\psi}$ of such a system is represented as an MPS Schollwöck (2011),

[TABLE]

where $\mathbf{M}_{s_{k}}^{[k]}$ are matrices of size $\chi_{k-1}\times\chi_{k}$. The total number of parameters in the MPS representation is $\mathcal{O}\left(LK\chi^{2}\right)$, where $\chi$ is the maximum matrix dimension across the system and is known as the bond dimension. In principle, any state in the many-body Hilbert-space $\mathcal{H}=\mathbb{C}^{K\,\otimes L}$, can be represented exactly as an MPS with $\chi\leq K^{\text{floor}(L/2)}$.

In general, the exact representation of a state as an MPS is not possible, as the number of parameters increases exponentially with $L$. Nonetheless, frequently, relevant quantum states belong to set of states which can be represented as an MPS with $\chi=\mathcal{O}\left(\text{poly}\left(L\right)\right)$, or approximated to some arbitrary accuracy by such an MPS Eisert (2013). These are said to be represented/approximated efficiently by MPSs.

As a trivial example, product states can be represented by MPSs with $\chi=1$. More generally, the $L$ dependence of the von Neumann entanglement entropy,

[TABLE]

where $\varrho_{A/B}$ is the reduced density matrix of subsystem $A/B$ generated by a bipartition of the system,

[TABLE]

can be used to characterise the efficiency of an MPS representation. In one spatial dimension, an entanglement area law, $S\sim L^{0}$, for a given state suggests that the state can be efficiently represented as an MPS. This is indeed true for the ground-states of gapped local Hamiltonians, Hastings (2007), or states with exponentially decaying correlations, Brandão and Horodecki (2015). However, generally the scaling of $S$ alone is not enough to establish the accuracy of an efficient MPS approximation (in fact all Renyi entropies with index $\alpha<1$ must be used Schuch et al. (2008); Huang (2019)). Nonetheless, the entanglement entropy is useful in practice, particularly as a number of physical systems have been shown to obey area laws Eisert et al. (2010).

Assuming the initial state of the system to be one represented efficiently by an MPS, then an MPS representation of the time-evolved state can be found by applying some set of operators Schollwöck (2011), typically those constituting a Trotter decomposition of a quantum dynamics Hatano and Suzuki (2005). When these operators are applied, the resulting MPS will generally have a larger bond dimension than before. In fact, over time this will lead to an exponential increase in the required value of $\chi$ and generally time-evolution can only be treated exactly with MPS for short times Osborne (2006). This can be linked to the scaling of $S$ with time: If the entanglement entropy is growing linearly with time, as is the case in common quantum quench scenarios Calabrese and Cardy (2007), an efficient MPS approximation of the exact state is impossible Schuch et al. (2008).

Given the build-up of bond-dimension in an MPS representation over time, the key to performing time-evolution with MPSs is to repeatedly approximate the time-evolved MPSs, thus keeping the total number of parameters under control.

To achieve this, consider the MPS approximation to the state at time $t$, $\ket{\psi_{t}}$, and assume it has a bond dimension $\chi$, which is the maximum we will allow. To approximate the state at time $t^{\prime}$, we then apply the operator, $\hat{O}$, so that $\ket{\psi_{t^{\prime}}}=\hat{O}\ket{\psi_{t}}$. The new state $\ket{\psi_{t^{\prime}}}$ will now have some higher bond-dimension $\chi^{\prime}$, beyond the maximum we allow in our simulation. To remedy this, we want to find an MPS approximation to $\ket{\psi_{t^{\prime}}}$ that has bond dimension $\chi$. Calling this approximation $\ket{\phi_{t^{\prime}}}$, we then want to solve the minimisation problem,

[TABLE]

where $|\cdot|$ indicates the Hilbert-space norm. One can iterate this procedure to produce an approximation to the time-evolution of the state, allowing one to calculate desired observables along the way Schollwöck (2011).

The TEBD algorithm is, in essence, a simple approximation to the solution of $\left(\ref{eqn:truncation_problem}\right)$ given by considering successive bipartitions of the system at $k=1,2,...,L$. At each cut, one performs a Schmidt decomposition of the state $\ket{\psi_{t^{\prime}}}$ and discards a sufficient number of the smallest Schmidt coefficients so as to reduce the bond-dimension to $\chi$. More specifically, across the cut at site-$k$ the Schmidt coefficients, $\lambda_{j}^{[k]}$, are placed into a non-ascending order $\lambda_{1}^{[k]}\geq\lambda_{2}^{[k]}\geq...>0$. The MPS approximation with bond-dimension $\chi$ at site $k$ is then constructed by retaining the $\chi$ largest Schmidt coefficients. In other words, all Schmidt coefficients with $j\geq\chi+1$ are discarded, and the corresponding discarded weight,

[TABLE]

measures the error in this approximation. At any given cut, this approximation is optimal and so overall this provides a simple approximation to the solution of (14). We will use this method throughout, though we remark that many more sophisticated variations exist Paeckel et al. (2019).

IV Simulation of Universal Dynamics in the Double-Space

IV.1 Double-space Representation of Lindblad Dynamics

Perhaps the most straightforward way to apply TEBD to the study of Lindblad dynamics, (1), is to represent the density matrix $\rho(t)$ as a vector, $\ket{\rho(t)}$, in the “double-space” defined via the Choi-isomorphism, $\ket{n}\bra{m}\to\ket{n}\otimes\ket{m}$, Choi (1975). One thus has,

[TABLE]

Mapped in this way, the evolution of the quantum state can then be shown to be generated by the following Schrödinger-like equation ,

[TABLE]

where $\mathbb{L}$ is the representation of the Lindblad map (1) in the double space. It is possible to show that one must have,

[TABLE]

where $H_{D}$ is Hermitian and has the form,

[TABLE]

while,

[TABLE]

where ∗ means complex conjugation and T matrix transposition.

Solutions to (16) give the evolved state up to time $t$, $\ket{\rho(t)}=e^{t\mathbb{L}}\ket{\rho(0)}$, where $\ket{\rho(0)}$ is the initial condition for the density matrix in the vectorized representation. By performing a Trotter decomposition of the time-evolution operator $e^{t\mathbb{L}}$, which we choose to be a second-order scheme, one can apply the TEBD algorithm naturally to approximate $\ket{\rho(t)}$ using an MPS with bond-dimension $\chi$. From the approximation of $\ket{\rho(t)}$, observables can then be calculated as

[TABLE]

where $\hat{O}_{D}=\hat{O}\otimes\mathds{1}$ and $\ket{\mathds{1}}$ is the double-space state representation of the identity operator.

When considering the MPS representation of $\ket{\rho(t)}$, a natural quantity to consider is the operator space entanglement,

[TABLE]

with

[TABLE]

The value of $\tilde{S}$ plays an analogous role to the von Neumann entropy, $S$, for closed quantum systems and provides a characterisation of the computational difficulty of TN simulations.

IV.2 Schrödinger Picture Results

The evolution of the total density, survival probability and seed-site density are shown for the CCP and QCP in Fig. 2. In both cases the same TEBD algorithm is used, with a fixed Trotter step of $\gamma\delta t=0.1$ and bond-dimensions of $\chi=16,32,64$ and $256,512,1024$ for the CCP and QCP respectively. The simulations are performed with $\gamma=1$ at the estimated critical points, $\Gamma=6.75\gamma$ and $\Omega=6\gamma$ for the CCP and QCP respectively. For the CCP the critical point was estimated by scanning various values of $\Gamma$ and finding where both the total density and the survival probability show little deviation from a straight line in a log-log plot. In the case of the QCP, we take the previously estimated critical value of $\Omega=6\gamma$, Carollo et al. (2019). Both sets of observables show the correct qualitative behaviour, in line with the expectations of the critical dynamics, (9) - (11).

To approximate the critical exponents, power-law fits were performed to $N_{\text{a}}(t)$, $P_{\text{sur}}(t)$ and $n_{\text{seed}}(t)$ thus estimating $\Theta,\delta$ and $z$. To provide a best-estimate of these values, the simulations with largest $\chi$ were used for the fits, while the absolute differences between these estimates and those obtained by fitting to a $\chi$ of half the maximum were used for the error estimates.

In both the CCP and QCP, while all the different bond-dimension simulations agree at early times, at later times the low bond-dimension runs deviate considerably. This suggests that finite-bond effects can lead to a significant build-up of errors in observables, as in the closed quantum system case, even for classical states. However, the bond-dimensions needed to reach convergence until $\gamma t=20$ are very modest for the CCP. Since the fits for the CCP performed over $\gamma t\in[5,10]$ lead to estimated exponents within a few percent of the true $1d$ DP values, we see that critical dynamics of the CCP can be accurately simulated with MPSs in the double-space, and critical exponents estimated with small errors.

Compared with the CCP, simulation of the QCP requires much larger values of $\chi$ to achieve convergence in the examined observables, with $\chi=1024$ showing large deviations from a power-law by $\gamma t=10$. As such, fitting from $\gamma t\in[2,4]$ was used to establish these exponents, where the $\chi=1024$ simulations closely follow a power-law. As this is a relatively early time for which to perform the fits, these exponents may well contain finite-time errors. However, for comparison, fitting the CCP from $\gamma t\in[1,5]$ as opposed to $\gamma t\in[5,10]$ only increases the errors in the estimated exponents from around $1\%$ to $5\%$.

IV.3 Operator Space Entanglement in the Schrödinger Picture

To understand why the QCP is much harder to simulate than the CCP, we can compare the evolution of the operator space entanglement entropy, $\tilde{S}$. Taking the maximum value of the entanglement entropy across all bipartitions throughout, the evolutions of $\tilde{S}$ for the CCP and QCP are shown in Fig. 3(a).

In both the CCP and QCP, $\tilde{S}$ shows a clear “barrier” behaviour, where initially it grows rapidly to a peak around $\gamma t=0.5$ before decaying to a lower final value. This is consistent with an initial period of branching/coagulation evolution, where correlations build up rapidly and spontaneous decay is irrelevant, followed by a period where the latter becomes relevant and removes correlations/excitations from the state. While the overall picture seems the same for both the CCP and QCP, in the classical case the barrier is clearly much lower than in the QCP. Given that the operator space entanglement entropy should characterise the error in simulations, the difference between $\tilde{S}$ in the CCP and QCP helps explain the difference in accuracies found in observables. The relationship between $\tilde{S}$ and errors is illustrated in Fig. 3(b), which shows the error in the simulations of the QCP over time, equal to the square root of the discarded weight defined in (15). As the value of $\chi$ is increased, the error at each time drops, approximately halving when the bond-dimension is doubled. As with the entropy, the error shows a peak-like structure, with the peak occurring shortly after that of the entropy.

V Heisenberg Picture in the Double-space

To study the Heisenberg picture dynamics with TEBD in the double-space, one takes a representation of an operator $O$ as a vector $\ket{O}$, then evolves it through the dual Lindbladian $\mathbb{L}^{\dagger}$, such that

[TABLE]

In the case of closed quantum systems, the Heisenberg picture can be used to extend the maximum time over which simulations are accurate by roughly a factor of two Paeckel et al. (2019). Intuitively, this is expected because the dual dynamics of the unitary evolution, $U^{\dagger}(t)$, is equal to the original dynamics but backwards in time, $U^{\dagger}(t)=U(-t)$. Thus, one might expect that the dual dynamics is not radically different in terms of computational difficulty than that of the usual dynamics, and both maybe accurately approximated up to the same time, thus doubling the maximum time when combined.

In contrast, for open quantum systems, the dynamics implemented by $\mathbb{L}^{\dagger}$ is in principle completely different from the one of $\mathbb{L}$. This is exemplified by the case of the QCP, where the dual dynamics does not have an absorbing state. It has previously been suggested that the absorbing state, which is the unique steady state for any finite-size system, makes the application of tensor network methods more challenging Carollo et al. (2019). Therefore, one might expect that the absence of the absorbing state in the dual dynamics will mean that applying the Heisenberg picture is more accurate than the Schrödinger picture for the QCP.

To explore this, we calculate the survival probability for the QCP in the Heisenberg picture using TEBD and $\chi=256,512$, along with the operator-space entanglement entropy, as shown in Fig. 4. As can be seen, the entropy displays a characteristic barrier as in the Schrödinger picture. However, the barrier is substantially lower for the Heisenberg picture than for the Schrödinger picture case (though it is also less sharply peaked). Correspondingly, the survival probability shows dramatically reduced finite bond-dimension effects, with the $\chi=512$ approximation showing reasonable power-law behaviour up until $\gamma t=10$, leading to an estimated exponent of $\delta=0.27\pm 0.04$ when fit over $\gamma t\in[2,4]$.

On a practical level, these results suggest that the Heisenberg picture might allow for a more accurate approximation of the survival probability, $P_{\text{sur}}$, and thus $\delta$, at cheaper computational cost (lower bond-dimension). Other observables such as the total density can also be approximated in the Heisenberg picture, e.g., in order to establish the exponent $\Theta$ with greater accuracy, though we do not discuss this direction further.

VI Trajectories for the QCP

VI.1 Entanglement Distribution

As an alternative to the Lindblad formalism and double-space approach, we now consider a stochastic unravelling of the master equation realised algorithmically by Quantum Jump Monte Carlo (QJMC) Plenio and Knight (1998); Gardiner and Zoller (2004). In this QTs approach, each individual quantum trajectory corresponds to a pure state evolution. Therefore, the standard TEBD method for closed quantum systems can be applied quite directly to simulate a given trajectory, with the sample means over all trajectories providing an approximation for the observables of the average (Lindblad) dynamics.

As before, we will be interested in quantifying the presence of entanglement in the dynamics and how this affects the accuracy of TN calculations of universal non-equilibrium physics. Since each individual trajectory in QJMC is governed by a pure state evolution, the von Neumann entanglement entropy of the state has the usual physical meaning. However, unlike the closed system case, the entanglement of a trajectory, $S_{\text{traj}}$, can now be considered a random variable and thus associated to a distribution/probability. While this might seem to be a complicating feature compared to the single operator space entanglement entropy in the double-space case, it is actually very helpful for building a picture of the accuracy of simulations using MPSs: If we associate to a given bond-dimension, $\chi$, some characteristic maximum entanglement “cutoff”, $\bar{S}(\chi)$, then we expect that the set of trajectories with $S_{\text{traj}}\ll\bar{S}(\chi)$ will be well approximated, while those with $S_{\text{traj}}\approx\bar{S}(\chi)$ will be subject to finite bond-dimension effects. In other words, for low-entanglement trajectories the entanglement cutoff will be irrelevant but for high-entanglement trajectories it will be relevant.

To explore these issues, we first consider the distribution of $S_{\text{traj}}$ in $1000$ trajectories simulated with TEBD for $\chi=16,128$ and $\chi=256$. The evolution of each trajectory was calculated up until $\gamma t=10$ using a Trotter-step of $\gamma\delta t=0.01$, chosen so that the order of the associated error in the QJMC (which is a first-order scheme) is comparable with the second-order scheme used in the double-space case. The entanglement distributions for $\gamma t=1,5$ and $10$ are shown in Fig. 5, with the rows illustrating the evolution of the entanglement distributions for $\chi=16,128$ and $\chi=256$ respectively. Details of the histogram construction are given in the caption.

In Fig. 5(a), when $\gamma t=1$, all the histograms displayed are very similar; they show the same bimodality with one peak at $S=0$ - attributable to the unentangled absorbing state - and are similar except for the fact that the lower bond-dimension distributions display slightly more weight near $S=0$. As such, the means of the distributions, $\langle S\rangle$, shown as vertical dashed lines, are very close. This similarity can be explained by comparing these distributions with the maximum value of $S$ found for any trajectory at any time. These values, plotted as dotted vertical lines, can be considered as a measurement of the cutoff, $\bar{S}$. For all $\chi=16,128$ and $256$, the values of $\bar{S}$ lie well above the support of the distribution, and we can expect the effect of the cutoff to be minimal. In fact, given the separation between the support and $\bar{S}$, we might expect there to be essentially no effect and it is interesting that there is still a clear discrepancy at $S=0$. This discrepancy again suggests that the presence of an absorbing state poses challenges for numerical simulations.

In Fig. 5(b), where $\gamma t=5$, the distribution of entropies for $\chi=16$ and $\chi=256$ differ significantly, with a large proportion of trajectories for $\chi=256$ displaying entropies larger than the entanglement cutoff for $\chi=16$, leading to an artificial build-up of weight around $\bar{S}$ for the $\chi=16$ simulations. Thus, we can conclude that the finite-entanglement cutoff is relevant at this time for $\chi=16$, and the effect on the mean is clearly visible. This is in contrast to the case of $\chi=128$, which retains a good agreement with the $\chi=256$ simulations.

In Fig. 5(c), where $\gamma t=10$, the support of the distribution for $\chi=256$ is close to the entanglement cutoff $\bar{S}(\chi=256)$, and above $\bar{S}(\chi=128)$. This corresponds to a noticeable difference in the entanglement distributions for $\chi=128$ and $\chi=256$, though the difference is only slight compared to that of $\chi=16$. In fact, since only a small weight appears in the $\chi=256$ distribution above the value of $\bar{S}(\chi=128)$, we might expect that the accuracy of the $\chi=128$ simulations for observables will still be reasonably good. Moreover, we might expect that the $\chi=256$ simulations themselves are accurate due to the fact there is no significant build-up of weight near $\bar{S}(\chi=256)$, which for the $\chi=16$ case provided a clear indication of the entanglement cutoff’s relevance.

VI.2 Universal Dynamics with Trajectories

To assess the accuracy of the trajectories approach for the QCP and investigate the relationship between the entanglement distributions and the approximation of observables, we repeat the analysis for the QCP performed in Section IV, as shown in Fig. 6 for $\chi=64,128$ and $256$. From analysis of the entanglement distributions, Fig. 5, we expect to find good accuracy up to $\gamma t=10$ for $\chi=128$ and $256$, and indeed the curves of $N_{\text{a}}(t),P_{\text{sur}}(t)$ and $n_{\text{seed}}(t)$ overlap closely for these bond-dimensions, with $\chi=64$ deviating more noticeably.

Compared to the double-space simulations, Fig. 2, the observables calculated from trajectories seem much more accurate; with $\chi=128$ and $256$ they converge up to $\gamma t=10$ within statistical errors (given as two standard errors from the sample-mean and indicated by the shaded region in plots). Furthermore, they display the expected power-law behaviours, (9) – (11), for this whole period. This allows for fits to be performed over a longer region of time, $\gamma t\in[2,10]$, thus helping to eliminate finite-time errors. Since the curves for $\chi=64,128$ and $\chi=256$ lie well within the same shaded region for each observable, we can conclude that the finite bond-dimension effects are essentially negligible relative to the statistical errors in this region. As such, we provide purely statistical error estimates on the estimated exponents, obtained via a statistical bootstrap (see Fig. 6 for details). While these error estimates are large (corresponding to an approximate $95\%$ confidence interval), they can easily be reduced by increasing the number of samples.

VII Conclusions and Outlook

In this paper, we have made use of MPSs to study the critical dynamics of the classical and the quantum contact processes, which display non-equilibrium absorbing state phase transitions. For the QCP we have shown how the Heisenberg picture, where the dynamics does not display an absorbing state, can be used to improve the accuracy of simulations in the Lindblad formalism, over and above the Schrödinger picture. Moreover, when combining TEBD with a Quantum Jump Monte Carlo approach, we find that the expected critical behaviour can be reproduced with much higher accuracy for longer times than that of the Schrödinger picture Lindblad approach. In all approaches, we show that the entanglement in the MPSs can be used to understand these differences clearly, providing both a useful diagnostic tool and physical picture that links the numerical method to the dynamics in question.

The difference in accuracies found between the Lindblad and QTs approaches in the case of the critical QCP emphasises, as has been mentioned previously Jaschke et al. (2019), that when considering the simulation of open quantum systems with TN, one should examine different approaches to simulating the dynamics carefully. Given the observed superiority of a QTs approach for capturing the critical QCP dynamics, it would be interesting to know if this result is more general and whether approaches such as QJMC can allow one to study critical dynamics in other systems at higher accuracies than possible otherwise.

Finally, in the process of investigating these primary issues, we have also provided several results for the critical physics of the QCP. The most convincing conclusion that we can draw from these is that the universality class of the QCP cannot be DP, as evidenced by the fact that the best estimate of the exponent $\delta=0.26$ is far from that of $1d$ DP and $2d$ DP, see Fig. 6 and Table 1. This was also confirmed by the results from the double-space calculations, shown in Fig. 2. However, since that finite bond-dimension errors are small relative to statistical errors when using QJMC, we can use the statistical error to quantify the difference between exponents more carefully. With a standard error of $0.02$ for the estimate $\delta=0.26$, we have that the $1d$ DP value of $\delta=0.16$ lies five standard errors from the QCP estimate, while the $2d$ DP value of $\delta=0.45$ lies $9.5$ standard errors away. As such, it seems that the QCP universality class is genuinely different to directed percolation, though the presence of finite-time errors prevents us from stating this with absolute certainty.

Given these results, it is of interest to understand whether the QCP can be associated to some other known universality class and to identify exactly what the relevant quantum contributions might be that push QCP away from 1$d$ DP. An interesting remark in this regard can be made concerning the rapidity reversal symmetry present in DP, Henkel et al. (2008), which leads to the relation $\alpha=\delta$ between the two exponents characterising the decay of density, see Table 1. While we have only investigated the value of $\delta$ in this work, the value of $\alpha$ has been estimated previously in Carollo et al. (2019) to give $\alpha=0.36\pm 0.08$. This value lies $5$ standard errors from our estimated $\delta=0.26$. This seems to suggest that rapidity reversal is indeed broken in the QCP, though confirmation of this will require a better determination of $\alpha$, with QJMC offering a promising approach given the results presented here.

Acknowledgement

We thank H. Weimer and M. Buchhold for useful discussions. The research leading to these results has received funding from the European Research Council under the European Unions Seventh Framework Programme (FP/2007-2013)/ERC [grant agreement number 335266 (ESCQUMA)], the Engineering and Physical Sciences Council [grant numbers EP/M014266/1, EP/N03404X/1, EP/R04421X/1], and the Leverhulme Trust [grant number RPG-2018-181]. IL gratefully acknowledges funding through the Royal Society Wolfson Research Merit Award. We acknowledge the use of Athena at HPC Midlands+, which was funded by the EPSRC on grant EP/P020232/1, in this research, via the EPSRC RAP call of December 2018. We are grateful for access to the University of Nottingham’s Augusta HPC service.