On reaction-subdiffusion equations

To analyze possible generalizations of reaction-diffusion schemes for the case of subdiffusion we discuss a simple monomolecular conversion A -->B. We derive the corresponding kinetic equations for local A and B concentrations. Their form is rather unusual: The parameters of reaction influence the diffusion term in the equation for a component A, a consequence of the nonmarkovian nature of subdiffusion. The equation for a product contains a term which depends on the concentration of A at all previous times. Our discussion shows that reaction-subdiffusion equations may not resemble the corresponding reaction-diffusion ones and are not obtained by a trivial change of the diffusion operator for a subdiffusion one.

To analyze possible generalizations of reaction-diffusion schemes for the case of subdiffusion we discuss a simple monomolecular conversion A → B. We derive the corresponding kinetic equations for local A and B concentrations. Their form is rather unusual: The parameters of reaction influence the diffusion term in the equation for a component A, a consequence of the nonmarkovian nature of subdiffusion. The equation for a product contains a term which depends on the concentration of A at all previous times. Our discussion shows that reaction-subdiffusion equations may not resemble the corresponding reaction-diffusion ones and are not obtained by a trivial change of the diffusion operator for a subdiffusion one. Many phenomena in situations out of equilibrium can be described using a picture based on reaction processes. Apart from chemical reactions, the examples are exciton quenching, recombination of charge carriers or radiation defects in solids, predator-pray relationships in ecology etc. Reactions in homogeneous media are often described by formal kinetic schemes. Thus, the concentrations C i (t) of the components follow the first-order differential equations dC i (t)/dt = f i {C 1 (t), ..., C N (t)} where the reaction terms typically have a form with the powers n j depending on the stoichiometry of the reaction and κ i denoting the corresponding reaction rates. In inhomogeneous situations (layered systems, fronts, etc.) the mesoscopic approach based on reaction-diffusion equations for the positiondependent concentrations C i (r, t) is often the appropriate way of description. In case of normal diffusion such equations are obtained by adding a diffusive term to classical reaction schemes and have the form with f i = f i {C 1 (r, t), ..., C N (r, t)} and K i being the diffusivity of the component i. This approach is applicable whenever characteristic scales of spatial inhomogeneities are much larger than the typical interparticle distances and particles' mean free paths (see e.g. [1]). As we proceed to show, the possibility to put down such schemes is due to the Markovian nature of normal diffusion. Nowadays more and more attention is paid to situations when the diffusion is anomalous, which are found to be abundant [2]. One of the most important situations here is the case of subdiffusion described within the continuous-time random walk (CTRW) scheme [3]. In this case the subdiffusive nature of motion stems from the fact that particles get trapped and have to wait for a time t (distributed according to power-law probability density function ψ(t) ∝ t −1−α ) until the next step can be performed. It was shown that the properties of the reaction under such subdiffusion might be vastly different from those in diffusive systems [4,5,6,7]. The microscopic approach of these works aims on the understanding of the situation when the particles performing CTRW react on encounter (and don't react as long as they do not move). Such situation is pertinent to exciton quenching in solids, or to transport in ion channels [8].
In many cases, however, a mesoscopic approach is desirable. Such an approach was adopted in the case of reactions under superdiffusion due to Lévy flights, Refs. [9,10], where the transport process involved is Markovian. The situation with subdiffusion is much more subtle due to strongly nonmarkovian character of subdiffusive transport [11]. Here two different situations can be encountered: either the reaction at small scales is also subdiffusion-controlled (like in the models discussed above, where particles can only react if a new step is made) or it locally follows normal, classical kinetics. This last case that we address here is physically relevant since it describes reactions in porous media. The situation is of extreme importance in hydrology, where the transport in catchments is hindered by trapping in stagnant regions of the flow, caves and pores on all scales. The transport at long times and large scales is adequately described by CTRW [12]. However on small scales reactions take place in normal aqueous solutions, so that particles trapped in stagnant regions still can react with each other. A mesoscopic approach to such a case was adopted in [13] within a probabilistic scheme, while [14] tackle this problem by using equations of the same form as our Eq.(1) where the diffusion operator is changed for a subdiffusion one, containing an additional fractional derivative in time [3,15]: In this equation α i is the exponent of the anomalous diffusion for the component i, 0 D 1−αi t is the operator of fractional (Riemann-Liouville) derivative, and K i,αi is the corresponding anomalous diffusion coefficient. Such equations with decoupled transport and reaction term were postulated based on the analogy with Eq.(1) and look quite plausible. In some cases also a reaction term has to be modified by applying a fractional derivative as suggested by a microscopic model in [5].
In what follows we derive the reaction-subdiffusion equations for the simplest reaction scheme (monomolecular conversion A → B) corresponding e.g. to radioactive decay of isotope A which is introduced into the ground water at some place at time t = 0 and is transported according to anomalous diffusion. We show that the corresponding equations do not follow a pattern of Eq.(2), so that the reaction and diffusion terms do not decouple.
Let us assume for the time being that all properties of A and B particles are the same, so that the reaction corresponds to a relabeling of A into B taking place at a rate κ. In what follows we will use one-dimensional notation, the generalization to higher dimensions is trivial. The Eqs.(1) for this case read: with K being the normal diffusion constant.
It evolves according to a diffusion equation: Both concentrations A and B follow To see this apply Laplace-transform to the equations (3). Note that the solution for A(x, t) in the Laplace domain isÃ(x, u) =C(x, u + κ). Eq.(5) reflects the fact, that the conversion is independent from the motion of particles, so that concentrations of As and of Bs are proportional to the overall concentration multiplied by the probability for a particle to survive as A or to become B. The same argument leads to the conclusion that Eq.(5) also holds in anomalous diffusion, whatever the evolution equation for C is. For subdiffusion so that in Fourier-Laplace domain for the initial condition However, neither the solution of Eq.(2) nor the solution of the fractional equation with the modified reaction term [5] reproduce this result: the simple reaction-subdiffusion schemes do not describe the conversion reaction correctly.
Let us now turn to deriving a correct reaction-diffusion equations for our case. Our derivation will follow the way of derivation of the generalized master equation for CTRW used in ref. [16] based on the ideas of [17]. We start from a discrete scheme and consider particles occupying sites of a one-dimensional lattice. The generalized reaction (sub-)diffusion equations follows from the balance conditions for particle numbers. A balance condition for the mean number A i of particles A on site i of the system reads where I − i (t) is the loss per unit time due to the particles' departure from the site (loss flux) at site i, I + i (t) is the gain flux, and κA i is the loss due to conversion. Particles' conservation for transitions between the two neighboring sites corresponds to where w i,j is a probability to jump to site j when leaving i. For unbiased walks one has w i−1,i = w i+1,i = 1/2. Thus: We now combine this continuity equation with the equation for I − i (t) following from the assumption about the distribution of sojourn times. The loss flux at time t is connected to the gain flux at the site in the past: the particles which leave the site i at time t either were at i from the very beginning (and survived without being converted into B), or arrived at i at some later time t ′ < t (and survived). A probability density that a particle making a step at time t arrived at its present position at time t ′ is given by the waiting time distribution ψ(t − t ′ ), the survival probability being p(t) = e −κt . Thus: where ψ s (t) = ψ(t)e −κt is the non-proper waiting time density for the actually made new step provided the particle survived. This approach can also be generalized to bimolecular reactions [19]. Changing to the Laplace domain and noting thatψ s (u) =ψ(u + κ) we get withΦ κ (u) given bỹ Returning to the time-domain we thus get Note that Φ κ (t) given by the inverse Laplace transform ofΦ κ (u) corresponds to Φ κ (t) = Φ 0 (t)e −κt where Φ 0 (t) obtained by taking κ = 0 is the usual memory kernel of the generalized master equation for CTRW. Combining Eq.(15) with Eqs. (8) and (9) we get: Transition to a continuum in space (x = ai) gives (17) a rather unexpected form, where the reaction rate explicitly affects the transport term. For the exponential waiting time distribution ψ(t) = e −λt corresponding toψ(u) = λ/(u + λ) the kernel reads Φ 0 (t) = λδ(t), and the existence of an additional exponential multiplier does not play any role: The reaction diffusion equation is perfectly exact.
In the case of slowly decaying Φ 0 (t) the exponential cutoff introduced by the reaction is crucial. For power-law waiting time distributions and for κ = 0 the integral operator t 0 Φ 0 (t − t ′ )f (t ′ )dt ′ is the operator of the fractional derivative: For such distributions ψ(u) ≃ 1 − (τ u) α Γ(1 − α) (where τ is the appropriate time scale) and (for u → 0) we haveΦ 0 (u) ≃ (1/τ α Γ(1 − α))u 1−α which is proportional to the operator of the Riemann-Liouville derivative of the order α: For κ > 0 however the reaction affects the diffusion part of the equation: the Laplace transform of the integral kernel Φ κ (t) reads and is no more a fractional derivative. The integral op- (19) turning to be a fractional derivative only for κ = 0. The equation for the A-concentration thus finally reads: Although our reaction does not depend on the particles' motion, the parameters of the reaction explicitly enter the transport operatorT of the equation. Analogously we will now derive an equation for the Bparticles. As for A one has a balance condition for the mean particle number B i of B-particles on site i where J + i denotes the gain flux and J + i the loss flux of particles B at site i. The continuity equation reads:

