Multispecies reaction–diffusion systems

A. Aghamohammadi, A. H. Fatollahi,
M. Khorrami, A. Shariati

Department of Physics, Alzahra University,
Tehran 19834, Iran.

Institute for Advanced Studies in Basic Sciences, P.O.Box 159, Gava Zang, Zanjan 45195, Iran.

Institute for Studies in Theoretical Physics and
Mathematics, P.O.Box 5531, Tehran 19395, Iran.

PACS numbers: 82.20.Mj, 02.50.Ga, 05.40.+j

Keywords: reaction–diffusion, multispecies

Multispecies reaction–diffusion systems, for which the time evolution equation of correlation functions become a closed set, are considered. A formal solution for the average densities is found. Some special interactions and the exact time dependence of the average densities in these cases are also studied. For the general case, the large time behaviour of the average densities has also been obtained.

## 1 Introduction

In recent years, reaction–diffusion systems have been studied by many people. As mean field techniques, generally do not give correct results for low dimensional systems, people are motivated to study stochastic models in low dimensions. Moreover, solving one dimensional systems should in principle be easier. Exact results for some models in a one–dimensional lattice have been obtained, for example in [1–10].

Different methods have been used to study these models, including analytical and asymptotic methods, mean field methods, and large-scale numerical methods. Systems with more than one species have also been studied [11–23]. Most of the arguments are based on simulation results. There are, however, some exact results as well ( [18, 20, 22] for example).

In [23], a 10–parameter family of stochastic models has been studied. In these models, the –point equal time correlation functions satisfy linear differential equations involving no higher–order correlations. These linear equations for the average density has been solved. But these set of equations can not be solved easily for higher order correlation functions. We have generalized the same idea to multi-species models. We have considered general reaction diffusion processes of multi-species in one dimension with two-site interaction. We have obtained the conditions the Hamiltonian should satisfy in order to give rise to closed set of time evolution equation for correlation functions. The set of equations for average densities can be written in terms of four matrices. The time evolution equation for more-point functions, besides these four matrices, generally depend explicitly on the elements of the Hamiltonian , and generally can not be solved easily. These matrices are not determined uniquely from the Hamiltonian: there is a kind of gauge transformation one can apply on them which of course, does not change the evolution equation. A formal solution for average densities of different species is found. For some special choices of the four matrices we also give the explicit form of interactions and the exact time dependence of average densities. At the end, we study the large time behaviour of the average densities of different species for the general case.

## 2 A brief review of linear stochastic systems

To fix the notation used in this article, here we briefly review the already well known formalism of linear stochastic systems. The master equation for is

(1) |

where is the transition rate from the configuration to . Introducing the state vector

(2) |

where the summation runs over all possible states of the system, one can write the above equation in the form

(3) |

where the matrix elements of are

(4) | |||

(5) |

The basis is dual to , that is

(6) |

The operator is called a Hamiltonian, and it is not necessarily hermitian. Conservation of probability,

(7) |

shows that

(8) |

where

(9) |

So, the sum of each column of , as a matrix, should be zero. As is a left eigenvector of with zero eigenvalue, has at least one right eigenvector with zero eigenvalue. This state corresponds to the steady state distribution of the system and it does not evolve in time. If the zero eigenvalue is degenerate, the steady state is not unique. The transition rates are non–negative, so the off–diagonal elements of the matrix are non–negative. Therefore, if a matrix has the following properties,

(10) |

then it can be considered as the generator of a stochastic process. It can be proved that the real part of the eigenvalues of any matrix with the above conditions is less than or equal to zero.

The dynamics of the state vectors (3) is given by

(11) |

and the expectation value of an observable is

(12) |

## 3 Models leading to closed set of evolution equations

The models which we address are multispecies reaction–diffusion models. That is, each site is a vacancy or has one particle. There are several kinds of particles, but at any time at most one kind can be present at each site. Suppose the interaction is between nearest neighbors, and the system is translationally invariant.

(13) |

The number of sites is and the number of possible states in a site is ; different states of each site are denoted by , where one of the states is vacancy. Introducing as the number operator of particle in the site , we have

(14) |

The average number density of the particle in the site at the time is

(15) |

where represents the state of the system at the time ,

(16) |

and

(17) |

So, the time evolution of is given by

(18) |

