Kinetics of subdiffusion-assisted reactions: non-Markovian stochastic Liouville equation approach

Anomalous specific features of the kinetics of subdiffusion-assisted bimolecular reactions (time-dependence, dependence on parameters of systems, etc) are analysed in detail with the use of the non-Markovian stochastic Liouville equation (SLE), which has been recently derived within the continuous-time random-walk (CTRW) approach. In the CTRW approach, subdiffusive motion of particles is modelled by jumps whose onset probability distribution function is of a long-tailed form. The non-Markovian SLE allows for rigorous describing of some peculiarities of these reactions; for example, very slow long-time behaviour of the kinetics, non-analytical dependence of the reaction rate on the reactivity of particles, strong manifestation of fluctuation kinetics showing itself in very slowly decreasing behaviour of the kinetics at very long times, etc.


Introduction
The kinetics of some physical and chemical diffusion-assisted processes in disordered media are known to exhibit anomalous specific features [1,2], which cannot be described by conventional kinetic models [3]. Specific examples of such processes (excitation quenching, electron-hole recombination, exciton annihilation, etc) as well as the results of a number of experimental and theoretical investigations are summarized and thoroughly discussed in a number of review articles [1,4,5].
These kinetic peculiarities, which are believed to be due to memory effects on migration of reacting particles, have been investigated thoroughly within the anomalous variant of the continuous-time random-walk (CTRW) approach [1,2]. The approach is based on the model of particle migration as a set of stochastically independent jumps whose kinetics is controlled by the distribution W(t) of times of waiting for a jump [2,6,7] (it is also often denoted as ψ(t) [2]). The anomalous variant of the CTRW approach assumes the long-tailed form of W(t): W(t) ∼ 1/t 1+α (0 < α < 1), modelling the above-mentioned long memory (non-Markovian) effects typical for disordered systems [5]. These effects are known to reveal themselves in anomalous diffusion (subdiffusion) for which the mean square of displacement, r 2 (t) , grows slower than that known for the conventional diffusion: r 2 (t) ∼ t α (α < 1). Recently, the CTRW approach has been applied to the analysis of the kinetics of geminate subdiffusion-assisted reactions [8]. Unfortunately, the results of this work are represented in a rather cumbersome form which can hardly be suitable for applications. In general, the significant drawback of all the above-mentioned CTRW-based studies consists in oversimplified treatment of reactive interaction of particles at short (contact) distances, i.e., interference of non-Markovian migration and reactivity.
Similar drawback is typical for another popular approach based on the fractional diffusion equation (FDE) [9]- [12] which can, actually, be derived within the CTRW approach [5]. The important advantage of the FDE approach (FDEA) consists in the convenient and transparent form, in which it allows for representing final results, by analogy with the conventional Smoluchowski equation. However, the non-Markovian effects, manifesting themselves in the memory term in the FDE, lead to serious difficulties in FDEA-based treatment of the abovementioned interference of non-Markovian migration with Markovian processes (similar to those observed in the CTRW approach), for example, reactions, classical or quantum evolution, etc. Often, the corresponding FDE for such processes are introduced semi-empirically [5,13] and are not always accurate enough.
A rigorous and efficient method of analysing the memory effects in subdiffusion-assisted reactions has been recently proposed in [14]. The method is also based on the CTRW approach but, unlike other methods discussed above, essentially applies the Markovian representation of CTRW processes [14,15] which allows one to obtain rigorous expressions for the evolution operator of the system correctly describing the aforementioned interference of non-Markovian migration (modelled by the CTRW approach) with Markovian processes. This operator can be considered as the resolvent form of non-Markovian generalization of the conventional (Markovian) stochastic Liouville equation (SLE) [16] for the probability distribution function (PDF) of the system.
In this work, the non-Markovian SLE is applied to the detailed analysis of memory effects in bimolecular subdiffusion-assisted reactions. The anomalous properties of diffusion are found to manifest themselves in a considerable slowing down of the kinetics of both geminate and bulk reactions. In addition, subdiffusive character of relative motion proves to result in nonanalytical dependence of the rate of kinetically controlled reactions on reactivity, which cannot be reproduced by the semi-empirical variants of FDE [5]. Simple expressions are derived for long-time asymptotic behaviour of the subdiffusion-assisted reaction kinetics. The subdiffusive relative motion of particles is also found to show itself in anomalously slow fluctuation kinetics [17] observed at very long times.

