Numerical simulation of quantum nonequilibrium phase transitions without finite-size effects

Classical $(1+1)$-dimensional (D) cellular automata, as for instance Domany-Kinzel cellular automata, are paradigmatic systems for the study of nonequilibrium phenomena. Such systems evolve in discrete time steps, and are thus free of time-discretization errors. Moreover, they display nonequilibrium phase transitions which can be studied by simulating the evolution of an initial seed. At any finite time, this has support only on a finite light cone. Thus, essentially numerically exact simulations free of finite-size errors or boundary effects are possible, leading to high-accuracy estimates of critical exponents. Here, we show how similar advantages can be gained in the quantum regime: The many-body critical dynamics occurring in $(1+1)\mathrm{D}$ quantum cellular automata with an absorbing state can be studied directly on an infinite lattice when starting from seed initial conditions. This can be achieved efficiently by simulating the dynamics of an associated one-dimensional, nonunitary quantum cellular automaton using tensor networks. We apply our method to a model introduced recently and find accurate values for universal exponents, suggesting that this approach can be a powerful tool for precisely studying nonequilibrium universal physics in quantum systems.

Figure 1

Seed evolutions in a $(1+1)\mathrm{D}$ (quantum) cellular automaton. (a) Seed evolution of the classical Domany-Kinzel cellular automaton at the critical site-DP point [3] performed directly on an infinite lattice. Occupied sites only fall inside the indicated light cone (dashed lines and shaded region). (b) The total number of occupied sites, $N\left(t\right)$, is shown averaged over 1000 runs (solid black line). Even for relatively short times, the universal power law can be seen (straight solid red line), and the observed exponent is in agreement with the expected value for 1D DP: ${\theta }_{\text{DP}}=0.314$. (c) In a $(1+1)\mathrm{D}$ QCA, a 2D lattice is initiated in a product state with all empty sites apart from the first row, which encodes the initial condition. The state $|\psi \left(1\right)\rangle$, obtained by updating the target sites in the subsequent row via the application of three-body unitary gates, is shown. The operator ${\mathcal{G}}_1$ is the product of the applied gates, and thus $|\psi \left(1\right)\rangle ={\mathcal{G}}_1|\psi \left(0\right)\rangle$. (d) When $(1+1)\mathrm{D}$ QCA on an infinite lattice display a strict light cone, the reduced dynamics of a row with seed initial conditions are fully captured by those of a finite-size reduced state, $\rho \left(t\right)$. The size of this state grows proportionally with time. This allows for the dynamics of the underlying unitary $(1+1)\mathrm{D}$ QCA to be studied via the corresponding nonunitary dynamics of $\rho \left(t\right)=\mathrm{\Lambda }\left[\rho \right(t-1\left)\right]$, free of finite-size effects.

Figure 2

Reduced dynamics of the $(1+1)\mathrm{D}$ QCA. (a) To evolve $\rho (t-1)\to \rho \left(t\right)$ through the map $\mathrm{\Lambda }$, we begin via an MPO representation of $\rho (t-1)$, shown here for $t=2$. Empty sites are then added at locations where gates act nontrivially. This operation defines the state ${\mathrm{\Xi }}_{t-1}$. We have depicted these states here using the standard diagrammatic notation for TNs [14, 22]. In this notation, tensors are represented by shapes with a number of legs corresponding to their order. Each tensor corresponding to operators in a local Hilbert space must have two legs for the “physical” indices, here represented by a solid black circle. Other “virtual” legs which join the shapes (solid black lines) encode correlations between these, and we denote the trivial legs (indicating no correlations) as dashed lines. (b) The gates are then applied in MPO form to update the state. By tracing out the sites of row $t-1$, indicated diagrammatically by removing physical legs, the exact TN representation of $\rho \left(t\right)$ is obtained. Approximating this by an MPO allows the scheme to be iterated.

Figure 3

Critical behavior. (a) The evolution of the total number of occupied sites, $N\left(t\right)$, is shown for $t\in \left[1,100\right]$ for various values of $\mathrm{\Gamma }$. This includes five intermediate values close to the critical point—as indicated by the almost linear behavior in the log-log plot—and two values (the bottom-most and top-most lines) further from the critical point, illustrating the NEPT from a state of zero particles to one with a diverging number. In ascending order from the lowest line, the values of $\mathrm{\Gamma }$ are $0.98,0.995,0.996,0.997,0.998,0.999,1.01$. The inset shows $\theta \left(t\right)$ for the central five values. From these, the flattest curve provides the estimate for the critical point, ${\mathrm{\Gamma }}_c=0.997\pm 0.01$, and exponent, $\theta =0.307\pm 0.017$. The shaded region represents the error [35]. The estimated value of $\theta$ is consistent with that of 1D DP, indicated by the dotted green line. $\chi$ is the maximum MPO bond dimension used. (b) Corresponding plots for $\omega =1$ with $\mathrm{\Gamma }=1.015,1.03,1.032,1.034,1.035,1.04,1.05$. These produce the estimates ${\mathrm{\Gamma }}_c=1.034\pm 0.02$ and $\theta =0.32\pm 0.03$, also consistent with 1D DP.