t). (22)
A B-particle that leaves the site i at time t either has come there as a B-particle at some prior time or was converted from an A-particle that either was at site i from the very beginning or arrived there later, at t ′ > 0. Thus: where the local balance equation (8) was used. Applying Laplace transform and solving forJ − i (u) we get: Here the initial condition B i (0) = 0 was explicitly used. Using Eq.(13) forĨ − i we get: After inverse Laplace transformation we get: We now substitute this into the continuity equation (22), perform the transition to a continuum and get For the exponential waiting time density for which Φ 0 (t) = λδ(t) the third term in (27)  : Note that the equation for the product contains the term depending on the concentration of the component A at all previous times. This term has to do with the fact that the products are introduced into the system later on in course of the reaction, and their motion therefore is described not by the normal CTRW (and the corresponding fractional diffusion equation) but by the aged one [18]. Note also that the sum of Eqs.(20) and (28) always yelds the "normal" subdiffusion equation for the overall concentration C(x, t), Eq.(6). Moreover the solutions of Eqs. (20) and (28) satisfy Eq. (5).
In summary, we derived here the equations describing the time evolution for the local A and B-concentrations in a simple monomolecular conversion reaction A → B taking place at a constant rate κ and under subdiffusion conditions. These equations do not have a usual form of reaction-diffusion equations with the transport term independent on the reaction one. This fact is due to nonmarkovian property of the subdiffusion process, and will persist for more complex reaction schemes as well.
IMS gratefully acknowledges the hospitality of the University of Barcelona and the financial support by the HPC-Europa program. Authors are thankful to J. Klafter and R. Metzler for stimulating discussions.