Open quantum reaction-diffusion dynamics: absorbing states and relaxation

Open quantum reaction-diffusion dynamics: absorbing states and relaxation

Merlijn van Horssen    Juan P. Garrahan School of Physics and Astronomy, University of Nottingham, Nottingham, NG7 2RD, UK
10th July 2019

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 and where particle hopping is coherent while reactions, such as pair annihilation or pair coalescence, are dissipative. These are quantum open generalisations of the and classical RD models. We characterise 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 non-equilibrium 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 and cases, may go asymptotically as with .

I Introduction

Non-equilibrium statistical mechanics is a field characterised by a richness of dynamical phenomena, many of which elude full theoretical understanding. One of the foremost theoretical challenges is the characterisation of non-equilibrium steady states (NESS) (and the relaxation towards them), which remains a topic of current research efforts Hinrichsen (2000); Zia and Schmittmann (2007); Peliti (2011). Several classes of systems (such as particle hopping models and directed percolation) exhibit phase transitions from fluctuating phases into a particular type of NESS, namely absorbing states; once reached, such absorbing states cannot be left. This behaviour is typical of reaction-diffusion (RD) models Toussaint and Wilczek (1983); Hinrichsen (2000); Tauber et al. (2005). If the diffusive mixing is not strong enough, asymptotic decay of global degrees of freedom can be slower than predicted by mean-field approximations Toussaint and Wilczek (1983); Hinrichsen (2000); Tauber et al. (2005), behaviour which is explained by fluctuation effects which have been confirmed experimentally Kroon et al. (1993); Allam et al. (2013), and continues to be a topic of current research Durang et al. (2014).

In this paper, we will consider a type of dynamics similar to RD models, extended to an open quantum spin chain. The theory of open quantum systems Buchleitner and Hornberger (2002); Gardiner and Zoller (2004) is the topic of ongoing current theoretical and experimental research in both quantum optics and cold atomic systems Diehl et al. (2008, 2010); Torre et al. (2013); Lee et al. (2012); Foss-Feig et al. (2013); Ates et al. (2012); Pichler et al. (2010). We consider a class of one-dimensional open quantum systems with dynamics analogous to that of classical RD models: particle propagation is coherent, through quantum hopping, while reactions between particles are dissipative. We show that, as in the classical RD models, these quantum systems exhibit non-exponential decay to absorbing stationary states, with mean-field approximations to the dynamics failing to predict the correct rate of decay. We connect this behaviour to the algebraic structure of the system Hilbert space Baumgartner and Narnhofer (2008) and present a classification of the dark states which generate the absorbing part of the dynamics – thus fully describing the stationary states (or NESS) for the associated quantum master operator. We also numerically study the dynamics of small systems via quantum jump Monte Carlo simulations Plenio and Knight (1998) and find evidence for a power law decay of the particle density towards the absorbing state, with an exponent which appears to be neither that predicted by mean-field analysis, nor that of the classical one-dimensional RD systems.

This article is organised as follows. We start by introducing the quantum reaction-diffusion models studied in this paper and discuss their immediate properties (explaining why they are sensible open quantum analogues of classical RD systems), by considering conservation of particle number and invariant subspaces, and by looking at their quantum jump trajectories. This first section concludes with a decomposition of the Hilbert space of these systems into transient (decaying) and recurrent (absorbing) subspaces. The next section contains our main analytical result as we consider the recurrent subspace in more detail. We provide a full classification of the dark states associated to our quantum RD models, and we consider the spectral properties of the quantum master operator to argue that these dark states generate the recurrent subspace. We will also explain the asymptotic behaviour of the evolution of the quantum state in terms of these dark states; we discuss how not all of the recurrent subspace is within reach from any initial state. We also provide an argument for the equivalence between annihilation and coagulation dynamics in our class of quantum RD models. The nontrivial asymptotic structure of the state resulting from the non-trivial collection of dark states is a fundamental feature of the quantum model, compared to classical RD dynamics. In the final section we characterise the decay of the density of particles. We show that a mean-field approximation to the dynamics, as in the classical RD case, predicts asymptotically a decay which is reaction-limited. In contrast, from quantum jump Monte Carlo simulations on small systems, we find evidence that in dimension one the dynamics is instead hopping-limited (cf. diffusion-limited) with a power law decay with an exponent smaller than the mean-field one, but apparently larger than that of the classical RD dynamics.