General formulation
In this work we analyse the kinetics of subdiffusion-assisted bimolecular physical and chemical processes. In the conventional pair approximation [3], the kinetics is described by the timedependent PDF ρ(t) for the pair of reacting particles satisfying the master equatioṅ in whichK(t) is the microscopic reaction rate and ρ i is the initial PDF. For simplicity in applications, we will restrict ourselves to the one-channel processes for which K(t) is the just parameter (rather than operator) fluctuating in time due to stochastic relative motion of reacting particles.
From a mathematical point of view, equation (2.1) is similar to the Liouville equation for the PDF of a classical system governed by the 'Liouville' operator, −K(t); therefore, the problem under study is close to that of noise-induced relaxation in the dynamical systems. This similarity enables us to apply some methods of treating the relaxation [18], for example, using a recently developed method based on the CTRW approach [14].
In accordance with the above general formulation, fluctuations ofK(t) are assumed to result from stochastic relative migration of particles, say a and b, in pairs over the states |r a r b ≡ |r a |r b , corresponding to the location of these particles at the points with coordinates r = r a and r = r b , respectively, in the isotropic (discrete or continuum) six-dimensional (6D) space {r a ⊗ r b } with different ratesK =K r a ,r b . In reality, however, in absolute majority of processesK r a ,r b depends only on relative distance r between particles. For this reason it is worth making all further consideration with the coordinates In these coordinates the set of different ratesK =K r a ,r b ≡K r can conveniently be combined into the matrix:K = r |r K r r|, (2.3) in which |r are the states in {r}-space.
Hereafter, the states in {r}-space and corresponding matrices will be denoted with the use of bra and ket vectors: |r and r|, for states in {r}-space.
The reaction kinetics is determined by the evolution operatorR(t) averaged overK(t)fluctuations and the distribution of pairs over volume V (over R-coordinates): (

2.4)
This operator is equivalently represented in terms of the conditional evolution operator G(r, R; r i , R i |t)R Our further analysis shows that the time evolution of the systems under study are conveniently represented in terms of the Laplace transforms in time. We will often use this transformation conventionally denoted as Z( ) = ∞ 0 dt Z(t)e − t for any function Z(t). It is instructive to start our discussion with simpler Markovian models in which the problem reduces to solving the conventional Markovian SLE.

Markovian approach
In the Markovian approach, relative migration of particles in {r ⊗ R}-space leading to K(t) fluctuations is described by the PDF P(r, R, t|r i , R i , t i = 0) for a reacting pair satisfying the equation [19]Ṗ = −LP with P(r, R, 0|r i , R i , 0) = δ rr i δ RR i , (3.1) in whichL is some linear operator in {r ⊗ R}-space. In the Markovian approach (3.1), thê G(r, t|r i , 0) obeys the SLE [16]: withĜ(r, R, 0|r i , R i , 0) = δ rr i δ RR i , so that the Laplace transformˆ G is determined by the equation (in operator form) whereδ rR is the Kronecker delta symbol in {r ⊗ R}-space and therefore: The idea of the SLE (3.2) is based on the simple observation that for Markovian dynamical motion (equation (2.1)) and stochastic evolution (equation (3.1)) the total change Ĝ during short time t consists of two contributions: dĜ = −(KĜ) t and fĜ = −(LĜ) t, so that one gets the relation Ĝ = fĜ + dĜ = −(L +K)Ĝ t, which is equivalent to equation (3.2).
Noteworthy is also the fact that in the Markovian model (3.1), the operatorsL can be of very different analytical forms describing the general Kolmogorov-Feller jump processes [19], diffusion processes [19], the Lévy processes [5], etc.