The only terms of the Hamiltonian which are relevant in the above equation are and . The result of acting any matrix on the ket is equivalent to acting the diagonal matrix on the same ket, provided each diagonal element of the matrix is the sum of all elements of the corresponding column in the matrix . So, the action of and on are equivalent to the action of two diagonal matrices on . We use the notation , for the equivalent action on .

(19) | |||||

(20) |

where and are as the following

(21) | |||||

(22) |

Then, equation (18) takes the following form

(23) |

Generally, in the time evolution equation of the two–point functions appear. Using (14), one can see that iff and satisfy the following equations, then the right hand side of the (23) can be expressed in terms of only one–point functions.

(24) | |||||

(25) |

These equations give constraints on the Hamiltonian, so adding the condition of stochasticity of , we have relations between the elements of . The constraints (24) mean

(26) | |||||

(27) |

So, (18) takes the form

(28) |

In the simplest case, the one-species, each site is vacant or occupied by only one kind of particles. Then, the matrices , , , and are two-dimensional. Using (14), The equation for is

(30) | |||||

This is a linear difference equation, of the kind obtained [23], and its solution can be expressed in terms of modified Bessel functions.

The time evolution equation for two-point functions also can be obtained.

(33) | |||||

(35) | |||||

For more–point functions, one can deduce similar results. In fact, it is easy to show that if the evolution equations of one–point functions are closed, the evolution equation of –point functions contain only - and less-point functions. However, generally these set of equations can not be solved easily.

## 4 Equivalent Hamiltonians regarding one-point functions, and gauge transformations

Knowing , , , and does not determine the Hamiltonian uniquely, but as it is seen from (28), the time evolution of one–point functions depends only on , , , and . The two– and more–point functions depend explicitly on the elements of . So, different Hamiltonians may give same evolutions for . Take two Hamiltonians and . Defining

(36) |

if

(37) |

these two Hamiltonians give rise to the same and . Regarding one–point functions , these models are the same. So, we call these models, regarding one–point functions, equivalent.

However, and do not determine , , , and uniquely. The stochastic condition

(38) |

results in some constraints on , , , and :

(39) | |||

(40) |

So, the sum of all elements of any column of () should be the same

(41) |

Then, the state is the left eigenvector of and , with the same eigenvalue . and have also the same property, of course with different eigenvalue .

Changing and according to the gauge transformation,

(42) | |||||

(43) |

does not change . With a suitable choice of :

(44) |

the sum of the elements of any column of or can be set to zero. In this gauge, the eigenvalues of and for the eigenvector will be zero.

## 5 One-point functions

To solve (28), we introduce the vector

(45) |

Equation (28) can then be written as

(46) |

Introducing the generating function ,

(47) |

one arrives at,

(48) |

the solution to which is

(49) |

’s are the coefficients of the Laurent expansion of , so

(50) |

This is the formal solution of the problem, which is of the form

(51) |

### 5.1 Some special cases

We now consider special choices for , , , and .

#### 5.1.1 The matrices , , , and are two–dimensional (the single–species case)

We can use the gauge transformation to make the simultaneous null left–eigenvector of , , , and . In this gauge, one has

(52) | |||||

(53) | |||||

(54) | |||||

(55) |

where

(56) |

This means that it is orthogonal to , and is a simultaneous right–eigenvector of , , , and . Using (52), one can easily calculate the exponential in (50):

(57) |

where

(58) |

and

(59) |

Now take and to be the left–eigenvector of dual to , and the right–eigenvector of dual to , respectively. One can normalize these, so that

(60) | |||||

(61) |

Of course, is orthogonal to . Then,

(62) |

Acting this on , and noting that

(63) |

it is seen that

(64) | |||||

(65) |

or

(66) |

This is equivalent to (30).

#### 5.1.2

Using (26)

(67) |

means that or , and or . If is not the left null eigenvector of and , then . So, we will have , and . Now, using the definition of

(68) | |||

(69) |

For , and , all the terms in the right hand side summations in the above equations are reaction rates and should be non–negative, but the sum of the left hand sides is zero. So,

(70) |

All the elements of each row except the diagonal elements of ( or ) are the same. That is,

(71) |

where is some diagonal matrix. The fact that is a left–eigenvector of , shows that it should be a left–eigenvector of as well. And this demands to be proportional to the unit matrix. One can do the same arguments for and . So, after gauge transformation,

(72) |

Although, the time evolution of average densities can be written in terms of , , , and , the Hamiltonian is not uniquely be determined by these matrices. There exist different Hamiltonians which are equivalent, regarding one–point functions.