Ii Models

In this section we introduce the family of physical models considered in this paper, along with prerequisite mathematical background. We discuss features of the dynamics of the models by considering associated quantum jump trajectories, and we characterise the structure of the system Hilbert space according to conservation of the number of particles.

We consider quantum models which are analogous to classical one-dimensional RD models. Each of the models consist of a one-dimensional lattice where a site can be either empty, denoted by , or occupied by a single particle, denoted by . Particles can hop between lattice sites, symbolically . A reaction may only occur when two adjacent sites are occupied. We will consider three particular types of reaction: analogous to the classical reaction, we define pair annihilation, where two neighbouring particles annihilate, symbolically denoted by . We also define two analogous reactions to the classical reaction: asymmetric coagulation, where two neighbouring particles coalesce into one, denoted by ; and symmetric coagulation, where two neighbouring particles coalesce in two possible ways, . The quantum nature of our models is given by the fact that we will take particle hopping to be coherent, while reactions will be dissipative.

With periodic boundary conditions the Hilbert space of these models is that of a quantum spin chain of sites, . The coherent part of the dynamics is quantum hopping with rate , described by a Hamiltonian


where are Pauli operators acting on site , and due to the choice of periodic boundaries. For any one-site operator , the notation is employed as the usual shorthand for the tensor product , where denotes identity matrix. We work in a basis for in which the raising and lowering operators and take the form

with eigenvectors denoted by ; we denote the site number operator by .

The reactions are dissipative and the resulting dynamics of our model are those of an open quantum system. In the Markov approximation, the corresponding time evolution of the density matrix is determined by a quantum master equation (QME),

We use the notation for the quantum dynamical semigroup giving the solution to the QME as . The choice of jump operators determines the specific RD model. For the case of the annihilation reaction we have one operator for each pair of nearest neighbours, so that in one dimension there is one operator per site and the label coincides with the site label ,


The same occurs for asymmetric coagulation,


For the case of symmetric coagulation there are two jump operators per site since coagulation can be with the left or right neighbour,


With these definitions these quantum models are intuitively comparable to their classical counterparts. In particular, as we will see below, the QMEs associated to these three models (and therefore the dynamical properties) are equivalent in a well-defined sense.

ii.1 Invariant subspaces

One of the main points of interest of this paper is the time evolution of the density of particles. We define the density operator to be the global observable , whose expectation value is the density of particles.

The operators and the satisfy the algebraic relations


where the superscript (m) indicates any of the three models (3)-(5) above. The unravelling of the dynamics in terms of jump operators which generates the quantum jump trajectories, is such that the evolution between jumps is governed by the effective Hamiltonian . This operator coincides for the three models and takes the form


Due to the fact that , at the level of quantum jump trajectories associated to the models above, the density is a conserved quantity between jumps, only decreasing whenever a jump occurs. This is the first indication of the dynamical equivalence of the annihilation and coalescence models: both the evolution of the wave-function and the distribution of times between quantum jumps (for example in the quantum jump Monte Carlo unravelling of the dynamics Plenio and Knight (1998)) is determined by the which is the same for the three models. The only differences are in the post-jump states due to the distinct jump operators. We will return to this equivalence between the three models in the next section.

Conservation of the particle density between jumps allows us to define a natural decomposition of the system Hilbert space into subspaces with a given number of particles (see Fig. 2). We write


where the -particle subspace is characterised by for all states ; we will later consider this useful decomposition in more detail.

ii.2 Quantum jump trajectories

Figure 1: Sample quantum jump trajectories associated to the asymmetric coagulation model, for system size sites. The top three panels show individual site occupation as a function of time. The bottom panel is the particle density . The ratio characterises the relative relevance of reactions to diffusion. In all cases we expect an eventual crossover to a diffusion limited regime at low enough densities, which is most obvious for the case of : here an initial fast decay which exhibits an almost cutoff-like transition to a much slower relaxation.