Non-Markovian processes
In this section we will consider CTRW-based models of non-Markovian particle migration, for which the above-mentioned simple arguments do not hold, and which, nevertheless, allow for a relatively simple analysis of memory effects in the reaction kinetics.
4.1. Single-particle processes 4.1.1. CTRW approach. The CTRW approach treats particle migration as a sequence of sudden jumps in 3D {r}-space [2,6,7]. The onset of any particular jump of number j is described by the probability P j−1 (in {x}-space) not to have any change during time t and its derivative W j−1 (t) = −dP j−1 (t)/dt, i.e., the PDF of times of waiting for the jumps. The functions P l and W l (t) are independent of l for l 1 whereas for l = 0 depend on the problem considered: For example, for non-stationary (n) fluctuations [2,6,7] W i (t) = W n (t) ≡ W(t).

Markovian representation of CTRW approach.
Here, we will outline the Markovian representation of the CTRW approach [14], which is very useful for the analysis of the problems under study. Suppose that the kinetics of (r → r )-jumps of the particle is controlled by the Markovian process in {q}-space (different from {r}-space) governed by the operatorˆ . The corresponding PDF σ(q, t) satisfies the equatioṅ describing evolution in {q}-space and equilibration if the operatorˆ has the equilibrium state |e q (ˆ |e q = 0) with the equilibrium population probabilities p e q : The q-process is assumed to control migration in {r}-space as follows: (r → r )-jumps occur with the rate κ r r whenever the system visits the transition state |t in {q}-space. Transitions can lead to the change in |q -state, i.e., to |t → |n -transition with |n = |t . For simplicity, we assume that |n = |t as well as thatˆ and |t -state are independent of the state in {r}-space.
The evolution of the system in {r ⊗ q}-space is described by the PDF matrix |ρ obeying the conventional SLE whereK µ (µ = d, o) are the transition matrices operating in {x ⊗ q}-space: with κ rr = r ( =r) κ r r andL = 1 −κ o /κ d is the operator which describes jump-like motion in {r}-space. The transition matricesK d ,K o , andL are assumed to correspond to isotropic migration in {r}-space. Equation (4.7) should be solved with the initial condition |ρ t=0 = |i q |q q|, where |i = q p i q |q and e q |i = q p i q = 1.

The function of interest is the Laplace-transformed PDF in {r}-spacê
It is determined by the Laplace transform |ˆ ρ( ) satisfying equation Hereafter, for brevity, we will sometimes omit the argumentˆ of the Laplace transforms of functions under study if this does not result in confusions.
The solution of equation (4.11) can be represented as a CTRW-like expression, which in what follows (by analogy with the Markovian model (3.1)) is considered as a solution of the so-called non-Markovian SLE [14].

