Open quantum reaction-diffusion dynamics: Absorbing states and relaxation

We consider an extension of classical stochastic reaction-diffusion (RD) dynamics to open quantum systems. We study a class of models of hard-core particles on a one-dimensional lattice whose dynamics is generated by a quantum master operator. Particle hopping is coherent while reactions, such as pair annihilation or pair coalescence, are dissipative. These are quantum open generalizations of the $A+A\to ⌀$ and $A+A\to A$ classical RD models. We characterize the relaxation of the state towards the stationary regime via a decomposition of the system Hilbert space into transient and recurrent subspaces. We provide a complete classification of the structure of the recurrent subspace (and the nonequilibrium steady states) in terms of the dark states associated to the quantum master operator and its general spectral properties. We also show that, in one dimension, relaxation towards these absorbing dark states is slower than that predicted by a mean-field analysis due to fluctuation effects, in analogy with what occurs in classical RD systems. Numerical simulations of small systems suggest that the decay of the density in one dimension, in both the open quantum $A+A\to ⌀$ and $A+A\to A$ systems, behaves asymptotically as $t^{-b}$ with $1/2<b<1$.

Figure 1

Sample quantum jump trajectories associated to the asymmetric coagulation model for system size $N=22$ sites. The top three panels show individual site occupation $\langle n_i\rangle$ as a function of time. The bottom panel is the particle density $\langle \Lambda \left(t\right)\rangle$. The ratio $\nu :=\kappa /\Omega$ characterizes the relative strength of reactions to diffusion. There is an eventual crossover to a diffusion limited regime at low-enough densities, which is most obvious for $\nu \gg 1$: here an initial fast decay which exhibits an almost cutoff-like transition to a much slower relaxation.

Figure 2

The decomposition of $\mathcal{H}$ in Eq. (8) can be thought of as a cascade of decay for the quantum jump trajectory which ultimately ends up in the recurrent subspace $P_0\mathcal{H}$. The corresponding evolution of a density matrix (due to the quantum transition operator $\mathcal{T}:={\mathcal{T}}_1)$ is sketched on the right; the stationary state will be a mixture of the states in $P_0\mathcal{H}$ (dark blocks) and stationary or nondecaying oscillating coherences between them (lighter blocks); shown in the inset are absolute values of matrix elements of the actual density operator $\rho \left(t\right)$ for large $t$ (obtained by integrating the QME), restricted to the one- and two-particle subspaces (for details see Fig. 5).

Figure 3

Spectrum of the transition operator ${\mathcal{T}}_1=exp\mathbb{W}$ for $N=6$; the numbers indicate the algebraic multiplicities of the peripheral eigenvalues (cf. Fig. 5).

Figure 4

Number of dark states on a chain of $N$ sites predicted in Tables 1 and 2 by our analytical results and found by numerical search.

Figure 5

Partial spectrum of ${\mathcal{T}}_1=exp\mathbb{W}$ restricted to $k=2$ lowest sectors with $N=10$; as predicted the total algebraic multiplicity of the peripheral eigenvalues is ${14}^2=196$. Each of the peripheral eigenvalues is of the form $e^{i({\lambda }_i-{\lambda }_j)}$, where ${\lambda }_i$ and ${\lambda }_j$ are eigenvalues of $H$ associated to the dark states found in Tables 1 and 2.

Figure 6

Structure of the recurrent subspace $P_0\mathcal{H}$ appearing in the matrix elements of the asymptotic state $\rho \left(t\right)$ obtained by numerical integration of the master equation for $N=8$. The nonvanishing entries in $\rho \left(t\right)$ are generated by the dark states found in Table 2, which are indicated in this figure. The off-diagonal blocks are coherences between dark states, which oscillate at a rate determined by the difference of the corresponding eigenvalues. (Note: The highlighted regions are only for intuitive purposes; the actual dark states are spread out in the position basis.)

Figure 7

Asymptotic behavior of quantum jump trajectories in terms of dark states, obtained by taking final states of 200 trajectories and collecting identical states. The numbers (and colors) in main image indicate inner products of final trajectory states (vertical axis) with dark states (horizontal axis); the bottom three rows correspond to repeated final states, with their multiplicities indicated on the arrows. The one- and two-particle dark subspaces are represented as convex sets with the dark states in Tables 1 and 2 as their extremal points. The reachable one-particle states $|{\theta }_i\rangle$ are indicated in the interior of ${\mathcal{H}}_1$; the reachable two-particle states appear randomly distributed in the interior of the two-particle dark subspace.

Figure 8

Comparison of density $\langle \Lambda \left(t\right)\rangle$ from quantum jump trajectory averages for the three models for $N=12$. The coagulation models ${\langle \Lambda \left(t\right)\rangle }_{\text{(ac)}},{\langle \Lambda \left(t\right)\rangle }_{\text{(sc)}}$ show a nonzero final density, while the density for the annihilation model ${\langle \Lambda \left(t\right)\rangle }_{\text{(ann)}}$ decays to zero. As discussed below, these expectation values are necessarily related by a transformation induced by the the similarity transformation between the master operators.

Figure 9

(a) Comparison of the spectrum of $\mathbb{W}$ for the three models, with overlapping symbols showing complete agreement $\left(N=4\right)$. (b) For clarity, instead of plotting all eigenvalues for larger $N$ we show the density of states as function of the real and imaginary part of the spectrum $\left(N=6\right)$; (c) $N=7$, truncated to three particles; (d) $N=8$, truncated to two particles.

Figure 10

(a) Density $\langle \Lambda \left(t\right)\rangle$ for $N=20$ from quantum jump Monte Carlo, compared to power-law fit. (b) Comparison of density $\langle \Lambda \left(t\right)\rangle$ for $N=4,...,21$. (c) Convergence of the power-law exponent $b$ with increasing system size $N$ [dependence on trajectory length $T$ (arb. units) is also indicated]. With increasing $N$, oscillations between subsequent values of $b$ become less pronounced, and the value of $b$ in the large system limit appears to stabilize in the range $0.7<b<0.9$. We conclude that in the large system limit $b_{\text{class}\text{1d}}<b<b_{\text{mf}}$, where $b_{\text{class}\text{1d}}=1/2$ is the exact classical reaction-diffusion exponent and $b_{\text{mf}}=1$ is the quantum mean-field prediction for $b$, respectively (see text).