We follow a quantum Markov chain Monte Carlo approach Plenio and Knight (1998) to simulate continuous-time quantum jump trajectories, which are realisations of the unraveling Wiseman and Diósi (2001) of the master equation in Eq. (II). In Fig. 1 we have plotted the site occupations and the density along such sample trajectories. Two important features, common to all three models, stand out in these trajectories. The first is that in between quantum jumps the wave function spreads due to the coherent hopping of particles, something that is more evident when the density is lower. The second feature is the decay of the density towards a stationary value. The rate of decay of is of particular interest to us, and one of our main results in this paper will be the characterisation of the power-law decay of this density.

We consider quantum jump trajectories with the completely filled initial state , with initial number of particles (so ). The number of particles subsequently decreases whenever a jump occurs, decreasing either by (in either of the coagulation models) or by (in the annihilation model). In each trajectory, the number of particles reaches a stationary value; as we will show below, the stationary values available to each trajectory are determined by the existence of stationary states within the sectors from Eq. (8). We show below that there are such dark states in the and sectors, and if is even, in the sector. The asymptotic values of are, depending on model and parity of , and ; in the limit of large system size, we of course have .

We note that in either of the coagulation models the sector is inaccessible when the trajectory is initialised outside of this sector; this is because at least two particles are required for a reaction to occur. The annihilation model conserves parity in the following sense: if is odd the reachable sector in a quantum jump trajectory is , whereas for even the sectors and are accessible.

A further comment relates to differences in the behaviour of the trajectories as we vary the diffusion rate relative to the reaction rate . The classical literature Hinrichsen (2000) distinguishes diffusion-limited and reaction-limited regimes depending on which of the two processes represents the limiting factor for the evolution of the particle density. This can refer both to the transient or to the asymptotic time behaviour. We already see something similar in the quantum jump trajectories of Fig. 1: starting from a fully occupied lattice, in all cases we see that there is a change in behaviour at higher densities where reactions are plentiful, to one at lower densities where particles need to propagate increasingly larger distances for reactions to occur. If we consider trajectories with different ratios , we see that the crossover between the two regimes occurs earlier for larger . As we will see below, there is a key distinction between the mean-field approximation for the dynamics of these systems and the actual behaviour in one dimension: the former predicts that the asymptotic relaxation is reaction-limited, while the latter is actually hopping-limited, as in the classical problem. (Aspects of the QME in the limit of low density are also investigated in Ref. Everest et al. (2014).) Having made these observations, in the remainder of this paper we will set the rates .

ii.3 Transient and recurrent subspaces

Figure 2: The decomposition of 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 . The corresponding evolution due to the quantum transition operator is sketched on the right; the stationary state will be a mixture of the states in (dark blocks) and stationary or non-decaying oscillating coherences between them (lighter blocks); shown in inset are absolute values of matrix elements of the actual density operator for large (obtained by integrating the QME), restricted to the one and two-particle subspaces (for details see Fig. 5).

The decomposition of the Hilbert space in Eq. (8) may be split up into a direct sum of a transient (or maximal decaying) subspace and a recurrent subspace. As indicated in Fig. 2, along a single quantum trajectory, the state jumps from higher into lower subspaces; collecting together all the subspaces that are eventually empty defines the transient subspace.

Following the terminology of Baumgartner and Narnhofer (2008), the higher levels can be regarded as a cascade of decay, with a natural subdivision into levels . As the state evolves, it moves strictly down the levels, that is, for any state and all ,

Indeed, following (Baumgartner and Narnhofer, 2008, Thm. 2) we write the system Hilbert space as a direct sum , where is the maximal decaying subspace. The latter is defined by requiring that for all initial states , and contains no further decaying subspace. In other words, the state of the system asymptotically becomes supported on the recurrent subspace .

Since the dissipative part of the dynamics consists of jump operators acting on two adjacent particles, it seems reasonable to suppose that the decaying subspace is given by the Hilbert space of states with two or more particles. However, as we will show below, there are dark states even in the two-particle subspace ; the recurrent subspace is generated by a fully known set of dark states, and in fact whenever is even (see Fig. 2).