Non-Markovian SLE.
The above-mentioned non-Markovian SLE is written in a quite compact resolvent matrix form [14]: (4.14) In equations (4.12) and (4.13),P is the matrix of jump probabilities with zeroth diagonal elements,L is the above-mentioned Kolmogorov-Feller-type operator describing isotropic jumplike motion in {r}-space, are the PDF-matrices of waiting times in whichĜ µ = t|Ĝ|µ , and (ˆ ) is the matrix ). Some important expressions forĜ µ (µ = n, i) andˆ can be found in [14]. Formulae (4.15) offer the representation ofŴ i (t) andŴ n (t) =Ŵ(t) in terms of the Green's functionĜ for the controlling {q}-process:Ŵ i (t) is the PDF matrix for the time of the first visit of |t -state, whileŴ n (t) =Ŵ(t) is that for the time of revisits of this state.
Notice that the representation (4.12) and (4.13) for the non-Markovian SLE is somewhat different from that proposed earlier [14]. However, both representations are equivalent and can equally be used for the analysis of the problem under study.
According to equations (4.12)-(4.15) the initial state |i (in {q j }) manifests itself only in the PDF matrixŴ i (t). In particular [2,6,7], in the non-stationary n-CTRW approach most interesting for our further consideration |i = |n , so thatˆ W i =ˆ W n ,ˆ P i =ˆ P n and Equation (4.16) can be rewritten as the non-Markovian SLE forˆ G: andδ rr is the Kronecker delta symbol in r-space. Comparison with equation (3.3) shows that the non-Markovian nature of migration manifests itself in the occurrence of the memory operatorˆ M in the SLE (in the right-hand side of equation (4.17)), which significantly affects the migration kinetics. In the Markovian caseˆ (ˆ ) ∼ˆ and, therefore, equation (4.17) reduces to the corresponding Markovian SLE (3.3).