(73) | |||

(74) |

All the elements of the column of are zero. For , the elements of satisfy

(75) | |||

(76) | |||

(77) | |||

(78) |

In general, these sets of equations have several solutions, but for the one–species case, the reaction rates are the following

(79) |

(80) |

The above system, with no diffusion, has been studied in [24]. There, the –point functions have been investigated. This solution can be generalized to the multispecies case. For

(81) | |||||

(82) | |||||

(83) |

The only constraint is the non–negativeness of the reaction rates:

(84) |

This model has free parameters. However, only the two parameters and appear in the time evolution equation of average densities:

(85) |

As it is seen, dynamics of average densities of different particles decouple, and despite the complex interactions of the model, ’s can be easily calculated. But in the time evolution of two-point functions ’s appear as well. So, although models with different exchanging rates () and same initial conditions have the same average densities, their two–point functions generally are not the same.

#### 5.1.3 commute

Generally, the gauge transformation do not preserve the commutation relation of and ( and that of and ). But if and commute, there is a gauge transformation which leaves the transformed and commuting. If we choose to be a right eigenvector of and dual to , that is

(86) |

then and commute. If

(87) |

then times and will be zero. So, if , , , and commute with each other, there exists a suitable gauge transformation that makes their eigenvalue corresponding to zero, while they remain commuting:

(88) |

Denote, the matrix which simultaneously diagonalize these four matrices by , diagonalized matrices by primes, and their eigenvalues by , , , and , respectively. We have

(89) | |||

(90) |

We take , and normalize and so that

(91) |

and

(92) |

will also diagonalize the exponential in (50). So we have

(93) |

where

(94) |

The matrix in the argument of the exponential in (93) is diagonal, so the integral can be easily calculated:

(95) |

Introducing , one arrives at

(96) |

which can be written in terms of modified Bessel functions

(97) |

Then,

(98) |

Note that the right–hand side of (97) is for , since the -th eigenvalue of , , , and is zero.

One can start with four special diagonal matrices, and then construct the Hamiltonians with different reaction–diffusion rates. Not all diagonal matrices lead to physical stochastic models: negative reaction rates may be obtained. Considering the large time behaviour of average number densities, one can show that

(99) |

which also shows that

(100) |

Now, we consider a special choice for :

(101) |

Then

(102) | |||||

(103) | |||||

(104) | |||||

(105) |

Now, consider

(106) |

For and ,

(107) |

So, taking and ,

(108) |

The same reasoning is true for and :

(109) |

Here too, similar to the previous example, the above choices for , , , and do not determine uniquely. One particular solution for the reaction rates is

For

(110) |

(111) |

and, for

(112) |

For , the following reactions may also occur. For

(113) |

and for

(114) |

The constraint of non–negativeness of the reaction rates leads to

(115) | |||||

(116) | |||||

(117) | |||||

(118) | |||||

(119) | |||||

(120) | |||||

(121) |

#### 5.1.4 type–change invariance

Suppose , , , and have the property

(122) |

and the same for the other three matrices. Note that the indices of these matrices are defined periodically, so that , as an index, is equivalent to . This is in fact a special case of commuting matrices discussed earlier. One can use (98). To do so, one should know the simultaneous eigenvectors of , , , and , and their corresponding eigenvalues. It is not difficult to see that the eigenvectors are

(123) |

The corresponding eigenvectors of , for example, are

(124) |

Finally, the matrix elements of the inverse of are

(125) |

These can be put directly in (98).

## 6 Large time behaviour of average densities

The large–time behaviour of the system is deduced through a steepest–descent analysis of the formal solution (50). One should consider the eigenvalues and the eigenvectors of the –dependent matrix

(126) |

Denote the eigenvalues of this matrix by . As for any value of , the matrix has as its left eigenvector corresponding to the eigenvalue zero, will have a right eigenvector dual to . is –dependent, but one can normalize it so that

(127) |

The fact that should not blow up at assures that the real–part of the eigenvalues of are non–positive (at least for ). If all of the other eigenvalues have negative real–parts, then at only survives. That is,

(128) | |||||

(129) |

We have used . This could also be obtained directly, using the evolution equation (28), by setting equal to zero and assuming independent of . So, the final state of the system is the eigenvector of , corresponding to the eigenvalue zero.

To investigate the next–to–leading term at