Analysis of the eigenvalues of the master operator reveals the structure of the recurrent space beyond the stationary states (see Fig. 3). The peripheral eigenvalues of (i.e. those with vanishing real part) correspond to eigenstates in the recurrent subspace: they are either stationary (if the eigenvalue is ) or oscillating without decay. We will show below that each of the peripheral eigenvalues is associated to a specific dark state, thus completely characterising the recurrent subspace.

Figure 3: Spectrum of the transition operator for ; the numbers indicate the algebraic multiplicities of the peripheral eigenvalues (cf. Fig. 5).

Iii Dark states and the recurrent subspace

This section contains our main analytical results, arranged as follows: we first provide a complete classification of the dark states for the class of reaction-diffusion models described in the previous section. We consider the implications of these dark states for the structure of the recurrent subspace and we end the section with an argument for the equivalence of the three models. If we briefly return our attention to classical RD processes, we note that typically the recurrent subspace is fairly straightforward, consisting of only localised one-particle states which hop indefinitely. As we will see in this section, our quantum models present a more richly structured collection of absorbing states.

iii.1 Dark states

Recall that a vector is called a dark state Kraus et al. (2008) if it is an eigenvector of the Hamiltonian and belongs to the nullspace of all of the jump operators:


Note that we may immediate conclude that the dark states are the same for all three models since the nullspace of is independent of choice of model. We now proceed with one of our main results, which is a classification of the dark states of for any of the three models (3)-(5). We treat each of the subspaces separately: it is clear that the singleton subspace is dark. We will now separately consider the subspaces and (where we refer the reader to Appendices A and B for some of the calculations).

iii.1.1 One-particle subspace

The subspace of states with a single particle is contained in the nullspace of all the jump operators: whenever , we have . This means that the task of finding dark states in is reduced to that of diagonalising the restriction of the Hamiltonian to the subspace . The restricted Hamiltonian is a translation invariant hopping Hamiltonian,


in the usual position basis, corresponding to the action . We find (see Appendix A) that the eigenvalues of are given by


where . The corresponding eigenvectors are generated by the states


in the following sense: the largest eigenvalue is which has the unique eigenvector . If is even, is the smallest eigenvalue, and has the unique eigenvector .

All other eigenvalues are degenerate and have a pair of orthogonal eigenvectors associated to them: for (if is even) or (if is odd), the eigenspace associated to is two-dimensional, and is spanned by the vectors and .

We note the special case when is a multiple of , in which case the ground states of in are the eigenvectors associated to the zero eigenvalue , which take the simple form

It is clear from this discussion that is diagonalisable when restricted to ; the eigenstates , define an alternative basis for , sometimes referred to as the quasi-momentum basis. The subspace is entirely contained in the recurrent subspace ; below we will return to the dynamical consequences of this observation.

iii.1.2 Two-particle subspace

We now classify the dark states in the two-particle subspace . We first consider how the basis vectors for are generated; in particular, we express the position basis vectors in terms of the translation operator in order to exploit the translation invariance of the Hamiltonian. For example, the basis vectors for are generated by the single vector as

For the two-particle subspace , we denote (for ) by the vector . The vectors of this form generate all two-particle basis vectors : for ,

In Appendix B we use this expression for the basis vectors of to derive the equations that any dark state must satisfy. We find that if is odd, there are no dark states in .

However, when is even, there are dark states in , the number of which increases linearly in . In particular, the dark states are linear combinations of the states


and, if is itself even,


As we argue in the final section of the Appendix, these are the only dark states in .

In Tables 2 and 2 we summarise our classification of dark states, using the notation for the one-particle eigenvalues and quasimomentum states from Eqs. (11) and (12) and the notation for the two-particle dark states from Eqs. (13) and (14). To argue that there are no dark states to be found in any of the higher sectors, we provide numerical evidence in the remainder of this section.

Subspace Eigenvectors Eigenvalue Dark states
Table 2: Dark states, .
Subspace Eigenvectors Eigenvalue Dark states
( even)
Table 1: Dark states, .

iii.1.3 Numerical analysis

To further confirm the results summarised in Tables 2 and 2 and argue that there are no dark states found in any of the higher subspaces , , we now discuss related numerical results. We approach the problem in two ways: we perform an exhaustive search for dark states within each of the subspaces , and we exactly diagonalise the master operator, accounting for all the non-decaying eigenvalues solely using the dark states in Tables 2 and 2.