Anomalous memory.
In the CTRW approach the model of the process is specified by the operatorˆ ( ). In the simplest model of anomalous migration [5] ( are the matrices of jump (fluctuation correlation) rates and anomalous powers, respectively, with 0 < α r < 1. Both matrices are diagonal in |r -basis. The model (4.18) describes anomalously slow decay of the matrixŴ(t) ∼ 1/t 1+α (very long memory effects in the system [5]), for which only the case of non-stationary (n) fluctuations is of physical sense (see equation (4.3)).
In what follows, for simplicity, we will assume that w r and α r are independent of r, i.e.,ŵ andα are just parameters: In this simple variant of the model (4.18) where ∇ r is the m-dimensional gradient operator in r-space, λ 2 is the average square of the jump length, assumed to be independent of r and u r is the external interaction potential. In our work we will mainly concentrate on the analysis of the free diffusion model, i.e., assume u r = 0, though some manifestations of the potential will also be discussed.
In general, the approximation (4.22) leads to the non-Markovian Smoluchowski-type equation [14,15] forĜ which for the Laplace transform is represented as In general, for multichannel processesK r will be a matrix; however, as we have already mentioned above, for simplicity we will consider only the one-channel problems in which K r are just functions.
Noteworthy is the fact that equation (4.17) can be treated as the Laplace transform representation (4.17) of the conservation law for number of particles:Ġ = −∇ r · J, where J is the corresponding flux whose Laplace transform in time J is given by (4.25) It is also worth noting that equations (4.23) and (4.25), the product of Laplace transforms M G with M given by equation (4.18), should be treated as [5] where 0 D 1−α t is the Riemann-Liouville fractional derivative operator.

Two-particle processes
4.2.1. General remarks. As we have mentioned in section 2, the diffusion-assisted reactions are quite accurately described within the pair approximation in which the reaction kinetics is determined by the space/time evolution of the two-particle PDF. Evidently, if one of the reacting particles is immobile the problem reduces to the singleparticle one discussed in section 4.1. In majority of migration-assisted reaction processes, however, both reacting particles are mobile. If the migration is Markovian and the PDF satisfies equation (3.1), the generalization to the two-particle processes is straightforward. To clarify this generalization we suppose, for definiteness, that the particles undergo isotropic diffusive motion in the potential u r , which depends only on the interparticle distance r = |r a − r b |. In such a case, operators L a and L b , governing spatial evolution of particles a and b respectively, are given by equation (4.22) with different parameters λ = λ a and λ = λ b for the two particles.
To describe time evolution of the two-particle PDF in the Markovian case it is sufficient to replace the single-particle migration operator L by the two-particle operator L ≡ L ab = L a + L b describing migration of the particles a and b. This means, in particular, that in the absence of interparticle interactions (u r = 0) the two-particle evolution operator is represented as a product of single-particle operators:Ĝ(t) =Ĝ a (t)Ĝ b (t), whereĜ µ (t) = exp(−L µ t). If interaction exists but depends only on interparticle coordinate r = r a − r b , same as the initial PDF ρ i : ρ i (R, r) ≡ ρ i (r), the problem essentially reduces to the analysis of evolution in r-space governed by L r = L a r + L b r .
For non-Markovian migration, however, situation is more complicated since the time dependence of the evolution operatorsĜ µ (t) is strongly non-exponential. This means that, in the presence of memory effects, the space/time evolution of the two-particle PDF is determined by both operators L a and L b , rather than by the sum L a + L b only. This can be easily demonstrated for the case of two non-interacting particles in which the two-particle evolution operator is still represented asĜ(t) =Ĝ a (t)Ĝ b (t). For ρ i (R, r) ≡ ρ i (r) depending only on r the averaging ofĜ(t) over the pair position R (see equation (2.2)) gives the following relation for averaged evolution operators Ĝ µ (t) R ≡Ĝ µ r (t): (4.27) Equation (4.27) is valid not only for non-interacting (freely diffusing) particles, but also in the presence of interparticle interaction, u r . This can easily be proved in the diffusion approximation for L a and L b (see equation (4.22)) with the use of the Markovian representation [14] outlined in section 4.1. It is sufficient to assume that relative migration of particles a and b in {r}-space (r = r a − r b ) is governed by two statistically independent Markovian processes in q a -and q b -spaces (as described in section 3) which lead to the jumps in r-space described by operators L a and L b (see equation (4.8)). The only thing one should take into account is that after averaging over R all terms containing derivatives ∇ R vanish. All other terms are combined into relation (4.27) in whichĜ µ −r =Ĝ µ r (because the operators L µ r are of the second order in the derivative ∇ r and the potential u r depends only on r = |r|).

Yield of geminate reaction. Equation (4.27
) allows one to make some general conclusions on steady-state properties of non-Markovian migration in the presence of a reaction. Here, we will calculate the yield Y 0 r = Y r (t → ∞) of the geminate reaction of ab-pairs in the limit of high reactivity which will be modelled by the infinitely deep box-type potential well u r = −u 0 θ(d − r), where u 0 → ∞.
We assume that ab-pairs are initially created isotropically at a distance r = r i > d. In this case, the geminate reaction yield Y r can conveniently be evaluated in terms of the yield Y e of escaping from the reaction It is important to note that, since the operators L a and L b commute with each other (differing only in the average square of the jump length λ 2 ), they have the same set of eigenvectors ψ k (r). For r > d, at which L µ = λ 2 µ ∇ 2 r , they are written as where In this equation we have used the relation G µ k=0 (t) = 1 (µ = a, b), which can easily be seen from equation (4.16) taking into account that in the Fourier transform G µ k→0 (t) is still given by this equation but with L µ replaced by the term ∼k 2 → 0. Expression (4.33) shows that under fairly general assumptions about the migration mechanism the conventional formula for the reaction yield Y 0 r = 1 − Y 0 e = d/r i is valid independently of the mechanism.

Anomalous diffusion.
The general formulas obtained above enable one to find some simple expression for the pair evolution operatorĜ =Ĝ ab in the case of anomalous diffusion (subdiffusion) of particles.
It is important to note that, in the case of free diffusion, in the region r λ µ t α/2 with high accuracy we can use the approximation It is clear that for long times the region r λ µ t α/2 covers all distances which determine the reaction process. This means that the approximate representation can be applied to the description of long time reaction kinetics. It is also clear that these representations predict the same asymptotic time dependence of the mean square displacement: r 2 (t) ∼ t α . In addition, in accordance with the above-mentioned general results, both representations reproduce the conventional expression for the reaction and escaping yields.
In other words, in the case of subdiffusive motion of particles, the expression (4.35) similar to the single-particle case (analogously to the case of conventional diffusion) seems to be fairly accurate at long times and can be applied to the analysis reaction kinetics in the two-particle approximation.
In our further discussion we will use the expression (4.35) for the description of both the cases of one and two mobile particles, and the situations when its accuracy is not very high will also be emphasized.

Applications and discussion
Here, we will consider some applications of general expressions obtained above. For simplicity, we will restrict ourselves to the analysis of isotropic 3D migration in the spherically symmetric potential u r . The migrating particles are assumed to be created with isotropic initial spatial distribution. In this analysis we will consider the simplest case of spherically symmetric onechannel reaction in which the microscopic rate K r .
The evolution of the system in such isotropic processes is expressed in terms of the spherically symmetric evolution operator (in the Laplace transform form) where r = |r|, It is seen that the problem of describing anomalous reaction kinetics reduces to solving the Schrödinger-type equation (5.2). In general, this can be done only numerically; however, we will mainly concentrate on the analytical analysis of the long time asymptotic behaviour of reaction kinetics, which can demonstrate the most interesting specific features of the kinetics. We will start this analysis with some steadystate properties of the processes under study.

Reaction yield in the case of localized reactivity
Very pronouncedly anomalous effects on subdiffusion-assisted reaction manifest themselves in the yield Y 0 r of reactions of particles with finite short-range reactivity. The specific features of these effects are conveniently illustrated by the case of free diffusion (for u r = 0). For simplicity, we will suggest the memory function M r (or r to be independent of r (see equation (4.20)).
In the considered case of spherically symmetric diffusion and isotropic initial condition, the corresponding SLE reduces to the 1D equation σ(r| ) = r ρ(r| ). (5.4) Suppose that the particles are created at r = r i : We also assume that the reactivity is highly localized at contact distance d: In this case, the 1D SLE for σ(r|t) is written as [(κ/w) α + λ 2 d 2 /dr 2 ]σ = 0, for d < r < d + l. (5.8) Solution of these equations by well-known methods give with Q r = (ld/λ 2 )(κ/w) α . In the limit of high reactivity, Q r 1, the process is diffusion controlled and Y 0 r = d/r i as should be. Equation (5.9) is similar to that known in the theory of processes controlled by conventional diffusion. However, the dependence Y 0 r (κ) determined by the function Q r (κ) ∼ (κ/w) α differs from the conventional one (corresponding to case α = 1) and cannot be reproduced by the semi-empirical FDE approach [13]. The dependence Q r (κ) ∼ (κ/w) α can easily be understood as follows: the subdiffusive system spends in the small reaction region l the time τ s ∼ w −1 (l/λ) 2/α , so that the effect of reactivity κ, expected to be strong when the reaction rate satisfies the relation κτ s ∼ 1, is similar to that predicted for the conventional diffusion with the term κ/w replaced by (κ/w) α , in agreement with equation (5.9). Since for anomalous diffusion the time τ s is much larger than that for the conventional diffusion, the obtained relation means that in the processes assisted by anomalous diffusion the effect of reactivity is stronger than in conventional diffusion-assisted reactions.
Recall that equation (5.9) is quite rigorous in the case of one mobile particle in the pair, but if both particles are mobile and undergo subdiffusive motion, it should be considered as approximate, in accordance with the results of section 4.2. In the high reactivity limit Q r 1, however, it is valid independently of the mobility of the second particle. As for finite reactivities, the expression for Y 0 r is also expected to be represented in the form (5.9) with Q r ∼ (κ/w) α independent of the mobility because of the above simple qualitative arguments essentially based on the dependence r 2 ∼ t α , which is valid both for one and two mobile particles.

Long time behaviour of reaction kinetics
The general formulas derived above enable one to fairly easily describe the kinetics of anomalous diffusion-assisted reaction processes. For brevity and simplicity, in this work, we will analyse only long-time behaviour of this kinetics. We will also restrict ourselves to the case of freely diffusing particles suggesting u r = 0.

Kinetics of geminate reactions.
In accordance with the above general analysis we assume that the pair of reacting particles is created with spherically symmetric distribution at a distance r = r i , so that the initial distance PDF is ρ i (r) given by equation (5.5). As in the steady-state case the reactivity is described within the model of highly localized reactivity (5.6).
The reaction kinetics, i.e., the time-dependent reaction yield Y r (t), can equivalently be characterized by the reaction flux J r (t) = dY r /dt, whose Laplace transform J r ( ) is written as where g(r, r i | ) is defined in equation (5.1). In the considered case of non-interacting particles (u r = 0) the expression for g(r, r i | ) and, thus, for J r can easily be obtained analytically: The asymptotic time dependence of J r (t) at long times t w −1 (d/λ) 2/α is determined by the behaviour of J r ( ) at small for which Y 0 r = l r /r i is the total reaction yield and l r is the effective reaction radius defined in equation (5.9). The inverse Laplace transformation of formula (5.12) yields Naturally, for α = 1, equation (5.13) reduces to that obtained earlier for reactions assisted by conventional diffusion [22]. It is interesting to note that formula (5.12) can also be represented as (5.14) where is the Laplace transform of the concentration of one reacting particle at the point of location of the other particle in the absence of reaction (the evolution operator for the PDF of the particle in this case is denoted as G 0 (r, r i | )). Formula (5.15) generalizes the corresponding expression derived in [22] for conventional diffusion-assisted reactions. The important modification of this expression in the case of anomalous diffusion consists in the replacement of the conventional bimolecular reaction rate K B = 4πl r λ 2 w by the anomalous one (see equation (5.14)) describing the long time-tailed memory effects which show themselves in the term M( ).

Kinetics of bulk reactions.
Here, we will analyse some specific features of the kinetics of irreversible anomalous diffusion-assisted bimolecular reaction of particles a and b within the mean-field (Smoluchowski) approximation [3]. In this approximation the equations for concentrations C a (t) and C b (t) are written aṡ where k ab (t) is the bimolecular reaction rate whose Laplace transform k ab ( ) k ab = d 3 rK r ρ(r| ) = 4πd 2 M∇ r ρ(r| )| r→d (5.17) is expressed in terms of the Laplace transform of the reactive-pair spatial distribution function The rate k ab (t) is similar to the reactive flux J r (t) introduced above for describing the kinetics of geminate reaction. The function ρ(r| ) satisfies the equation which can easily be solved analytically. In the limit → 0, corresponding to the considered limit of long times, . (5.20) Substitution of this solution into equation (5.14) yields where K B ( ) is the Laplace-transformed effective anomalous bimolecular rate defined in equation (5.14). After inverse Laplace transformation, one gets In particular, if the concentration of particles b is much larger than that of a-particles, i.e., if C b C a then equations (5.10) and (5.22) predict the reaction kinetics is the initial concentration of b-particles. This formula generalizes the expression obtained earlier in the diffusion-controlled limit, when l r = d [1,13].
The brief analysis of the kinetics of bulk anomalous diffusion-assisted reactions presented above shows that it can be described by expressions similar to those known in the theory of reactions assisted by conventional diffusion in which, the bulk reaction rate K B = 4πl r D, however, should be replaced by the Laplace-transformed effective anomalous rate K B ( ) (see equations (5.14) and (5.21)).

Fluctuation limit of bulk reaction kinetics.
The mean field (Smoluchowski) approximation applied in section 5.2.2 to the analysis of the kinetics of bimolecular diffusionassisted reactions is known to be correct at relatively small times (which, however, in the case of conventional diffusion cover the whole region of practically interesting times). At very long times this kinetics is changed for the so-called fluctuation case [17].
Here, we will briefly discuss the specific features of the fluctuation kinetics of subdiffusionassisted reaction. For the sake of rigor of further analysis we also assume particles b to be immobile: λ b = 0. In this case, the anomalous equations for the PDF are exact.
For certainty, we will consider the reaction in which the concentration of one type of particles, say b, is much larger than a: C a C b . In this limit, the difference between mean field and asymptotic kinetics is expected to be most pronounced. In particular, in the case of conventional diffusion-assisted reactions, the mean-field exponential kinetics C a (t) ∼ exp [ − A m (wt)] is replaced by a somewhat slower one C a (t) ∼ exp [ − A f (wt) 3/5 ] [17], where A m and A f are some constant parameters.
The mean-field kinetics of the reaction under study has already been evaluated above (see equations (5.23)). This kinetics proves to be of stretched-exponent type.
As for the fluctuation kinetics, it can be obtained using simple arguments applied earlier for the analysis of the conventional diffusion-assisted reactions [17]. The idea is that the long-time fluctuation kinetics is determined by the reaction of particles a which survive in relatively large regions (holes), where particles b are absent. For a semiquantitative consideration, it is sufficient to suggest these holes to be of spherical shape.
Suppose that the particle a is initially located at the centre of the spherical hole of radius l which does not contain particles b. The effect of b-particles outside this hole is modelled by the highly reactive (absorptive) boundary condition for the particle a at r = l. In accordance with the theory proposed above, the survival kinetics for the particle a, i.e., the survival probability P a (t), is determined by lowest eigenvalue k 0 of the diffusion equation Simple calculation gives k 0 = π/ l so that P a (t) ∼ l 2 . The observed fluctuation reaction kinetics C a (t) ∼ C 0 a P a (t) is obtained by averaging the probability P a over distribution P h (l) of radii of spherical holes: . Comparison of equations (5.23) and (5.27) shows that, in the case of anomalous (subdiffusive) motion of particles, the difference between the mean-field and fluctuation kinetics is much stronger than in the case of conventional diffusion: in the mean-field regime, the kinetics is exponential: C a (t)/C 0 a ∼ exp [ − Z a (t)], where Z a (t) = φ b −1 (1 + α)(wt) α , while in the fluctuation regime it is of very slow inverse power type: C a (t)/C 0 a ∼ (l 3 r C 0 b ) 1/3 /Z a (t).

Concluding remarks
The theoretical analysis of bimolecular subdiffusion-assisted reaction kinetics presented in this paper shows that the anomalous specific features of relative subdiffusive migration strongly manifest themselves in the kinetics of both geminate and bulk reactions. Despite fairly extensive investigations of these processes [1,5], there are still some problems to be discussed. In our work the analysis is made with the use of the non-Markovian SLE [14] rigorously derived within the CTRW approach. In the case of fractional diffusion (subdiffusion), the rigorous non-Markovian SLE looks similar to that obtained in the FDEA [13], but with some important differences which result in a number of different predictions of both equations, for example, on the form of the dependence of the reaction radius l r on reactivity κ (see equation (5.9)).
Our study demonstrates that the most important specific features of the subdiffusion-assisted reaction kinetics can be described by simple anomalous kinetic equations similar to those known in the theory of diffusion-controlled reactions, but with the conventional flux and bimolecular reaction rate replaced by the corresponding parameters, containing the anomalously long-tailed memory terms. These anomalous equations predict that subdiffusive motion of reacting particles manifests itself in the substantial slowing down of the reaction kinetics of subdiffusion-assisted reactions both geminate and bulk. It is also shown that, in subdiffusion-assisted reactions, the long-time fluctuation limit of the reaction kinetics is much more pronounced than in reactions assisted by conventional diffusion.
It is worth noting that the results on two-particle processes obtained in section 4.2 can be applied to general analysis of the statistical properties of relative random motion of particles undergoing anomalous diffusion. In particular, with the use of the expression (4.34) for the evolution operator G(r, r i |t), one can easily calculate the moments of the distribution of interparticle distance. Such an analysis will be the subject of a separate publication.