We search for dark states satisfying (9) in each of the subspaces . That is with


where and are the operators restricted to . We fully diagonalise and for each eigenvalue of , we look for a linear combination of the associated eigenvectors which is in the nullspace of all the jump operators. This means solving the system of linear equations with matrix of coefficients

where is the total number of jump operators ( for the pair-annihilation and asymmetric coagulation processes, for the symmetric coagulation process). The number of linearly independent dark states in is then given by the number of linearly independent solutions this system of equations, which is given by . This quantity is plotted in Fig. 4 for various , showing full agreement with the predicted number of dark states from Tables 2 and 2.

Figure 4: Number of dark states on a chain of sites predicted in Tables 2 and 2 by our analytical results, and found by numerical search.

For further confirmation of our results, and to clarify their role in the asymptotic behaviour of the dynamics, we consider the spectral properties of the master operator. Mixtures of dark states are eigenstates of corresponding to its peripheral eigenvalues – that is, eigenvalues of which are either or purely imaginary. To make this point clearer, suppose and are two dark states for corresponding to eigenvalues and of , respectively. Then , that is, is an eigenmatrix of corresponding to the peripheral eigenvalue : in general, if there are dark states, the sum of the algebraic multiplicities of the peripheral eigenvalues associated to dark states of must be . This, of course, does not rule out the existence of other stationary states of which are not pure; however, we argue along numerical lines that the peripheral eigenvalues are entirely accounted for by the eigenvalues obtained in this way from dark states.

Figure 5: Partial spectrum of restricted to lowest sectors with ; as predicted the total algebraic multiplicity of the peripheral eigenvalues is . Each of the peripheral eigenvalues is of the form where and are eigenvalues of associated to the dark states found in Tables 2 and 2.

Although, computationally speaking, full diagonalisation of for larger becomes prohibitively expensive, we are able to extract the peripheral eigenvalues by restricting the master operator to the lowest particle subspaces. We thus obtain a list of peripheral eigenvalues of along with their multiplicities; for example, for , the results are shown in Fig. 5. We find that, for values of up to , the total algebraic multiplicity of the peripheral eigenvalues is indeed , where is the number of dark states. Furthermore, if is one of the resulting peripheral eigenvalues of , there is a pair of eigenvalues associated to dark states such that .

This confirms that the peripheral eigenvalues are exactly determined by the dark states listed in Tables 2 and 2. Conversely, all peripheral eigenvalues are accounted for by exact diagonalisation in the sectors up to two particles; extending the diagonalisation scheme to higher sectors results in eigenvalues of with larger negative real part, as expected from the cascade of decay discussed in Fig. 2.

We briefly expand on the structure of the recurrent subspace . If is odd, the recurrent subspace consists only of , and the asymptotic dynamics is purely unitary on this subspace: as we have where This remains true when is even, but the recurrent subspace has an additional component in the subspace spanned by the two-particle dark states in Table 2.

In either case, dark states corresponding to the same eigenvalue of lead to pure stationary states for ; dark states corresponding to different eigenvalues of lead to oscillating coherences in the state , appearing in the off-diagonal blocks. In Fig. 6 we show how the matrix elements of the state for large reflects this structure. Note that these observations about the structure of serve to illustrate the general theory concerning the asymptotic structure of the time evolution of quantum dynamical semigroups found in e.g. Refs. Baumgartner and Narnhofer (2008); Albert and Jiang (2014).

Figure 6: Structure of the recurrent subspace appearing in the matrix elements of the asymptotic state obtained by numerical integration of the master equation for . The non-vanishing entries in 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).

iii.2 Reducibility of the dynamics and reachability of dark states

Our results concerning the dark states for this class of reaction-diffusion models allow us to conclude that the recurrent subspace is generated by the set of dark states. The question remains if it is possible to reach any state in the recurrent subspace by choosing an appropriate initial state – a process of interest in quantum control Mabuchi and Khaneja (2005), where it is known as quantum state engineering Verstraete et al. (2009). Aside from this, the question also relates to what occurs in the classical pair-annihilation model Hinrichsen (2000) where there is a conservation of parity, based on which we can always predict whether the stationary state will be or a -particle state.

Figure 7: Asymptotic behaviour of quantum jump trajectories in terms of dark states, obtained by taking final states of 200 trajectories and collecting identical states. The one- and two-particle dark subspaces are represented as convex sets with the dark states in Tables 2 and 2 as their extremal points. The reachable one-particle states are indicated in the interior of ; the reachable two-particle states appear randomly distributed in the interior of the two-particle dark subspace.

Although a full treatment of engineering of dark states is beyond the scope of this paper, we will briefly consider reachability Ying et al. (2013) of the dark states we have previously identified. In the current context, the question of reachability is as follows: given a state in the recurrent subspace (which is necessarily a mixture of dark states), is it possible to find an initial state in the transient subspace such that as ?

Our results indicate a negative answer to this question; that is, it appears that for pure states starting outside of the recurrent subspace , only a few superpositions of the one-particle dark states are available as stationary states. We approach this problem as follows: starting with the full initial state, we compute an ensemble of trajectories, each of which until it has reached stationarity – that is, until the quantum state is fully supported on . We then take each final state and compute its inner product with each of the dark states found in Tables 2 and 2. Recalling that the dark states form an orthogonal basis of , these inner products uniquely determine the final trajectory states.

We have plotted the magnitudes of the resulting coefficients in Fig. 7, for the case , taking an ensemble of trajectories. It appears that, whenever the initial state does not have any dark components, the final state state component in the one-particle space is limited to one of only a few possible states, while the component in the two-particle space is chosen randomly. Rephrasing in terms of reachability, the few one-particle dark states that appear asymptotically are reachable, while the others are not – their coefficients need to be present initially. In the case of Fig. 7 with , the reachable states in are

As indicated in Fig. 7, out of the final states, these occur , and times, respectively; the remaining states are randomly distributed elements of the two-particle dark space.

Repeating this experiment with a randomly chosen initial state vector (outside of ) confirms these results: the only states reached in are the . Based on this, we conjecture that not all dark states are reachable from a state which starts outside of the dark space. We conclude our observations on reachable states by noting that the coefficients which appear in the vectors (for general ) are exactly the coefficients of the dark states in Eq. (12); this is certainly an interesting property of the reachable states, and the question remains open whether this is an indication of some deeper symmetry in the dynamics.

iii.3 Equivalence of models

Figure 8: Comparison of density from quantum jump trajectory averages for the three models for . The coagulation models show a non-zero final density, while the density for the annihilation model decays to zero.

The classical reaction-diffusion processes have been shown to be equivalent in the sense that their Liouvillian operators are related by similarity transformations Hinrichsen (2000). This notion of equivalence may be extended to quantum systems by letting two quantum systems, defined by their master operators (or quantum Liouvillians), be equivalent if their master operators are related by a similarity transformation Burgarth and Yuasa (2012, 2014). This relation, in turn, may be defined via the Choi-Jamiolkowski isomorphism on the level of the operator-sum representations of the master operators, by requiring that their matrix representations are similar in the usual sense. That is, since the operator-sum representation is an isomorphism, the notion of matrix similarity between the matrix representations of linear maps immediately defines similarity of master operators: if and are two master operators we may say that and are similar if there exists an invertible linear map such that .

This allows us to call the quantum systems defined by master operators equivalent if their master operators are similar in this well-defined sense. Expected properties for equivalent systems are immediate from this definition: in particular, two master operators have equivalent if and only if they have exactly the same spectrum. If this is the case, the eigenstates are related by the similarity transformation.

For the three master operators , and defined in Eq. (II)-(5) we have plotted the density in Fig. 8. On the level of the expectation value of observables such as , a similarity transformation between the master operators would need to likewise transform between the densities . For evidence that the master operators are indeed equivalent, we have computed their eigenvalues and compared them for increasing system size; for , for example, see Fig. 9. Indeed, it appears that, as far as we are able to numerically diagonalise the master operators, the eigenvalues (and their algebraic multiplicities) agree, and the models are indeed equivalent in the sense discussed above.

Figure 9: (a) Comparison of the spectrum of for the three models, with overlapping symbols showing complete agreement (). (b) For clarity, instead of plotting all eigenvalues for larger we show the density of states as function of the real and imaginary part of the spectrum (); (c) , truncated to three particles; (d) , truncated to two particles.

Iv Decay of particle density

Having established the asymptotic behaviour of the three equivalent quantum reaction-diffusion models by classifying the dark states and the structure of the recurrent space, we now turn to the dynamics of the system and consider how the density decays. We start by working out the mean-field approximation of the site densities , which results in a power-law prediction for the density of the form . We will then perform an analysis of the actual behaviour of the density as obtained from simulations of the dynamics; it turns out that a power law decay still holds, but with a slower rate than predicted by the mean-field approximation.

iv.1 Mean-field approximation

Using the Gutzwiller product state ansatz we approximate the equations of motion for the expectation values of observables and :

where . Assuming a homogeneous state, so that one-site expectation values are all equal, we arrive at the mean-field equation of motion for the average density ,

This equation is readily solved to obtain a power-law decay with exponent ,


where denotes the initial mean site density. This is the same type of behaviour obtained in the mean-field predictions for the analogous classical reaction-diffusion processes Hinrichsen (2000). Note that the mean-field result also predicts that the asymptotic regime is reaction-limited as it only depends on . In the classical case, however, the mean-field law is only accurate above two dimensions Hinrichsen (2000); Tauber et al. (2005): in one dimension spatial fluctuations dominate, the asymptotic regime is diffusion-limited, and the density decays as . A natural question is whether an analogous departure from mean-field behaviour is observed in our one-dimensional quantum models.

iv.2 Numerics in one dimension

Figure 10: (a) Density for from quantum jump Monte Carlo, compared to power law fit. (b) Convergence of density for . (c) Convergence of the power law exponent (the plot shows ) with increasing system size (dependence on trajectory length (arb. units) is also indicated.)

From the quantum jump simulations we obtain the time evolution of the particle density . As noted in the previous section, the three models defined in Eq. (II) are equivalent; in particular, if a power law decay holds for the density in one of the models, it is true for all of them. We therefore mainly consider the pair-annihilation model in the following discussion, unless otherwise specified.

In Fig. 10(a) we have plotted the density for sites, obtained by taking an average of trajectories, on a log-log scale. After an initial transient period, the density seems to follow a power law. If we fit a function of the form to this data, we obtain , with a vanishing coefficient , which is expected in the pair-annihilation model for even and a fully occupied initial state: the two-particle dark states in the subspace are highly unlikely to be visited, and hence any trajectory almost always ends up in the zero-particle state. If is odd, the trajectory will always converge to the one-particle subspace, resulting in a non-vanishing coefficient .

This distinction between odd and even is clearly visible in Fig. 10(b), where we have plotted the trajectory average density for ranging from to . The even-numbered sites decay to while the odd-numbered sites approach a non-zero asymptotic value determining the coefficient . However, since the density computed on the subspace is , in the limit of large chain lengths the difference between odd and even vanishes.

In order to estimate the decay exponent we have analysed fits as the one shown in Fig. 10(a) for varying trajectory lengths for times between and , for system sizes up to . As shown in Fig. 10(c), the distinct final state reached for even or odd sites impinges on the fitted value of , with even giving exponents and odd giving exponents . This disparity between odd and even fits appears to diminish as the number of sites increases. Extrapolation to large of our simulations would suggest a coefficient in the range . As mentioned previously, the difference between the asymptotic values for , computed for odd and even , should vanish based on the dynamics alone. A better indication of the actual behaviour in this large system limit, however, is beyond our computational abilities due to the usual exponential-in-size complexity of quantum dynamics.

Our numerical results suggest that dynamics of the quantum RD models studied here are fluctuation dominated in one dimension, and therefore disagree with a mean-field prediction, in analogy to what occurs in the classical case. In particular, the decay exponent is smaller than the mean-field prediction, . The numerical results, albeit for small systems, suggest however that is larger than the one-dimensional classical exponent, . One can speculate that this is due to the fact that in the quantum case particle propagation is via a quantum random walk Aharonov et al. (1993) which explores space somewhat more efficiently than classical diffusion. This potential discrepancy between the classical and quantum decay exponents is an interesting point that would require more solid confirmation.

V Conclusions and outlook

Our aim in this paper was to understand which features, characteristic of single-species classical reaction-diffusion systems, remain when particle propagation becomes coherent. The open quantum systems we considered were natural generalisations of the classical and reaction-diffusion models. We have found that many of the features of the classical systems are also present in the quantum models. As expected the quantum systems were shown to display a slow relaxation towards absorbing states. In particular, we were able to provide a thorough classification of the recurrent (i.e.  absorbing) subspace in terms of the dark states of the quantum master operator that generates the dynamics. These absorbing subspaces have a much richer structure in the quantum models than in the classical ones, something that may be of potential interest from the point of view of quantum state preparation. We have also provided strong evidence for the dynamical equivalence of the annihilation and coalescence quantum models, in analogy with the equivalence of the corresponding classical models. This was done through the coincidence of the spectra of the dynamical super-operators from direct numerical diagonalisation of finite systems. It would be interesting to find a similarity transformation that maps the quantum models, to bring this equivalence to an exact footing as in the classical case. Our final result is the observation that in one dimension the relaxation of the density seems to follow a power law with an exponent which is smaller than that of the mean-field analysis. This is again analogous to what occurs in the classical models. In the quantum case, however, our numerics would suggest that the exponent is not the same as the classical one, which would be an indication of a distinctive quantum feature in the non-equilibrium dynamics. This interesting possibility deserves further investigation.

We thank E. Levi, M. Hush and J. Schick for helpful discussions. This work was supported by EPSRC grant no. EP/L50502X/1 and grant no. EP/J009776/1.


Appendix A One-particle dark states

In this Appendix we derive the eigenvectors and eigenvalues of the Hamiltonian restricted to the one-particle subspace , which, since the entire subspace is in the nullspace for all of the jump operators, are also dark states for the master operator .

Starting with the restricted Hamiltonian in Eq. (10), we note that it takes a particular form, which we exploit to characterise its eigenvectors. Recall that a complex matrix is called a circulant matrix Gray (2005) if it is of the form

a matrix completely specified by its initial row . It is easy to see that is a circulant matrix with initial row . The eigenvalues of the general circulant matrix are Gray (2005)

with corresponding eigenvectors

where for we denote by the -th roots of unity. Therefore, the eigenvalues of are given by

note that for all , resulting in distinct eigenvalues.

For the eigenvectors of , we start with the eigenvalue , which has a unique corresponding eigenvector . For each pair of coinciding eigenvalues , with , there exists a two-dimensional eigenspace generated by the vectors .

We note that for we have the identities and similarly . By taking the difference and the sum we obtain an alternative pair of vectors spanning the eigenspace,

We finally note that if is even, the eigenvalue is simple with unique eigenvector . The identification results in the expression for the eigenvectors in Eq. (12).

Appendix B Two-particle dark states

In this Appendix we derive the dark states in the two-particle subspace . We start by finding a convenient notation for the basis vectors for in terms of the translation operator . As discussed in the main text, the vectors of the form generate the basis vectors of .

To see how the basis vectors are generated, first suppose is odd, say . Then there are generating vectors : for we have

for , generating the basis vectors of .

Now suppose is even, say . Then there are generating vectors : for we have

for and with

for , generating the basis vectors of .

Note that the basis vectors generated by are exactly the pure states of the form with two adjacent particles. Since consists only of scalar multiples of , dark states in are the eigenvectors of which are superpositions of all basis states except those generated by . The remainder of this Appendix concerns the characterisation of these eigenvectors; we consider separately the cases when is odd and when is even. As we are only considering elements of in the following, we will omit the subscript on .

b.1 Odd

Let be odd with and let ; then using the position basis introduced above we may write

where . To determine if is an eigenvector of we first consider the action of on the basis vectors: we find that is given by


Using translation invariance of we obtain the expression

for the first term, we use Eq. (17) and, by relabelling indices and using the periodicity condition , find

Similarly, the second and third terms are given by


respectively. We thus arrive at the expression