- Email: [email protected]
Phenomenological model for reaction kinetics coupled to a relaxing environment
Chemical Physics ELSEVIER
Chemical Physics 220 (1997) 25-41
Phenomenological model for reaction kinetics coupled to a relaxing environment Yuri A. Berlin a, * ,1, A l e x a n d e r L. B u r i n b, Sighart F. F i s c h e r a a Technische Universiti~t Miinchen, Physik Department T-38, D-85748 Garching, Germany b Department of Physics, CiO' College of the CiO' Universi~ of New York, 138-th Street at ConventAvenue, 10031 New York, NY. USA
Received 17 March 1997
Abstract Reasoning from the correlation between the height of the reaction barrier and the energy depth of substates attributed to different configurations of the reactant surroundings, the phenomenological approach to the rate processes in complex media has been proposed. The approach allows the formulation of the model based exclusively on the energetic consideration of the rate process rather than on the introduction of a particular configurational coordinate. This enables us to make conclusions concerning the most general features of rate processes in complex systems. In particular, we have identified three kinetic phases typical for chemical reactions proceeding in the complex molecular environment, including polychromatic (inhomogeneous), dispersive, and homogeneous phases. The expressions for the temporal behavior of the survival probability in each of these phases have been derived and compared with the empirical stretched exponential law. The temperature and scaling properties of survival probability has also been discussed. © 1997 Elsevier Science B.V.
1. Introduction It has become universally accepted that the classical phenomenology of chemical kinetics breaks down in complex systems like glasses, polymers and biomolecular objects (see e.g. Refs. [ 1-10]). A convincing body of experimental evidences shows that contrary to the classical phenomenology, the reaction patterns for complex systems cannot be adequately described in terms of the familiar time-independent rate constants (for review, see Refs. [2,4-6]); instead, the decay curves can be reproduced only by introducing the specific reaction rates which were found to decrease with time, t, following the near power-law dependence. The latter feature of rate processes in complex molecular media gave the reaction kinetics the name dispersive. The physical origin of the kinetic dispersion has received much attention in recent decades. At the early stage of investigations the dispersive reaction patterns were thought to arise only due to statistical reasons [11,12]. In particular, it has been assumed that a chemical reaction is characterized by a spectrum of activation energies,
* Corresponding author. I On leave from N. N. Semenov Institute of Chemical Physics, Russian Academy of Sciences, ulitsa Kossygina 4, 117334 Moscow, Russian Federation. 0301-0104/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S0301-0104(97)00120-1
26
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
rather than by a single value of this kinetic parameter (the so-called polychromatic concept). Then the decay curves, say, for a first-order thermally activated process, is easy to obtain by averaging the exponential time dependence for the survival probability of an individual reactant over the spectrum of reaction barrier heights (see e.g. Refs. [13-15]). If, in addition, the tunneling mechanism of reaction is considered [16-18], one should take into account the distribution of reaction barriers versus not only heights but also widths (for details, see e.g. [19-21]). As a result, decay curves calculated for first-order reactions turn out to be non-exponential in time both for thermally activated and for tunneling processes. This can formally be interpreted as a manifestation of the temporal evolution of the specific reaction rate, which was shown to proceed in agreement with the dispersive kinetic law, in the cases of certain distribution functions [2,21-25]. The key point in all early formulations of the polychromatic concept is the assumption that the distribution of reaction barriers over their heights and/or widths remains unchanged throughout the whole process. As has been recognized slowly, in view of more recent experimental data, such an assumption is restricted in its validity. For instance, the detailed studies of small ligand binding to heme proteins [7-9] lead to the conclusion that the barrier height distribution does not evolve over the time range (10 ns to 1 ks) as long as the temperature is below 160 K, while above 160 K it shifts as a whole to larger values of barrier heights without markedly changing the shape. Many other experimental evidences supporting the temperature variations of the barrier height distribution for different rate processes can be found in the review [6]. The experimental findings mentioned above indicate that statistics alone cannot provide the exhausted explanation of the kinetic dispersion observed for chemical reactions in complex molecular systems. Therefore significant efforts were made from the theoretical side to gain a deeper insight into the factors responsible for the temporal evolution of the specific reaction rates. Fairly fast progress has been achieved for elementary processes controlled by transport of active species in glasses and solid matrices. A short list of examples, which illustrate a success of theory, includes geminate and homogenous recombination of charged species [26-28], their scavenging by various acceptors [29-31], reactions of radicals [2] and H atoms [32], as well as the excitation transfer [33-35]. Reactants involved in these processes approach each other by successive transitions between structural traps randomly located in the bulk of the system (the bond disorder) or distributed in energy depths E (the so-called site disorder). Since the time of jump exponentially increases with E and with the 'trap-to-trap' distance, the transition rates for different pairs of trapping sites vary within a very wide range, making the resulting diffusion motion highly dispersive [36]. In turn, the latter circumstance leads to the dispersion in the specific reaction rate as verified by results of computer simulations [26,33] and analytical theory [27,28,31,37]. This enables one to conclude that the inherent structural disorder of complex molecular systems is a principle factor determining the dispersive kinetic behavior of active species in the course of diffusion-controlled processes. However available experimental data provides reason enough to believe that the temporal evolution of the specific reaction rate is also typical for the kinetic regime of chemical reactions, where reactants approach each other faster than a chemical event occurs [7-9,38-41]. It has been suggested that in the kinetic regime the dispersion can result from the coupling of the barrier crossing to the stochastic dynamics of the reactant surroundings (see e.g. [40,42-44]; for review, see [45] and references therein). A general framework for treating kinetic consequences of the coupling is provided by a theory of irreversible random transitions [46], which has been invoked earlier for studying randomly affected unimolecular reactions [47], elementary steps of hopping transport [48], the description of hierarchically constrained dynamics of configurational rearrangement in complex systems [49], and other rate processes [50]. Nevertheless, it should be borne in mind that the application of the theory to the situation under discussion calls for a concrete model of the coupling. A conventional approach used for the formulation of the model required is based on the introduction of the one-dimensional (i-d) stochastic process along the reaction coordinate distinct from the distance between reactive species (see e.g. [40,42,43,49,51-54]. If the coordinate is attributed to a configurational variable of the reactant surroundings x [49,51-54], the 1-d motion can be described in terms of sequential transitions between trapping sites corresponding to local minima of potential energy profile along the x-axis. Such minima are
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
27
associated with configurational substates of the chemically inert complex medium [49,52-54]. Taking into account that the substates should have a spread in their energy, again the problem can be formulated in terms of diffusion in disordered systems. Moreover, since a specific reaction rate depends on x due to the dependence of the reaction barrier on the configuration of the reactant surroundings, while x evolves with time because of diffusion, the rate process turns out to be coupled to the dynamics of configurational rearrangements of the molecular environment. It should be mentioned that the model of the stochastic impact considered above implies the knowledge of the specific degree of freedom which controls the reactivity of species in each particular system and over each particular time scale. So far such information remains unavailable. This calls for a more 'flexible' model for studying the effect of configurational rearrangements on the rate processes in complex molecular environment. A demand for the model increases in the light of recent experimental observations [9,55-57], molecular dynamic simulations [58-61], and theoretical findings [49,62-64] which clearly demonstrate that rate processes in complex systems proceeding over a wide range of time scales are in fact controlled by several degrees of freedom rather than by a single configurational coordinate. In this paper we present a phenomenological approach to the rate processes in complex media based on the analogy between the dynamics of their structural rearrangements and the transport phenomena in site-disordered systems. Our approach presumes the correlation between the height of the reaction barrier and the energy depth of substates (traps) attributed to different configurations of the reactant surroundings. This assumption has recently been justified theoretically for a particular case of relaxation in the system of interacting dipoles [65]. Our approach allows the formulation of the model based exclusively on the energetic consideration of the rate process rather than on the introduction of a particular configurational coordinate. This advantage enables us to draw the most general conclusions concerning the kinetic features of rate processes in complex systems and to show that the temporal dispersion in specific reaction rates indeed can arise at the kinetic stage of chemical reaction before the relaxation of the reactant surroundings is completed. The paper is organized as follows. Section 2 contains a brief discussion of our approach to the description of rate processes in relaxing complex media. In this section we also formulate the model by treating the population of configurational substates in terms of the master equation which represents a discrete version of the Feller integral equation commonly used for solving various kinetic problems (see e.g. [66,67]). The solution of the master equation presented in Section 2.2 provides the expression for the Laplace transform of the observable, i.e. of the survival probability. Using the expression derived, we demonstrate in Section 2.3 that this quantity has distinct temporal behaviors within different time and temperature ranges. The explicit form of the behaviors are specified in Section 3 using the exponential and Gaussian distributions as examples. The results obtained are discussed in Section 4, where conclusions followed from our findings are also summarized.
2. Model 2.1. A p p r o a c h and mathematical f o r m u l a t i o n
Here we consider a situation where the ensemble of complex objects (e.g. proteins) is subjected to the excitation leading to the formation of reactants inside each object. This event is supposed to trigger the rearrangement of the chemically inert part of a complex object proceeding through a large number of intermediate configurational substates. Any one of these substates, say i, is modeled as the trap with the energy depth E i, so that the lifetime of the system with respect to the thermally activated release of the i-th substate is given by
ri = Z w °
~B ~ / ,
28
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
where Z is the number of substates distinct from those associated with the i-th trapping site, w 0 is the frequency factor, k B is the Boltzmann constant and T is the temperature. Note that in Eq. (1) the subscript i labels the substates rather than their energies. These latter may have the same value for several substates with different i due to their energetic degeneracy. Qualitatively it is apparent that thermally activated transitions between sites of distinct energies lead to the gradual increase of the population of traps with the larger energy depth. As results, the reactant surroundings relaxes towards more stable configurations. Reasoning from the theoretical results obtained in the particular case of interacting dipoles [65], let us assume now that the deeper is the trap, the higher is the reaction barrier s i. Then for reactions proceeding due to the thermal activation, the rate coefficient can be written as
ki = koexp( _ +
) = ko exp( _ k--~ rlEi ) '
(21
where k 0 is the preexponential factor and r/ is the numerical coefficient. In Eq. (2) the particular value r / = 1 corresponds to the situation where the barriers for the reaction and for configurational rearrangements of the molecular environment have the same heights. In what follows 7/ is assumed to be greater or equal to unity for the sake of definiteness. The opposite case (i.e. "O< 1) can be analyzed in a similar way if one puts T1 = T/~7. Thus, the present model suggests that the relaxation of the molecular environment has influence on the reaction. Such an effect arises since the specific reaction rate depends on the reaction barrier e according to the Arrhenius law (2), while e evolves with time t due to the increase in the population of the deeper substates. Our main concern is the description of reaction kinetics coupled to the relaxation mentioned above. In order to achieve this objective, it is useful to note that the quantity recovered in typical experiments is the probability P(t) of the system remaining in the reactive state at time t (the so-called survival probability). By definition, P(t) is given by
P(t) = ~,pi(t),
(3)
i
where Pi(t) is the probability of the i-th substate to be populated at time t. The latter quantity satisfies the equation
dpi (t) dt
kiPi(t) + ~"(wiiPj(t) -wiJPi(t)) (4) J and the initial condition p~(t = O) = pO. Here Wi~ stands for the conditional probability, per unit time, of the transition from the state j to state i, while the symbol )Z' denotes the summation over all j 4: i. Eq. (4) represents the discrete version of the Feller integral equation commonly used for solving various kinetic problems (see e.g. [66,67]). For uncorrelated transitions, one can assume that 1
1
wii = wi = Z'(i and wi i = w i - Z'ci ,
(5)
where 7 t is given by Eq. (1). For such an assumption, Eq. (4) can be rewritten as ! dPi(t) dt = - ( Z w i + ki)Pi(t) + y" wjpj(t).
(6)
J
The same equation follows from the general analysis of the uncorrelated hopping in the particular case where the waiting time distribution is the exponential function of time with the relaxation rate dependent on the random variable [67]. Once the solution of Eq. (6) is known, the calculation of the observable, i.e. P(t), for any initial probability density pO is straightforward and reduces to the proper summation in Eq. (3).
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
29
2.2. Solution of the kinetic equation and the Laplace transform of the survival probability The problem formulated in Section 2.1 can be solved by introducing the Laplace transform
(7)
p,( s) = fo+~Pi( t) exp( - st)dt , which enables us to write the population of the i-th substate in the following form [68] pO 1 fii(S)
s+Zw
+
i4-k i
E'wj~i(,).
s4-Zw i4-k i j
(8)
Then the Laplace transform of the survival probability can be written as pO 1
fi( s) = E f i i ( s) = E i
+ E
i s4-Zwi4-ki
E'wjfij( s) •
i S j + Z w i 4- k i
(9)
In order to derive the final expression for F'(s), all we need to do is to calculate the sum in the second term of Eq. (9). For Z >> 1, this can be done multiplying both sides of Eq. (8) by w i and by a subsequent summation. Such a procedure gives
E L, J
pA'-"w/Is)=
wi pO
i s + Zw~ + k~ 1 -- ~_~ wi
(10)
i s4-Zwi4-ki
The combination of Eq. (9) with Eq. (10) and the subsequent replacement of the summation over i with the summation over energy of substates allow the continuous representation of the final result for the Laplace transform of the survival probability, viz.
/~(S)
~,w°(O
1/fgl(S) 4- I/F0(S ) , S + - -
(11)
,/,~(s)
where ~,°(s), ~ ( s ) , ~ l ( s ) and U/to(S) are defined in terms of the initial population distribution G(E) and the density of substates g(E) as ~ O ( s ) = fo~ Z w ( E ) G ( E ) d E
s + Z w ( E ) 4- k ( E ) ' g(EldE o s+Zw(E)+k(E)
(12)
,/tg(s) = f ~
qtxl(s)f~
'
k(E)g(E)dE
o ~ + Zw(E)
+ k(E)
'
(13)
(14)
G( E ) d E
~'o(,) = fo
s + Zw(E) + k ( E ) "
(15)
In order to avoid misunderstanding, we note that in mathematical literature the distribution G(E) is referred to as the frequency function, the probability density for the initial population or as the differential distribution function (see e.g. Ref. [69]).
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
30
The above result can be simplified further for several special cases of physical interest. First we consider the situation of very rigid media that corresponds to the limit where w(E) -- O. In this case Eq. (11) reduces to the integral G( E ) d E /~(s) = f0 s + k'(E--) ' (16) so that the subsequent inverse Laplace transform [68] immediately yields the expression for the survival probability derived earlier in the framework of the polychromatic concept (see e.g. Refs. [47] and [51]) oc
P ( t ) = £ G(E)exp( - k ( E ) t ) d E .
(17)
Another interesting situation relevant to the proteinquake hypothesis proposed by Frauenfelder et al. [70] becomes evident if we introduce the probability ~o(E) of the initial population for a single substate with the energy depth E and rewrite the initial population distribution as G(E) = ¢ ( E ) g ( E ) . In the case discussed in the framework of the above hypothesis the excess of energy dissipated into the media after the formation of reactants is assumed to be so high that q~(E) is independent of E. Under this condition, leg'(s) ~g(s)
~,(s) ~'0(s) '
(18)
grl(S) = ~ k(E)G(E)dE o s + Zw(E) + k(E)
(19)
where
Therefore, from Eq. (11) we have gto(S ) 1- f ~ Zw(E)G(E)dE
1
t5(s)
s + qe'(s)
S---S
Jo s + zw( e) +
(20)
e)
Eqs. (11)-(15) as well as Eq. (20) allow the derivation of the expression for the survival probability in several time limits. For the sake of simplicity, here we shall restrict our consideration to the situation, where Eq. (20) is valid. The analysis of the more general expression (11) will be presented elsewhere [71]. 2.3. Asymptotic time behaviors of the survival probability
As can be shown by expanding Eqs. (15), (19) and (20) in a power series in l / s , the Laplace transform of the survival probability for s >> Zw o > k o can be approximated as fi( s) =
1-
k
,
(21)
where the symbol ( . . . ) signifies the averaging over G(E). For the rate coefficient given by Eq. (2), ( k ) exists for any energy depth distribution which is normalized for E ~>0. Since the inverse Laplace transform of Eq. (21) is straightforward [68], we conclude that for t << 1 / ( Z w o) and for Zw o > k o, the survival probability exhibits the temporal evolution described by the linear function of t e(t) = 1-(k)t. (22) It is worth noting that Eq. (22) coincides with the time behavior predicted by Eq. (17) for t < 1/k o. Such a coincidence is not surprising since Eq. (17) results from the first equality in Eq. (20) as long as s >> Zw o. Thus, according to our model, P(t) follows the polychromatic kinetic law if t << 1/(Zwo). This can be expected since at the very early stage of the chemical reaction the reactant surroundings has no time for significant configurational changes and appears to be 'kinetically' frozen in the initial substates.
Y.A. Berlin et al. / ChemicalPhysics 220 (1997)25-41
31
Another asymptotic behavior of P(t) is expected at very long times, where s << min{geo(O)/(1/k2), gel(0)/(1/k)}. Now Eq. (20) becomes 1
fi(s)
gel(0 ) , s + - -
(23)
geo(O)
so that the inverse Laplace transform yields [68] P ( t ) = exp(
ge0(0 )gel(0) t) .
(24)
The analysis of the borderline case where 1 / ( Z w o) < t < 1/min{geo(O)/(1/k2), gel(O)/(1/k)} can be performed by invoking a method proposed in Ref. [72]. Following this method, we introduce the energy-dependent parameter 1
e(E)=
dlnG(E)
"
(25)
dE Then one can verify that the integrand in Eq. (15) will have a maximum if E satisfies the equation
e) - -
+ Zw( e) + k( e) Zw(e) + nk(E)
-
(26)
Based on physical arguments, G(E) is usually assumed to decrease with E as the energy depth of configurational substates becomes equal to or higher than a certain value E. >~0. Therefore e ( E ) > 0 for E >~E . . Besides, in the most commonly encountered situations e ( E ) does not increase with the energy depth when E >/E.. As a consequence, Eq. (26) has a solution only if T is below a certain value Tcr
e( E. ) Zw o + Tlko kB
(27)
Zwo + k o
The latter circumstance allows the conclusion that Eq. (15) gives distinct dependencies ge0 versus s in the highand low-temperature limits. Hence the time behaviors of the survival probability within these two limits should also be distinct. If the temperature is high and it exceeds Tcr, the integrals (15) and (19) approximately equal ge0(0) and gel(0), respectively. Therefore at high temperatures our result for P(t) is identical to Eq. (24). Otherwise (i.e. for T < e ( E . ) / k s ) , the values of these integrals are mainly determined by the energy depths E -= Emax(S) adjacent to the position of the maximum prescribed for the integrand in Eq. (15) by the solution of Eq. (26). Following Ref. [72], ge(s) and gel(s) in this limit can be evaluated by approximating G(E) in the vicinity of Emax(S) as G ( E ) = G(Emax(s)) exp
(28)
kaTmax
with AE = E - Emax(S) and kBTm,x = e( Em,x( S)). Such an approximation together with the introduction of the new variable z = A E / ( k B T ) enable one to show (for details, see [721) that for 77= 1 and for temperatures below Tmax
a(Emax(s))kBZf + ~ e x p ( - T z / T m a x ) d z ge°(s) =
s
.1_~
]'+ exp(-z)
=
G(Emax(s))k.T rr s sin(TrT/Tm.x)
(29)
32
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
Here the second equality sign corresponds to the result of integration given in Ref. [73], and Emax(S) is roughly given by Emax(s)
kaT In
~
s
t
.
(30)
Since in the particular case r/= 1 the integral (19) can be rewritten as ap,(s)
k0 - - [ 1
--saP(s)],
k o + Zw o
(31)
the substitution of Eqs. (29) and (31) in Eq. (20) yields "B'kBZa( Emax( S) )
[
P( s) =
Zw ° ]. ko sin('rrT/Tm,x) + _ _ r r k J G ( Emax( S) ) s ko + Zwo ko + ZWo
(32)
The same ansatz is applicable to the situation where 77> 1 and the reaction barrier exceeds the energy depth of configurational substates. Taking into account that for Z >> 1 the preexponential factor k 0 is much less than Zw o, one may neglect k(E) in the denominator of Eq. (15). Besides, it is evident that s < k o + Zw o within the time range under consideration and therefore aPt(s) can be replaced by aPl(0). As a consequence, for 7/> 1 we obtain
"l'¢kBZa( Emax( S) ) P( s) = s[rrG( Emax( s) )ksT + aPl(O)sin(rrT/Tmax) ] .
(33)
Note that this expression reduces to Eq. (32) as 7/ tends to 1 and k 0 remains less than Zw o. Thus, the inverse Laplace transform of Eq. (33) together with Eqs. (17), (22) and (24) allows the description of the asymptotic behavior of the survival probability within three characteristic time ranges 0 < t << 1/(Zwo), 1/(Zw o) < t < 1/min{aP(0)/( 1/k 2 ), aPl(0)/( 1/k)}, and t >> 1/min{aP(0)/(1/k 2 ), aPl(0)/( 1/k)}, respectively, if the initial population distribution of configurational substates G(E) is known. The explicit form of P(t) for these time ranges will be specified in the next section using the exponential and Gaussian distributions as examples.
3. Examples 3.1. Exponential energy distribution of substates In this case one should put G(E)=
{l
"~0exp -
1
= kRToexp -
0
forE>/0,
(34)
for E < 0,
where E 0 is the parameter of the distribution expressed here in terms of the effective temperature T0. Then the substitution of Eq. (34) into Eq. (17) demonstrates (see also Ref. [73]) that in the short-time limit (t << 1/Zw o) the survival probability can be expressed in terms of the incomplete gamma function T(a,x) and varies with t non-exponentially, namely
P ( t ) = --~o(kot)
~o3'
~7---~o,kot •
(35)
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
33
Using Eq. (21) and properties of the incomplete gamma function [73], one can verify that for k 0 < Zw o and for the very early stage of the reaction, where t << 1/(Zwo), Eq. (35) reduces to Eq. (22) with the mean rate coefficient given by T ( k ) = k0 r
(36)
+ tiT0
In order to establish the temporal behavior of the survival probability at t > 1/Zw o in the high-temperature limit T > T0, we need now to calculate the ratio ~ 1 ( 0 ) / ~ ( 0 ) . For k 0 < Zw o, this ratio coincides with the ratio ~ l ( 0 ) / q t ( 0 ) as long as T exceeds T0. This implies that except for very short times t << 1 / Z w o, the reaction kinetics remains exponential and follows the familiar exponential law P ( t ) = exp( - ( k ) t ) .
(37)
However the results for T < To are different. The substitution of Eq. (34) into Eqs. (15), (19) and (20) demonstrates that at low temperatures for s >> Zw o ,
qto(s) fi( S) =
(38)
S - I ( oIS) fl
for s < Zw o ,
1+
with parameters a and /3 being defined as
ol = Zwo
koTo sin('rrT/To)
1 + ( 7 7 - 1)
,
T /3 = - - < 1.
(39) (40)
ro
The inverse Laplace transform of Eq. (38) (see Ref. [68]) makes evident that in the low-temperature limit ( T < T0) the temporal behavior of the survival probability is non-exponential. For t << 1 / Z w o, P(t) can be described by Eq. (35) as in the high-temperature limit, while for t > 1/Zw o this quantity becomes equal to the Mittag-Leffler function E t 3 ( - ( t / a ) t J ) . By definition [74], the latter can be written as
( - 1 ) k ( t / a ) ~k E t 3 ( - ( t / c t ) t3) = • k=0
F(1 +ilk)
,
(41)
where F ( x ) is the gamma function. Since in the long-time limit E t ~ ( - ( t / o t ) ~ ) has the asymptotic expansion [74] =
-
E n=0
+o
we conclude that
1 r(l +/3) P(t) = I (t/ot)_t3 r(1 - / 3 )
for t<< c~, (42) for t>> a .
34
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
It is interesting to compare the time behavior of the survival probability following from Eq. (41) with the empirical stretched exponential law
commonly used for the description of the reaction patterns in terms of two fitting parameters o/1 and ¢tl (for review, see e.g. Refs. [2,4], and [6]). Numerical data presented in Fig. 1 allows the comparison between Eq. (41) and Eq. (43) for different values of the dispersion parameter /3. As can be seen, in the case where a 1 = a [ F ( 1 +/3)] - l / ~ and ~ =/3 the temporal behaviors predicted by these equations coincide only at the initial part of the decay curves as long as t is less than a few tenths of a. This can be expected reasoning from the short-time asymptotic behavior of the survival probability given by the upper line of Eq. (42). If, however, oq ¢ ~ and /3t 4=/3, the stretched exponential provides a reasonable approximation of Eq. (41) in a wider time range. So far our consideration has been focused on the cases, where the temperature is either far above or much below T0. The reaction kinetics in the vicinity of the characteristic temperature To merits a closer examination. 1.0 A 0.8 ¸
0.6 ¸ 0.4 ¸ 0.2 ¸
0.0 1.0 ~ ~
0.8
B
•
0.2
~
0.0 1.0
"o0.61, t , " i ~ ,
, C
0.8 0.6
~
0.4 0.2 0.0 0.0
2.0
4.0
6.0
8.0
10.0
t/a Fig. 1. Temporal behavior of the survival probability P calculated for different values of the dispersion parameter/3. Curves shown in panels A, B, and C are plotted for/3 = 0.2, 0.4, and 0.6, respectively. Solid curves are dependencies given by the Mittag-Leffier function, see Eq. (41). Dashed curves correspond to the stretched exponential given by Eq. (43) with a I = a and /3~ = ft. D o t - d a s h e d curves are the fits to the Mittag-Leffler time behavior. The fitting parameters are a I = 6 a and /31 = 0.75/3 (panel A), oq = 1.9or and fll = 0.75/3 (panel B), oq = 1.3or and /3t = 0.83/3 (panel C).
35
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
The calculations carried out for t << 1/Zw o show that in the short-time limit the time evolution of the survival probability at T = TO can also be described by Eq. (35). For the very early stage of the reaction where t << l / ( Z w o ) , the incomplete gamma function can be approximated as y(1/T1, kot) = 71(kot)t/n[1 - kot/(~l + 1)] if k 0 < Zw o [73]. Hence, in this case Eq. (35) becomes kot - -
P(t) = 1
- 1 -(k)t,
1+77
(44)
in agreement with Eq. (22). For longer times t >1 l / ( Z w o) and for Zw o > k o, the evaluation of the integral (15) yields 1
(
~ ° ( s ) = oZ w - l n
1+
~Wo) s '
(45)
while ~l(s) may still be approximated by ~1(0) = ko/(~lZwo). As a consequence, from Eq. (20) one obtains fi(s)---~ln
l+
s
"
(46)
Since the inverse Laplace transform of Eq. (46) is straightforward [68], we arrive at the asymptotic expression for the survival probability ~/~ (1 - e x p ( - Z w o t ) ) . P ( t ) = ko---
(47)
Eq. (47) is valid within the time interval 1 / ( Z w o) <~t <~( 1 / Z w o ) e x p ( T / [ T - To[) as long as the temperature is close to To. In addition, it reproduces the temporal behavior of the survival probability given by the lower line of Eq. (42) which holds for t > a even in the case where I T - To] << To. 3.2. Gaussian distribution
Another energy distribution of configurational substates frequently used in the literature [19,26,28,37,38,51] is defined as 1 G(E)
(E-
e0) 2
4tr z
=
0
for E>~ 0, (48) for E < 0 .
This is the so-called Gaussian distribution dependent on two parameters 60 and tr which predetermine the position of the maximum and the spread of substates in their energy depths, respectively. The preexponential factor in Eq. (48) arises due to the normalization condition o~
fo G( E ) d E = 1
and contains the error function [73] 2
eft(x)
=
x
,rrl/2 fo exp(-~:Z)d~:"
In order to derive the explicit form of P ( t ) for the above distribution in different time and temperature
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
36
ranges, we introduce the dimensionless parameter a = c0/(2 tr) and express the width of the Gaussian tr in terms of the characteristic temperature To = cr/k B. Then the mean rate coefficient becomes
(k)=k0ex p
1-
whereas the ratio ~(0)/1//1(0)
~(0)
k°exp
can
,+e~+_ ,
(49)
be written as
t1
r
,
2AT
]]
[( '
1 + erf A 1 + ~-~
(50)
It should be mentioned that for temperatures exceeding a certain value rtTo
for a = 0 ,
nTo
(51)
for a > 0,
a
Eqs. (49) and (50) provide the identical result, namely
,,,,
~(0)
(
for,,=O,
kBT]
k0ex p
(52)
for A > 0 .
This implies that the exponential and Gaussian distributions give the same results qualitatively in the high-temperature limit: The survival probability exhibits the exponential temporal behavior with the exception of the very early stages of the process, where P(t) decays non-exponentially following either Eq. (17) or its particular form (22) as long as t << 1/(Zwo). In the opposite case (i.e. at T < TG) the survival probability exponentially decreases with time only in the very late stage of the reaction (t>> 1/min{qt(O)/(1/k2), ~l(O)/(l/k)}). According to our theoretical analysis (see Section 2.3), the rate coefficient characterizing the temporal behavior of P(t) in this stage is given by Eq. (50). The latter expression allow the further simplification for two particular positions of the Gaussian distribution at the E-axis. If the Gaussian peaks at E = 0, one should put A = 0 in Eq. (50). This leads to the following expressions for qtl(0)/~(0) in the low-temperature limit
qt,(O)
~(0)
--}- e x p , - - ~ ]
k°T
for r/= 1,
7"°2
. . . . 7 (__~)] exp - -~-g "rrl/2(l - r / ) T o 1 + e f t
(53) for'q > 0.
EA. Berlin et al. / Chemical Physics 220 (1997) 25-41
37
In the other situation which corresponds to the relatively narrow Gaussian distribution with A/+/> 1 and with the maximum centered at E = e0 4=0, we get
( ( 2a7"0 // 0exp/
1/fl ( 0 )
-/_
T2 )
2"rrl/2(1 -- ~)T0 exp
f o r , = 1,
(
T
T2
for ~7> 0.
The comparison of Eq. (53) and Eq. (54) shows that the shift of the Gaussian distribution towards higher E values results in the increase of the reaction activation energy in the late stage of chemical process if it proceeds at T < TC. In addition, Eq. (53) and Eq. (54) evidence that the reaction kinetics in the long-time limit does not obey the Arrhenius law. The latter circumstances can be explained by the attaining of the dynamic equilibrium between the population of configurational substates and their release in the course of the relaxation process. The similar explanation has been proposed earlier for the non-Arrhenius temperature dependence of the diffusion coefficient for the anomalous transport in site-disordered systems [75-77]. In order to complete the case under consideration, we now turn our attention to the intermediate time range 1/(Zw o) < t < 1/min{~(O)/(1/k2), qtl(O)/(1/k)}. As follows from the general results presented in Section 2.3, the temporal evolution of the survival probability within this time range can be deduced by taking the inverse Laplace transform of Eq. (33). For the Gaussian distribution (48), the use of Eqs. (25) and (30) allows the conclusion that parameters G( Emax(S)), Tmax, and qq(0) appearing in Eq. (33) can be approximated as exp( -A2 )
( s )~
G(Emax(S))= "rrl/2o'[1Terf(A)] ~ w0
(55)
'
½T Tmax
I
(56)
AT '
/3+s-7. °
qtt(0)=
k°
[
Zwo exp
2(r/-1)AT°( ~ 1
(+/-1)T°)] 2-AT
t,
(57)
1 + eft(A)
with /3 only slightly dependent on s, namely /3=
In
To
The substitution of Eqs. (55)-(58) into Eq. (33) and the subsequent inverse Laplace transform [68] allow the derivation of asymptotic expressions for P(t) similar in form to Eq. (42), but with the characteristic time c~ and the dispersion parameter/3 given by 1 [
"rrl/2Texp(-A 2)
o~= ZWo a/Zl(0)T0(1 +erf(A))sin(.rrT//Tmax ) ( T):
/3= -~o
ln(Zw°t)-z
AT
]1/3 .
(59) (60)
This implies that the dispersive reaction patterns can be expected both for the exponential and for the Gaussian energy distribution of configurational substates, as long as kinetic measurements are performed within the time
38
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
range l/(Zw o) < t < 1/min{q~(O)/(1/k2), x/rj(0)/(1/k)} and at temperatures lower than a critical value T = T *. The latter is related to parameters of distribution functions. In the case of the exponential distribution (34) T * = T0, while for the Gaussian (48) the critical temperature T * equals T~ and has the explicit form given by Eq. (51).
4. Discussion and conclusions
In the present paper we propose a phenomenological model for rate processes proceeding in the chemically inert environment relaxing towards the dynamic equilibrium. The equilibrium arises as a result of two competing processes, i.e. the repopulation of trapping sites associated with the configurational substates of the reactant surroundings and their release due to the thermal activation. The model suggests the correlation between the height of the reaction barrier and the energy depth of traps corresponding to the configurational substates: the deeper traps are, the higher is the barrier. Such an assumption was justified theoretically for the particular case of relaxation in the system of interacting dipoles [65]. The main physical consequence of the correlation mentioned above is the coupling of the chemical reaction to the relaxation process. Starting with the pioneering work of Agmon and Hopfield [42,51], the influence of the coupling on the reaction kinetics is theoretically treated in the framework of the approach based on the adequate selection of the specific configurational coordinate controlled the rate process in each particular system. Since the experimental information needed for making the proper choice is inaccessible in most cases, the application of the approach to the description of reaction kinetics in concrete complex systems is usually hard to justify. By contrast, a model presented here treats the chemical reaction coupled to the relaxation of the molecular surroundings in terms of the only variable, namely the energy of configurational substates E. This enables us to avoid difficulties related to the proper choice of the configurational coordinate and to draw the most general conclusions concerning the kinetic features of rate processes in complex systems. According to our theoretical findings, such features include (i) the existence of the polychromatic phase of rate processes [ 11-15] known in the biophysical literature as inhomogeneous kinetics (see e.g. Refs. [78] and [79]) arising as long as t << 1/(Zw o) or if physical conditions prevent conformational rearrangements of the surroundings (very low temperatures, rigid media, etc.); (ii) the transition from the polychromatic to the dispersive phase; the advent of the latter signals the onset of the relaxation process associated with configurational changes in the reactant surroundings; (iii) the crossover to the exponential kinetics in the long-time limit associated with the completion of relaxation; (iv) the non-Arrhenius temperature dependence of mean rate coefficient for the exponential phase of the rate process (cf. Eqs. (53) and (54)); (v) the existence of the critical temperature T * dependent on the parameters of the energy distribution of the configurational substates, serving as a demarcation line between exponential and dispersive kinetics in experiments performed in a wide temperature range and within the fixed time window of the appropriate position and width; (vi) the scaling properties of the survival probability in the two asymptotic time limits which are shown to exist at temperatures below T * prior to the completion of the relaxation process (see the upper and the lower lines of Eq. (42)). In addition, the proposed model enables us to derive expressions for the observable, viz. for the survival probability P(t), which describe its temporal evolution during each of kinetic phases mentioned above. In particular, we show that in the dispersive phase P(t) can be expressed in terms of the Mittag-Leffler function E~(-(t/a)~), c.f. Eq. (41). Our results shown in Fig. 1 clearly demonstrate that the empirical stretched exponential law (see Eq. (43)) commonly employed for fitting kinetic data in the case of the first-order reactions (for review, see e.g. Refs. [2,4,6]), exhibits practical the same temporal behavior as Et~(-(t/a)~) within the
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
39
restricted time range 1/(Zw o) < t <~ 10a. In order to attain the coincidence, one should assume, however, that parameters of the stretched exponential decay are distinct in their values from those used in numerical calculations of Et3(- (t/oz)13). This, of course, does not provide arguments for deciding between the stretched exponential and the Mittag-Leffier functions. Nonetheless, our numerical findings clearly show that the use of the stretched exponential, for the description of rate processes in complex systems, can be justified only by comparison with kinetic curves measured within a wide time window and reflecting the temporal evolution of the process up to high degrees of the chemical conversion ( P ~< 0.1). Otherwise, the fitting procedure invoking this function may lead to the values of dispersive parameter and the characteristic relaxation time being far from real. There are grounds to believe that chemical reactions in real complex systems indeed proceed in three phases expected from the proposed theory. This is especially true for processes in biological objects, in particular for ligand binding to heme proteins. In the latter case multiple pulse experiments [78-80], kinetic hole burning [7], and flash photolysis with monitoring in the near infra red [81] strongly support the existence of these kinetic phases. Data obtained by these methods provide additional information concerning phase durations at different temperatures. Finally, we would like to note that a theoretical analysis presented above is restricted to the case, where the excess of energy dissipated into the media after the formation of reactants is much higher than the mean energy depth of substates. This corresponds to the proteinquake concept proposed by Frauenfelder [70]. Nevertheless, the expression (12) for the Laplace transform of the survival probability allows a more general analysis without invoking special assumptions concerning the amount of dissipated energy. The results for the latter case which will be published elsewhere [71] show that the main conclusions followed from our study remain unchanged, but the temporal behavior of the reactant becomes more diverse in comparison with kinetic features deduced in the present work.
Acknowledgements This work is supported by the Deutsche Forschungsgemeinschaft. YuAB greatly acknowledges the DFG grant No.436 RUS 1 7 / 1 6 / 9 7 . ALB is indebted to the Alfred Toepfer Foundation and to the Netherlands Organization for Scientific Research. The authors would like to thank Dr. A. Zharikov for discussions and Prof. Dr. H. B~issler for reprints of his works.
References [1] Yu.A. Berlin, J.R. Miller, A. Plonka (Eds.), Rate Processes with Kinetic Parameters Distributed over Time and Space, (Chem. Phys. 212 (1996)). [2] A. Plonka, in Lecture Notes in Chemistry, Vol. 40, Springer, Berlin, 1986. [3] H. Frauenfelder, G.U. Nienhaus, J.B. Johnson, Ber. Bunsenges. Physik. Chem. 17 (1988) 451. [4] A. Plonka, Prog. React. Kinetics 16 (1991) 157. [5] A. Plonka, Annu. Rep. Prog. Chem., Sec. C, Phys. Chem. 89 (1992) 37. [6] A. Plonka, Annu. Rep. Prog. Chem., Sec. C, Phys. Chem. 91 (1994) 107. [7] P.J. Steinbach, A. Ansari, J. Berendzen, D. Braunstein, K. Chi, B.R. Cowen, D. Ehrenstein, H. Frauenfelder, J.B. Johnson, D.C. Lamb, S. Luck, J.R. Moutant, G.U. Neinhaus, P. Osmos, R. Philipp, A. Xie, R.D. Young, Biochemistry 30 (1991) 3988. [8] R.D. Young, H. Frauenfelder, J.B. Johnson, D.C. Lamb, G.U. Neinhaus, R. Philipp, R. Scholl, Chem. Phys. 158 (1991) 315. [9] H. Frauenfelder, P. Wolynes, Physics Today 47 (1994) 58. [10] V.I. Goldanskii, F. Parak, Chem. Phys. Lett. 239 (1994) 379. [11] W. Primack, Phys. Rev. 100 (1955) 1677. [12] W. Primack, J. Appl. Phys. 31 (1960) 1525. [13] A.I. Mikhailov, Ya.S. Lebedev, N.Ya. Buben, Kinet. Catal. 5 (1964) 1020. [14] A.I. Mikhailov, A.I. Bolshakov, V.I. Goldanskii, Ya.S. Lebedev, Fiz. Tverd. Tela 14 (1972) 1172.
40
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
[15] V.I. Goldanskii, M.A. Kozhushner, Soviet Phys. Chem. Dokl. 274 (1984) 83. [16] J. Jortner, B. Pullman, Tunneling, Reidel, Dordrecht, 1986. [17] R.F. Khairutdinov, K.I. Zamaraev, V.P. Zhdanov (Eds.), Comprehensive Chemical Kinetics, Vol. 30, Elsevier, Amsterdam, 1989, p. 85. [18] V.I. Goldanskii, L.T. Trachtenberg, V.N. Flerov, Tunneling Phenomena in Chemical Physics, Gordon and Breach, New York, 1989. [19] Yu.A. Berlin, V.V. Kulakov, V.G. Nikol'skii, Soviet Phys. Chem. Dokl. 215 (1974) 887. [20] K.I. Zamaraev, R.F. Khairutdinov, A.I. Mikhailov, V.I. Goldanskii, Soviet Phys. Chem. Dokl. 199 (1971) 640. [21] M. Tachiya, A. Mozumder, Chem. Phys. Lett. 34 (1975) 77. [22] A. Plonka, Yu.A. Berlin, N.I. Chekunaev, Chem. Phys. Lett. 158 (1989) 380. [23] S.W. Provencher, Comput. Phys. Commun. 27 (1982) 213. [24] S.W. Provencher, Comput. Phys. Commun. 27 (1982) 229. [25] M. Dishon, J.T. Bendler, G.H. Weiss, J. Res. Natl. Inst. Stand. Technol. 95 (1990) 433. [26] B. Ries, H. B~issler, J. Mol. Electronics 3 (1987) 15. [27] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, J. Chem. Phys. 92 (1990) 7540. [28] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, Radiat. Phys. Chem. 37 (1990) 407. [29] W.H. Hamill, K. Funabashi, Phys. Rev. B 16 (1977) 5523. [30] W.H. Hamill, J. Phys. Chem. 82 (1978) 2073. [31] N.I. Chekunaev, Yu.A. Berlin, V.N. Flerov, J. Phys. C: Solid State Phys. 16 (1983) 1569. [32] V.L. Vyazovkin, B.V. Bolshakov, V.A. Tolkachev, Chem. Phys. 95 (1985) 93. [33] M. Scheidler, B. Cleve, H. Blissler, P. Thomas, Chem. Phys. Lett., in press. [34] R. Richert, H. B~sler, J. Chem. Phys. 85 (1986) 3567. [35] K. Hensel, H. B~issler, Adv. Materials for Optics and Electronics 2 (1992) 179. [36] H. Scher, M.F. Shlesinger, J.T. Bendler, Physics Today 44 (1991) 26. [37] Yu.A. Berlin, Mol. Cryst. Liq. Cryst. 228 (1993) 93. [38] D. Kleinfeld, M.Y. Okamura, G. Feher, Biochemistry 23 (1984) 5780. [39] A.B. Rubin, K.V. Shaitan, A.A. Kononenko, S.K. Chamorovsky, Photosynthesis Research 22 (1989) 219. [40] H. Sumi, in: N. Mataga, T. Okada, H. Masuhara (Eds.), Dynamics and Mechanisms of Photoinduced Transfer and Related Phenomena, Elsevier, Amsterdam, 1992, p. 177. [41] K. Cosstick, T. Asano, N. Ohno, High Pressure Research 11 (1992) 37. [42] N. Agmon, J.J. Hopfield, J. Chem. Phys. 78 (1983) 6947. [43] H. Sumi, R.A. Marcus, J. Chem. Phys. 84 (1986) 4894. [44] W. Nadler, R.A. Marcus, Isr. J. Chem. 30 (1990) 69. [45] H. Heitele, Angew. Chem. Int. Ed. Engl. 32 (1993) 359. [46] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, J. Phys. A: Math. Gen. 26 (1993) 791. [47] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, J. Chem. Phys. 100 (1994) 3163. [48] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, Synth. Metals 64 (1994) 171. [49] Yu.A. Berlin, Chem. Phys. 212 (1996) 29. [50] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, Russian Chem. Phys. 12 (1993) 791. [51] N. Agmon, J.J. Hopfield, J. Chem. Phys. 79 (1983) 2042. [52] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, Chem. Phys. Lett. 197 (1992) 81. [53] E. Gudowska-Novak, J. Phys. Chem. 98 (1994) 5257. [54] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, S.F. Fischer, Chem. Phys. 200 (1995) 369. [55] I. Schlichting, J. Berendzen, G.N. Phillips, R.M. Sweet, Nature 371 (1994) 808. [56] M. Lim, T.A. Jackson, P.A. Anfinrud, Science 269 (1995) 962. [57] S.J. Hagen, J. Hofrichter, W.A. Eaton, Science 269 (1995) 959. [58] D.A. Case, M. Karplus, J. Mol. Biol. 132 (1979) 343. [59] R.F. Tilton Jr., U.C. Sihgh, S.J. Weiner, M.L. Connolly, I.D. Kuntz Jr., P.A. Kollman, N. Max, D.A. Case, J. Mol. Biol. 192 (1986) 443. [60] R.F. Tilton Jr., U.C. Sihgh, I.D. Kuntz Jr., P.A. Kollman, J. Mol. Biol. 199 (1988) 195. [61] R. Elber, M. Karplus, J. Am. Chem. Soc. 112 (1990) 9161. [62] R.G. Palmer, D.L. Stein, E. Abrahams, P.W. Anderson, Phys. Rev. Lett. 53 (1984) 958. [63] M.F. Shlesinger and J. Klafter, in: L. Pietronero, E. Tosatti (Eds.), Fractals in Physics, Elsevier, Amsterdam, 1986, p. 393. [64] N. Eizenberg, J. Klafter, J. Chem. Phys. 104 (1996) 6796. [65] D.K. Lottis, R.M. White, E. Dan Dahlberg, Phys. Rev. Lett. 67 (1991) 362. [66] A.A. Zharikov, S.I. Temkin, A.I. Burshtein, Chem. Phys. 103 (1986) 1. [67] A.A. Zharikov, S.F. Fischer, Chem. Phys. Lett. 249 (1996) 459. [68] A. Erd61yi (Ed.), Tables of Integral Transforms, Vol. 1, McGraw-Hill, New York, 1954.
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
41
[69] G.A. Korn, T.M. Korn, Mathematical Handbook for scientists and engineers, McGraw-Hill, New York, 1968, p. 594. [70] A. Ansari, J. Berendzen, S.F. Bowne, H. Frauenfelder, I.E.T. Iben, T.B. Sauke, R. Shyamsunder, R.D. Young, Proc. Natl. Acad. Sci. USA 82 (1985) 5000. [71] Yu.A. Berlin, S.F. Fischer, Phys. Rev. E, to be submitted. [72] Yu.A. Berlin, A.L. Burin, Chem. Phys. Lett. 267 (1997) 234. [73] I.S. Gradshtein, I.M. Ryzhik, Tables of Integrals, Series and Products, Academic Press, New York, 1980. [74] A. Erd61yi (Ed.), Higher transcendental functions, Vol.3, McGraw-Hill, New York, 1955, p. 206. [75] H. B~issler, Phys. Stat. Sol. (b) 107 (1981) 9. [76] R. Zwanzig, Proc. Natl. Acad. Sci. USA 85 (1988) 2029. [77] Yu.A. Berlin, A.L. Burin, Chem. Phys. Lett. 257 (1996) 665. [78] F. Post, W. Doster, G. Karvounis, M. Settles, Biophys. J. 64 (1993) 1833. [79] N. Agmon, W. Doster, F. Post, Biophys. J. 66 (1994) 1612. [80] W.D. Tian, J.T. Sage, V. Srajer, P.M. Champion, Phys. Rev. Lett. 68 (1992) 408. [81] J.B. Johnson, D.C. Lamb, H. Frauenfelder, J.D. Miiller, B. McMahon, G.U. Nienhaus, R.D. Young, Biophys. J. 71 (1996) 1563.
Chemical Physics 220 (1997) 25-41
Phenomenological model for reaction kinetics coupled to a relaxing environment Yuri A. Berlin a, * ,1, A l e x a n d e r L. B u r i n b, Sighart F. F i s c h e r a a Technische Universiti~t Miinchen, Physik Department T-38, D-85748 Garching, Germany b Department of Physics, CiO' College of the CiO' Universi~ of New York, 138-th Street at ConventAvenue, 10031 New York, NY. USA
Received 17 March 1997
Abstract Reasoning from the correlation between the height of the reaction barrier and the energy depth of substates attributed to different configurations of the reactant surroundings, the phenomenological approach to the rate processes in complex media has been proposed. The approach allows the formulation of the model based exclusively on the energetic consideration of the rate process rather than on the introduction of a particular configurational coordinate. This enables us to make conclusions concerning the most general features of rate processes in complex systems. In particular, we have identified three kinetic phases typical for chemical reactions proceeding in the complex molecular environment, including polychromatic (inhomogeneous), dispersive, and homogeneous phases. The expressions for the temporal behavior of the survival probability in each of these phases have been derived and compared with the empirical stretched exponential law. The temperature and scaling properties of survival probability has also been discussed. © 1997 Elsevier Science B.V.
1. Introduction It has become universally accepted that the classical phenomenology of chemical kinetics breaks down in complex systems like glasses, polymers and biomolecular objects (see e.g. Refs. [ 1-10]). A convincing body of experimental evidences shows that contrary to the classical phenomenology, the reaction patterns for complex systems cannot be adequately described in terms of the familiar time-independent rate constants (for review, see Refs. [2,4-6]); instead, the decay curves can be reproduced only by introducing the specific reaction rates which were found to decrease with time, t, following the near power-law dependence. The latter feature of rate processes in complex molecular media gave the reaction kinetics the name dispersive. The physical origin of the kinetic dispersion has received much attention in recent decades. At the early stage of investigations the dispersive reaction patterns were thought to arise only due to statistical reasons [11,12]. In particular, it has been assumed that a chemical reaction is characterized by a spectrum of activation energies,
* Corresponding author. I On leave from N. N. Semenov Institute of Chemical Physics, Russian Academy of Sciences, ulitsa Kossygina 4, 117334 Moscow, Russian Federation. 0301-0104/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S0301-0104(97)00120-1
26
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
rather than by a single value of this kinetic parameter (the so-called polychromatic concept). Then the decay curves, say, for a first-order thermally activated process, is easy to obtain by averaging the exponential time dependence for the survival probability of an individual reactant over the spectrum of reaction barrier heights (see e.g. Refs. [13-15]). If, in addition, the tunneling mechanism of reaction is considered [16-18], one should take into account the distribution of reaction barriers versus not only heights but also widths (for details, see e.g. [19-21]). As a result, decay curves calculated for first-order reactions turn out to be non-exponential in time both for thermally activated and for tunneling processes. This can formally be interpreted as a manifestation of the temporal evolution of the specific reaction rate, which was shown to proceed in agreement with the dispersive kinetic law, in the cases of certain distribution functions [2,21-25]. The key point in all early formulations of the polychromatic concept is the assumption that the distribution of reaction barriers over their heights and/or widths remains unchanged throughout the whole process. As has been recognized slowly, in view of more recent experimental data, such an assumption is restricted in its validity. For instance, the detailed studies of small ligand binding to heme proteins [7-9] lead to the conclusion that the barrier height distribution does not evolve over the time range (10 ns to 1 ks) as long as the temperature is below 160 K, while above 160 K it shifts as a whole to larger values of barrier heights without markedly changing the shape. Many other experimental evidences supporting the temperature variations of the barrier height distribution for different rate processes can be found in the review [6]. The experimental findings mentioned above indicate that statistics alone cannot provide the exhausted explanation of the kinetic dispersion observed for chemical reactions in complex molecular systems. Therefore significant efforts were made from the theoretical side to gain a deeper insight into the factors responsible for the temporal evolution of the specific reaction rates. Fairly fast progress has been achieved for elementary processes controlled by transport of active species in glasses and solid matrices. A short list of examples, which illustrate a success of theory, includes geminate and homogenous recombination of charged species [26-28], their scavenging by various acceptors [29-31], reactions of radicals [2] and H atoms [32], as well as the excitation transfer [33-35]. Reactants involved in these processes approach each other by successive transitions between structural traps randomly located in the bulk of the system (the bond disorder) or distributed in energy depths E (the so-called site disorder). Since the time of jump exponentially increases with E and with the 'trap-to-trap' distance, the transition rates for different pairs of trapping sites vary within a very wide range, making the resulting diffusion motion highly dispersive [36]. In turn, the latter circumstance leads to the dispersion in the specific reaction rate as verified by results of computer simulations [26,33] and analytical theory [27,28,31,37]. This enables one to conclude that the inherent structural disorder of complex molecular systems is a principle factor determining the dispersive kinetic behavior of active species in the course of diffusion-controlled processes. However available experimental data provides reason enough to believe that the temporal evolution of the specific reaction rate is also typical for the kinetic regime of chemical reactions, where reactants approach each other faster than a chemical event occurs [7-9,38-41]. It has been suggested that in the kinetic regime the dispersion can result from the coupling of the barrier crossing to the stochastic dynamics of the reactant surroundings (see e.g. [40,42-44]; for review, see [45] and references therein). A general framework for treating kinetic consequences of the coupling is provided by a theory of irreversible random transitions [46], which has been invoked earlier for studying randomly affected unimolecular reactions [47], elementary steps of hopping transport [48], the description of hierarchically constrained dynamics of configurational rearrangement in complex systems [49], and other rate processes [50]. Nevertheless, it should be borne in mind that the application of the theory to the situation under discussion calls for a concrete model of the coupling. A conventional approach used for the formulation of the model required is based on the introduction of the one-dimensional (i-d) stochastic process along the reaction coordinate distinct from the distance between reactive species (see e.g. [40,42,43,49,51-54]. If the coordinate is attributed to a configurational variable of the reactant surroundings x [49,51-54], the 1-d motion can be described in terms of sequential transitions between trapping sites corresponding to local minima of potential energy profile along the x-axis. Such minima are
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
27
associated with configurational substates of the chemically inert complex medium [49,52-54]. Taking into account that the substates should have a spread in their energy, again the problem can be formulated in terms of diffusion in disordered systems. Moreover, since a specific reaction rate depends on x due to the dependence of the reaction barrier on the configuration of the reactant surroundings, while x evolves with time because of diffusion, the rate process turns out to be coupled to the dynamics of configurational rearrangements of the molecular environment. It should be mentioned that the model of the stochastic impact considered above implies the knowledge of the specific degree of freedom which controls the reactivity of species in each particular system and over each particular time scale. So far such information remains unavailable. This calls for a more 'flexible' model for studying the effect of configurational rearrangements on the rate processes in complex molecular environment. A demand for the model increases in the light of recent experimental observations [9,55-57], molecular dynamic simulations [58-61], and theoretical findings [49,62-64] which clearly demonstrate that rate processes in complex systems proceeding over a wide range of time scales are in fact controlled by several degrees of freedom rather than by a single configurational coordinate. In this paper we present a phenomenological approach to the rate processes in complex media based on the analogy between the dynamics of their structural rearrangements and the transport phenomena in site-disordered systems. Our approach presumes the correlation between the height of the reaction barrier and the energy depth of substates (traps) attributed to different configurations of the reactant surroundings. This assumption has recently been justified theoretically for a particular case of relaxation in the system of interacting dipoles [65]. Our approach allows the formulation of the model based exclusively on the energetic consideration of the rate process rather than on the introduction of a particular configurational coordinate. This advantage enables us to draw the most general conclusions concerning the kinetic features of rate processes in complex systems and to show that the temporal dispersion in specific reaction rates indeed can arise at the kinetic stage of chemical reaction before the relaxation of the reactant surroundings is completed. The paper is organized as follows. Section 2 contains a brief discussion of our approach to the description of rate processes in relaxing complex media. In this section we also formulate the model by treating the population of configurational substates in terms of the master equation which represents a discrete version of the Feller integral equation commonly used for solving various kinetic problems (see e.g. [66,67]). The solution of the master equation presented in Section 2.2 provides the expression for the Laplace transform of the observable, i.e. of the survival probability. Using the expression derived, we demonstrate in Section 2.3 that this quantity has distinct temporal behaviors within different time and temperature ranges. The explicit form of the behaviors are specified in Section 3 using the exponential and Gaussian distributions as examples. The results obtained are discussed in Section 4, where conclusions followed from our findings are also summarized.
2. Model 2.1. A p p r o a c h and mathematical f o r m u l a t i o n
Here we consider a situation where the ensemble of complex objects (e.g. proteins) is subjected to the excitation leading to the formation of reactants inside each object. This event is supposed to trigger the rearrangement of the chemically inert part of a complex object proceeding through a large number of intermediate configurational substates. Any one of these substates, say i, is modeled as the trap with the energy depth E i, so that the lifetime of the system with respect to the thermally activated release of the i-th substate is given by
ri = Z w °
~B ~ / ,
28
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
where Z is the number of substates distinct from those associated with the i-th trapping site, w 0 is the frequency factor, k B is the Boltzmann constant and T is the temperature. Note that in Eq. (1) the subscript i labels the substates rather than their energies. These latter may have the same value for several substates with different i due to their energetic degeneracy. Qualitatively it is apparent that thermally activated transitions between sites of distinct energies lead to the gradual increase of the population of traps with the larger energy depth. As results, the reactant surroundings relaxes towards more stable configurations. Reasoning from the theoretical results obtained in the particular case of interacting dipoles [65], let us assume now that the deeper is the trap, the higher is the reaction barrier s i. Then for reactions proceeding due to the thermal activation, the rate coefficient can be written as
ki = koexp( _ +
) = ko exp( _ k--~ rlEi ) '
(21
where k 0 is the preexponential factor and r/ is the numerical coefficient. In Eq. (2) the particular value r / = 1 corresponds to the situation where the barriers for the reaction and for configurational rearrangements of the molecular environment have the same heights. In what follows 7/ is assumed to be greater or equal to unity for the sake of definiteness. The opposite case (i.e. "O< 1) can be analyzed in a similar way if one puts T1 = T/~7. Thus, the present model suggests that the relaxation of the molecular environment has influence on the reaction. Such an effect arises since the specific reaction rate depends on the reaction barrier e according to the Arrhenius law (2), while e evolves with time t due to the increase in the population of the deeper substates. Our main concern is the description of reaction kinetics coupled to the relaxation mentioned above. In order to achieve this objective, it is useful to note that the quantity recovered in typical experiments is the probability P(t) of the system remaining in the reactive state at time t (the so-called survival probability). By definition, P(t) is given by
P(t) = ~,pi(t),
(3)
i
where Pi(t) is the probability of the i-th substate to be populated at time t. The latter quantity satisfies the equation
dpi (t) dt
kiPi(t) + ~"(wiiPj(t) -wiJPi(t)) (4) J and the initial condition p~(t = O) = pO. Here Wi~ stands for the conditional probability, per unit time, of the transition from the state j to state i, while the symbol )Z' denotes the summation over all j 4: i. Eq. (4) represents the discrete version of the Feller integral equation commonly used for solving various kinetic problems (see e.g. [66,67]). For uncorrelated transitions, one can assume that 1
1
wii = wi = Z'(i and wi i = w i - Z'ci ,
(5)
where 7 t is given by Eq. (1). For such an assumption, Eq. (4) can be rewritten as ! dPi(t) dt = - ( Z w i + ki)Pi(t) + y" wjpj(t).
(6)
J
The same equation follows from the general analysis of the uncorrelated hopping in the particular case where the waiting time distribution is the exponential function of time with the relaxation rate dependent on the random variable [67]. Once the solution of Eq. (6) is known, the calculation of the observable, i.e. P(t), for any initial probability density pO is straightforward and reduces to the proper summation in Eq. (3).
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
29
2.2. Solution of the kinetic equation and the Laplace transform of the survival probability The problem formulated in Section 2.1 can be solved by introducing the Laplace transform
(7)
p,( s) = fo+~Pi( t) exp( - st)dt , which enables us to write the population of the i-th substate in the following form [68] pO 1 fii(S)
s+Zw
+
i4-k i
E'wj~i(,).
s4-Zw i4-k i j
(8)
Then the Laplace transform of the survival probability can be written as pO 1
fi( s) = E f i i ( s) = E i
+ E
i s4-Zwi4-ki
E'wjfij( s) •
i S j + Z w i 4- k i
(9)
In order to derive the final expression for F'(s), all we need to do is to calculate the sum in the second term of Eq. (9). For Z >> 1, this can be done multiplying both sides of Eq. (8) by w i and by a subsequent summation. Such a procedure gives
E L, J
pA'-"w/Is)=
wi pO
i s + Zw~ + k~ 1 -- ~_~ wi
(10)
i s4-Zwi4-ki
The combination of Eq. (9) with Eq. (10) and the subsequent replacement of the summation over i with the summation over energy of substates allow the continuous representation of the final result for the Laplace transform of the survival probability, viz.
/~(S)
~,w°(O
1/fgl(S) 4- I/F0(S ) , S + - -
(11)
,/,~(s)
where ~,°(s), ~ ( s ) , ~ l ( s ) and U/to(S) are defined in terms of the initial population distribution G(E) and the density of substates g(E) as ~ O ( s ) = fo~ Z w ( E ) G ( E ) d E
s + Z w ( E ) 4- k ( E ) ' g(EldE o s+Zw(E)+k(E)
(12)
,/tg(s) = f ~
qtxl(s)f~
'
k(E)g(E)dE
o ~ + Zw(E)
+ k(E)
'
(13)
(14)
G( E ) d E
~'o(,) = fo
s + Zw(E) + k ( E ) "
(15)
In order to avoid misunderstanding, we note that in mathematical literature the distribution G(E) is referred to as the frequency function, the probability density for the initial population or as the differential distribution function (see e.g. Ref. [69]).
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
30
The above result can be simplified further for several special cases of physical interest. First we consider the situation of very rigid media that corresponds to the limit where w(E) -- O. In this case Eq. (11) reduces to the integral G( E ) d E /~(s) = f0 s + k'(E--) ' (16) so that the subsequent inverse Laplace transform [68] immediately yields the expression for the survival probability derived earlier in the framework of the polychromatic concept (see e.g. Refs. [47] and [51]) oc
P ( t ) = £ G(E)exp( - k ( E ) t ) d E .
(17)
Another interesting situation relevant to the proteinquake hypothesis proposed by Frauenfelder et al. [70] becomes evident if we introduce the probability ~o(E) of the initial population for a single substate with the energy depth E and rewrite the initial population distribution as G(E) = ¢ ( E ) g ( E ) . In the case discussed in the framework of the above hypothesis the excess of energy dissipated into the media after the formation of reactants is assumed to be so high that q~(E) is independent of E. Under this condition, leg'(s) ~g(s)
~,(s) ~'0(s) '
(18)
grl(S) = ~ k(E)G(E)dE o s + Zw(E) + k(E)
(19)
where
Therefore, from Eq. (11) we have gto(S ) 1- f ~ Zw(E)G(E)dE
1
t5(s)
s + qe'(s)
S---S
Jo s + zw( e) +
(20)
e)
Eqs. (11)-(15) as well as Eq. (20) allow the derivation of the expression for the survival probability in several time limits. For the sake of simplicity, here we shall restrict our consideration to the situation, where Eq. (20) is valid. The analysis of the more general expression (11) will be presented elsewhere [71]. 2.3. Asymptotic time behaviors of the survival probability
As can be shown by expanding Eqs. (15), (19) and (20) in a power series in l / s , the Laplace transform of the survival probability for s >> Zw o > k o can be approximated as fi( s) =
1-
k
,
(21)
where the symbol ( . . . ) signifies the averaging over G(E). For the rate coefficient given by Eq. (2), ( k ) exists for any energy depth distribution which is normalized for E ~>0. Since the inverse Laplace transform of Eq. (21) is straightforward [68], we conclude that for t << 1 / ( Z w o) and for Zw o > k o, the survival probability exhibits the temporal evolution described by the linear function of t e(t) = 1-(k)t. (22) It is worth noting that Eq. (22) coincides with the time behavior predicted by Eq. (17) for t < 1/k o. Such a coincidence is not surprising since Eq. (17) results from the first equality in Eq. (20) as long as s >> Zw o. Thus, according to our model, P(t) follows the polychromatic kinetic law if t << 1/(Zwo). This can be expected since at the very early stage of the chemical reaction the reactant surroundings has no time for significant configurational changes and appears to be 'kinetically' frozen in the initial substates.
Y.A. Berlin et al. / ChemicalPhysics 220 (1997)25-41
31
Another asymptotic behavior of P(t) is expected at very long times, where s << min{geo(O)/(1/k2), gel(0)/(1/k)}. Now Eq. (20) becomes 1
fi(s)
gel(0 ) , s + - -
(23)
geo(O)
so that the inverse Laplace transform yields [68] P ( t ) = exp(
ge0(0 )gel(0) t) .
(24)
The analysis of the borderline case where 1 / ( Z w o) < t < 1/min{geo(O)/(1/k2), gel(O)/(1/k)} can be performed by invoking a method proposed in Ref. [72]. Following this method, we introduce the energy-dependent parameter 1
e(E)=
dlnG(E)
"
(25)
dE Then one can verify that the integrand in Eq. (15) will have a maximum if E satisfies the equation
e) - -
+ Zw( e) + k( e) Zw(e) + nk(E)
-
(26)
Based on physical arguments, G(E) is usually assumed to decrease with E as the energy depth of configurational substates becomes equal to or higher than a certain value E. >~0. Therefore e ( E ) > 0 for E >~E . . Besides, in the most commonly encountered situations e ( E ) does not increase with the energy depth when E >/E.. As a consequence, Eq. (26) has a solution only if T is below a certain value Tcr
e( E. ) Zw o + Tlko kB
(27)
Zwo + k o
The latter circumstance allows the conclusion that Eq. (15) gives distinct dependencies ge0 versus s in the highand low-temperature limits. Hence the time behaviors of the survival probability within these two limits should also be distinct. If the temperature is high and it exceeds Tcr, the integrals (15) and (19) approximately equal ge0(0) and gel(0), respectively. Therefore at high temperatures our result for P(t) is identical to Eq. (24). Otherwise (i.e. for T < e ( E . ) / k s ) , the values of these integrals are mainly determined by the energy depths E -= Emax(S) adjacent to the position of the maximum prescribed for the integrand in Eq. (15) by the solution of Eq. (26). Following Ref. [72], ge(s) and gel(s) in this limit can be evaluated by approximating G(E) in the vicinity of Emax(S) as G ( E ) = G(Emax(s)) exp
(28)
kaTmax
with AE = E - Emax(S) and kBTm,x = e( Em,x( S)). Such an approximation together with the introduction of the new variable z = A E / ( k B T ) enable one to show (for details, see [721) that for 77= 1 and for temperatures below Tmax
a(Emax(s))kBZf + ~ e x p ( - T z / T m a x ) d z ge°(s) =
s
.1_~
]'+ exp(-z)
=
G(Emax(s))k.T rr s sin(TrT/Tm.x)
(29)
32
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
Here the second equality sign corresponds to the result of integration given in Ref. [73], and Emax(S) is roughly given by Emax(s)
kaT In
~
s
t
.
(30)
Since in the particular case r/= 1 the integral (19) can be rewritten as ap,(s)
k0 - - [ 1
--saP(s)],
k o + Zw o
(31)
the substitution of Eqs. (29) and (31) in Eq. (20) yields "B'kBZa( Emax( S) )
[
P( s) =
Zw ° ]. ko sin('rrT/Tm,x) + _ _ r r k J G ( Emax( S) ) s ko + Zwo ko + ZWo
(32)
The same ansatz is applicable to the situation where 77> 1 and the reaction barrier exceeds the energy depth of configurational substates. Taking into account that for Z >> 1 the preexponential factor k 0 is much less than Zw o, one may neglect k(E) in the denominator of Eq. (15). Besides, it is evident that s < k o + Zw o within the time range under consideration and therefore aPt(s) can be replaced by aPl(0). As a consequence, for 7/> 1 we obtain
"l'¢kBZa( Emax( S) ) P( s) = s[rrG( Emax( s) )ksT + aPl(O)sin(rrT/Tmax) ] .
(33)
Note that this expression reduces to Eq. (32) as 7/ tends to 1 and k 0 remains less than Zw o. Thus, the inverse Laplace transform of Eq. (33) together with Eqs. (17), (22) and (24) allows the description of the asymptotic behavior of the survival probability within three characteristic time ranges 0 < t << 1/(Zwo), 1/(Zw o) < t < 1/min{aP(0)/( 1/k 2 ), aPl(0)/( 1/k)}, and t >> 1/min{aP(0)/(1/k 2 ), aPl(0)/( 1/k)}, respectively, if the initial population distribution of configurational substates G(E) is known. The explicit form of P(t) for these time ranges will be specified in the next section using the exponential and Gaussian distributions as examples.
3. Examples 3.1. Exponential energy distribution of substates In this case one should put G(E)=
{l
"~0exp -
1
= kRToexp -
0
forE>/0,
(34)
for E < 0,
where E 0 is the parameter of the distribution expressed here in terms of the effective temperature T0. Then the substitution of Eq. (34) into Eq. (17) demonstrates (see also Ref. [73]) that in the short-time limit (t << 1/Zw o) the survival probability can be expressed in terms of the incomplete gamma function T(a,x) and varies with t non-exponentially, namely
P ( t ) = --~o(kot)
~o3'
~7---~o,kot •
(35)
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
33
Using Eq. (21) and properties of the incomplete gamma function [73], one can verify that for k 0 < Zw o and for the very early stage of the reaction, where t << 1/(Zwo), Eq. (35) reduces to Eq. (22) with the mean rate coefficient given by T ( k ) = k0 r
(36)
+ tiT0
In order to establish the temporal behavior of the survival probability at t > 1/Zw o in the high-temperature limit T > T0, we need now to calculate the ratio ~ 1 ( 0 ) / ~ ( 0 ) . For k 0 < Zw o, this ratio coincides with the ratio ~ l ( 0 ) / q t ( 0 ) as long as T exceeds T0. This implies that except for very short times t << 1 / Z w o, the reaction kinetics remains exponential and follows the familiar exponential law P ( t ) = exp( - ( k ) t ) .
(37)
However the results for T < To are different. The substitution of Eq. (34) into Eqs. (15), (19) and (20) demonstrates that at low temperatures for s >> Zw o ,
qto(s) fi( S) =
(38)
S - I ( oIS) fl
for s < Zw o ,
1+
with parameters a and /3 being defined as
ol = Zwo
koTo sin('rrT/To)
1 + ( 7 7 - 1)
,
T /3 = - - < 1.
(39) (40)
ro
The inverse Laplace transform of Eq. (38) (see Ref. [68]) makes evident that in the low-temperature limit ( T < T0) the temporal behavior of the survival probability is non-exponential. For t << 1 / Z w o, P(t) can be described by Eq. (35) as in the high-temperature limit, while for t > 1/Zw o this quantity becomes equal to the Mittag-Leffler function E t 3 ( - ( t / a ) t J ) . By definition [74], the latter can be written as
( - 1 ) k ( t / a ) ~k E t 3 ( - ( t / c t ) t3) = • k=0
F(1 +ilk)
,
(41)
where F ( x ) is the gamma function. Since in the long-time limit E t ~ ( - ( t / o t ) ~ ) has the asymptotic expansion [74] =
-
E n=0
+o
we conclude that
1 r(l +/3) P(t) = I (t/ot)_t3 r(1 - / 3 )
for t<< c~, (42) for t>> a .
34
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
It is interesting to compare the time behavior of the survival probability following from Eq. (41) with the empirical stretched exponential law
commonly used for the description of the reaction patterns in terms of two fitting parameters o/1 and ¢tl (for review, see e.g. Refs. [2,4], and [6]). Numerical data presented in Fig. 1 allows the comparison between Eq. (41) and Eq. (43) for different values of the dispersion parameter /3. As can be seen, in the case where a 1 = a [ F ( 1 +/3)] - l / ~ and ~ =/3 the temporal behaviors predicted by these equations coincide only at the initial part of the decay curves as long as t is less than a few tenths of a. This can be expected reasoning from the short-time asymptotic behavior of the survival probability given by the upper line of Eq. (42). If, however, oq ¢ ~ and /3t 4=/3, the stretched exponential provides a reasonable approximation of Eq. (41) in a wider time range. So far our consideration has been focused on the cases, where the temperature is either far above or much below T0. The reaction kinetics in the vicinity of the characteristic temperature To merits a closer examination. 1.0 A 0.8 ¸
0.6 ¸ 0.4 ¸ 0.2 ¸
0.0 1.0 ~ ~
0.8
B
•
0.2
~
0.0 1.0
"o0.61, t , " i ~ ,
, C
0.8 0.6
~
0.4 0.2 0.0 0.0
2.0
4.0
6.0
8.0
10.0
t/a Fig. 1. Temporal behavior of the survival probability P calculated for different values of the dispersion parameter/3. Curves shown in panels A, B, and C are plotted for/3 = 0.2, 0.4, and 0.6, respectively. Solid curves are dependencies given by the Mittag-Leffier function, see Eq. (41). Dashed curves correspond to the stretched exponential given by Eq. (43) with a I = a and /3~ = ft. D o t - d a s h e d curves are the fits to the Mittag-Leffler time behavior. The fitting parameters are a I = 6 a and /31 = 0.75/3 (panel A), oq = 1.9or and fll = 0.75/3 (panel B), oq = 1.3or and /3t = 0.83/3 (panel C).
35
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
The calculations carried out for t << 1/Zw o show that in the short-time limit the time evolution of the survival probability at T = TO can also be described by Eq. (35). For the very early stage of the reaction where t << l / ( Z w o ) , the incomplete gamma function can be approximated as y(1/T1, kot) = 71(kot)t/n[1 - kot/(~l + 1)] if k 0 < Zw o [73]. Hence, in this case Eq. (35) becomes kot - -
P(t) = 1
- 1 -(k)t,
1+77
(44)
in agreement with Eq. (22). For longer times t >1 l / ( Z w o) and for Zw o > k o, the evaluation of the integral (15) yields 1
(
~ ° ( s ) = oZ w - l n
1+
~Wo) s '
(45)
while ~l(s) may still be approximated by ~1(0) = ko/(~lZwo). As a consequence, from Eq. (20) one obtains fi(s)---~ln
l+
s
"
(46)
Since the inverse Laplace transform of Eq. (46) is straightforward [68], we arrive at the asymptotic expression for the survival probability ~/~ (1 - e x p ( - Z w o t ) ) . P ( t ) = ko---
(47)
Eq. (47) is valid within the time interval 1 / ( Z w o) <~t <~( 1 / Z w o ) e x p ( T / [ T - To[) as long as the temperature is close to To. In addition, it reproduces the temporal behavior of the survival probability given by the lower line of Eq. (42) which holds for t > a even in the case where I T - To] << To. 3.2. Gaussian distribution
Another energy distribution of configurational substates frequently used in the literature [19,26,28,37,38,51] is defined as 1 G(E)
(E-
e0) 2
4tr z
=
0
for E>~ 0, (48) for E < 0 .
This is the so-called Gaussian distribution dependent on two parameters 60 and tr which predetermine the position of the maximum and the spread of substates in their energy depths, respectively. The preexponential factor in Eq. (48) arises due to the normalization condition o~
fo G( E ) d E = 1
and contains the error function [73] 2
eft(x)
=
x
,rrl/2 fo exp(-~:Z)d~:"
In order to derive the explicit form of P ( t ) for the above distribution in different time and temperature
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
36
ranges, we introduce the dimensionless parameter a = c0/(2 tr) and express the width of the Gaussian tr in terms of the characteristic temperature To = cr/k B. Then the mean rate coefficient becomes
(k)=k0ex p
1-
whereas the ratio ~(0)/1//1(0)
~(0)
k°exp
can
,+e~+_ ,
(49)
be written as
t1
r
,
2AT
]]
[( '
1 + erf A 1 + ~-~
(50)
It should be mentioned that for temperatures exceeding a certain value rtTo
for a = 0 ,
nTo
(51)
for a > 0,
a
Eqs. (49) and (50) provide the identical result, namely
,,,,
~(0)
(
for,,=O,
kBT]
k0ex p
(52)
for A > 0 .
This implies that the exponential and Gaussian distributions give the same results qualitatively in the high-temperature limit: The survival probability exhibits the exponential temporal behavior with the exception of the very early stages of the process, where P(t) decays non-exponentially following either Eq. (17) or its particular form (22) as long as t << 1/(Zwo). In the opposite case (i.e. at T < TG) the survival probability exponentially decreases with time only in the very late stage of the reaction (t>> 1/min{qt(O)/(1/k2), ~l(O)/(l/k)}). According to our theoretical analysis (see Section 2.3), the rate coefficient characterizing the temporal behavior of P(t) in this stage is given by Eq. (50). The latter expression allow the further simplification for two particular positions of the Gaussian distribution at the E-axis. If the Gaussian peaks at E = 0, one should put A = 0 in Eq. (50). This leads to the following expressions for qtl(0)/~(0) in the low-temperature limit
qt,(O)
~(0)
--}- e x p , - - ~ ]
k°T
for r/= 1,
7"°2
. . . . 7 (__~)] exp - -~-g "rrl/2(l - r / ) T o 1 + e f t
(53) for'q > 0.
EA. Berlin et al. / Chemical Physics 220 (1997) 25-41
37
In the other situation which corresponds to the relatively narrow Gaussian distribution with A/+/> 1 and with the maximum centered at E = e0 4=0, we get
( ( 2a7"0 // 0exp/
1/fl ( 0 )
-/_
T2 )
2"rrl/2(1 -- ~)T0 exp
f o r , = 1,
(
T
T2
for ~7> 0.
The comparison of Eq. (53) and Eq. (54) shows that the shift of the Gaussian distribution towards higher E values results in the increase of the reaction activation energy in the late stage of chemical process if it proceeds at T < TC. In addition, Eq. (53) and Eq. (54) evidence that the reaction kinetics in the long-time limit does not obey the Arrhenius law. The latter circumstances can be explained by the attaining of the dynamic equilibrium between the population of configurational substates and their release in the course of the relaxation process. The similar explanation has been proposed earlier for the non-Arrhenius temperature dependence of the diffusion coefficient for the anomalous transport in site-disordered systems [75-77]. In order to complete the case under consideration, we now turn our attention to the intermediate time range 1/(Zw o) < t < 1/min{~(O)/(1/k2), qtl(O)/(1/k)}. As follows from the general results presented in Section 2.3, the temporal evolution of the survival probability within this time range can be deduced by taking the inverse Laplace transform of Eq. (33). For the Gaussian distribution (48), the use of Eqs. (25) and (30) allows the conclusion that parameters G( Emax(S)), Tmax, and qq(0) appearing in Eq. (33) can be approximated as exp( -A2 )
( s )~
G(Emax(S))= "rrl/2o'[1Terf(A)] ~ w0
(55)
'
½T Tmax
I
(56)
AT '
/3+s-7. °
qtt(0)=
k°
[
Zwo exp
2(r/-1)AT°( ~ 1
(+/-1)T°)] 2-AT
t,
(57)
1 + eft(A)
with /3 only slightly dependent on s, namely /3=
In
To
The substitution of Eqs. (55)-(58) into Eq. (33) and the subsequent inverse Laplace transform [68] allow the derivation of asymptotic expressions for P(t) similar in form to Eq. (42), but with the characteristic time c~ and the dispersion parameter/3 given by 1 [
"rrl/2Texp(-A 2)
o~= ZWo a/Zl(0)T0(1 +erf(A))sin(.rrT//Tmax ) ( T):
/3= -~o
ln(Zw°t)-z
AT
]1/3 .
(59) (60)
This implies that the dispersive reaction patterns can be expected both for the exponential and for the Gaussian energy distribution of configurational substates, as long as kinetic measurements are performed within the time
38
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
range l/(Zw o) < t < 1/min{q~(O)/(1/k2), x/rj(0)/(1/k)} and at temperatures lower than a critical value T = T *. The latter is related to parameters of distribution functions. In the case of the exponential distribution (34) T * = T0, while for the Gaussian (48) the critical temperature T * equals T~ and has the explicit form given by Eq. (51).
4. Discussion and conclusions
In the present paper we propose a phenomenological model for rate processes proceeding in the chemically inert environment relaxing towards the dynamic equilibrium. The equilibrium arises as a result of two competing processes, i.e. the repopulation of trapping sites associated with the configurational substates of the reactant surroundings and their release due to the thermal activation. The model suggests the correlation between the height of the reaction barrier and the energy depth of traps corresponding to the configurational substates: the deeper traps are, the higher is the barrier. Such an assumption was justified theoretically for the particular case of relaxation in the system of interacting dipoles [65]. The main physical consequence of the correlation mentioned above is the coupling of the chemical reaction to the relaxation process. Starting with the pioneering work of Agmon and Hopfield [42,51], the influence of the coupling on the reaction kinetics is theoretically treated in the framework of the approach based on the adequate selection of the specific configurational coordinate controlled the rate process in each particular system. Since the experimental information needed for making the proper choice is inaccessible in most cases, the application of the approach to the description of reaction kinetics in concrete complex systems is usually hard to justify. By contrast, a model presented here treats the chemical reaction coupled to the relaxation of the molecular surroundings in terms of the only variable, namely the energy of configurational substates E. This enables us to avoid difficulties related to the proper choice of the configurational coordinate and to draw the most general conclusions concerning the kinetic features of rate processes in complex systems. According to our theoretical findings, such features include (i) the existence of the polychromatic phase of rate processes [ 11-15] known in the biophysical literature as inhomogeneous kinetics (see e.g. Refs. [78] and [79]) arising as long as t << 1/(Zw o) or if physical conditions prevent conformational rearrangements of the surroundings (very low temperatures, rigid media, etc.); (ii) the transition from the polychromatic to the dispersive phase; the advent of the latter signals the onset of the relaxation process associated with configurational changes in the reactant surroundings; (iii) the crossover to the exponential kinetics in the long-time limit associated with the completion of relaxation; (iv) the non-Arrhenius temperature dependence of mean rate coefficient for the exponential phase of the rate process (cf. Eqs. (53) and (54)); (v) the existence of the critical temperature T * dependent on the parameters of the energy distribution of the configurational substates, serving as a demarcation line between exponential and dispersive kinetics in experiments performed in a wide temperature range and within the fixed time window of the appropriate position and width; (vi) the scaling properties of the survival probability in the two asymptotic time limits which are shown to exist at temperatures below T * prior to the completion of the relaxation process (see the upper and the lower lines of Eq. (42)). In addition, the proposed model enables us to derive expressions for the observable, viz. for the survival probability P(t), which describe its temporal evolution during each of kinetic phases mentioned above. In particular, we show that in the dispersive phase P(t) can be expressed in terms of the Mittag-Leffler function E~(-(t/a)~), c.f. Eq. (41). Our results shown in Fig. 1 clearly demonstrate that the empirical stretched exponential law (see Eq. (43)) commonly employed for fitting kinetic data in the case of the first-order reactions (for review, see e.g. Refs. [2,4,6]), exhibits practical the same temporal behavior as Et~(-(t/a)~) within the
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
39
restricted time range 1/(Zw o) < t <~ 10a. In order to attain the coincidence, one should assume, however, that parameters of the stretched exponential decay are distinct in their values from those used in numerical calculations of Et3(- (t/oz)13). This, of course, does not provide arguments for deciding between the stretched exponential and the Mittag-Leffier functions. Nonetheless, our numerical findings clearly show that the use of the stretched exponential, for the description of rate processes in complex systems, can be justified only by comparison with kinetic curves measured within a wide time window and reflecting the temporal evolution of the process up to high degrees of the chemical conversion ( P ~< 0.1). Otherwise, the fitting procedure invoking this function may lead to the values of dispersive parameter and the characteristic relaxation time being far from real. There are grounds to believe that chemical reactions in real complex systems indeed proceed in three phases expected from the proposed theory. This is especially true for processes in biological objects, in particular for ligand binding to heme proteins. In the latter case multiple pulse experiments [78-80], kinetic hole burning [7], and flash photolysis with monitoring in the near infra red [81] strongly support the existence of these kinetic phases. Data obtained by these methods provide additional information concerning phase durations at different temperatures. Finally, we would like to note that a theoretical analysis presented above is restricted to the case, where the excess of energy dissipated into the media after the formation of reactants is much higher than the mean energy depth of substates. This corresponds to the proteinquake concept proposed by Frauenfelder [70]. Nevertheless, the expression (12) for the Laplace transform of the survival probability allows a more general analysis without invoking special assumptions concerning the amount of dissipated energy. The results for the latter case which will be published elsewhere [71] show that the main conclusions followed from our study remain unchanged, but the temporal behavior of the reactant becomes more diverse in comparison with kinetic features deduced in the present work.
Acknowledgements This work is supported by the Deutsche Forschungsgemeinschaft. YuAB greatly acknowledges the DFG grant No.436 RUS 1 7 / 1 6 / 9 7 . ALB is indebted to the Alfred Toepfer Foundation and to the Netherlands Organization for Scientific Research. The authors would like to thank Dr. A. Zharikov for discussions and Prof. Dr. H. B~issler for reprints of his works.
References [1] Yu.A. Berlin, J.R. Miller, A. Plonka (Eds.), Rate Processes with Kinetic Parameters Distributed over Time and Space, (Chem. Phys. 212 (1996)). [2] A. Plonka, in Lecture Notes in Chemistry, Vol. 40, Springer, Berlin, 1986. [3] H. Frauenfelder, G.U. Nienhaus, J.B. Johnson, Ber. Bunsenges. Physik. Chem. 17 (1988) 451. [4] A. Plonka, Prog. React. Kinetics 16 (1991) 157. [5] A. Plonka, Annu. Rep. Prog. Chem., Sec. C, Phys. Chem. 89 (1992) 37. [6] A. Plonka, Annu. Rep. Prog. Chem., Sec. C, Phys. Chem. 91 (1994) 107. [7] P.J. Steinbach, A. Ansari, J. Berendzen, D. Braunstein, K. Chi, B.R. Cowen, D. Ehrenstein, H. Frauenfelder, J.B. Johnson, D.C. Lamb, S. Luck, J.R. Moutant, G.U. Neinhaus, P. Osmos, R. Philipp, A. Xie, R.D. Young, Biochemistry 30 (1991) 3988. [8] R.D. Young, H. Frauenfelder, J.B. Johnson, D.C. Lamb, G.U. Neinhaus, R. Philipp, R. Scholl, Chem. Phys. 158 (1991) 315. [9] H. Frauenfelder, P. Wolynes, Physics Today 47 (1994) 58. [10] V.I. Goldanskii, F. Parak, Chem. Phys. Lett. 239 (1994) 379. [11] W. Primack, Phys. Rev. 100 (1955) 1677. [12] W. Primack, J. Appl. Phys. 31 (1960) 1525. [13] A.I. Mikhailov, Ya.S. Lebedev, N.Ya. Buben, Kinet. Catal. 5 (1964) 1020. [14] A.I. Mikhailov, A.I. Bolshakov, V.I. Goldanskii, Ya.S. Lebedev, Fiz. Tverd. Tela 14 (1972) 1172.
40
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
[15] V.I. Goldanskii, M.A. Kozhushner, Soviet Phys. Chem. Dokl. 274 (1984) 83. [16] J. Jortner, B. Pullman, Tunneling, Reidel, Dordrecht, 1986. [17] R.F. Khairutdinov, K.I. Zamaraev, V.P. Zhdanov (Eds.), Comprehensive Chemical Kinetics, Vol. 30, Elsevier, Amsterdam, 1989, p. 85. [18] V.I. Goldanskii, L.T. Trachtenberg, V.N. Flerov, Tunneling Phenomena in Chemical Physics, Gordon and Breach, New York, 1989. [19] Yu.A. Berlin, V.V. Kulakov, V.G. Nikol'skii, Soviet Phys. Chem. Dokl. 215 (1974) 887. [20] K.I. Zamaraev, R.F. Khairutdinov, A.I. Mikhailov, V.I. Goldanskii, Soviet Phys. Chem. Dokl. 199 (1971) 640. [21] M. Tachiya, A. Mozumder, Chem. Phys. Lett. 34 (1975) 77. [22] A. Plonka, Yu.A. Berlin, N.I. Chekunaev, Chem. Phys. Lett. 158 (1989) 380. [23] S.W. Provencher, Comput. Phys. Commun. 27 (1982) 213. [24] S.W. Provencher, Comput. Phys. Commun. 27 (1982) 229. [25] M. Dishon, J.T. Bendler, G.H. Weiss, J. Res. Natl. Inst. Stand. Technol. 95 (1990) 433. [26] B. Ries, H. B~issler, J. Mol. Electronics 3 (1987) 15. [27] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, J. Chem. Phys. 92 (1990) 7540. [28] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, Radiat. Phys. Chem. 37 (1990) 407. [29] W.H. Hamill, K. Funabashi, Phys. Rev. B 16 (1977) 5523. [30] W.H. Hamill, J. Phys. Chem. 82 (1978) 2073. [31] N.I. Chekunaev, Yu.A. Berlin, V.N. Flerov, J. Phys. C: Solid State Phys. 16 (1983) 1569. [32] V.L. Vyazovkin, B.V. Bolshakov, V.A. Tolkachev, Chem. Phys. 95 (1985) 93. [33] M. Scheidler, B. Cleve, H. Blissler, P. Thomas, Chem. Phys. Lett., in press. [34] R. Richert, H. B~sler, J. Chem. Phys. 85 (1986) 3567. [35] K. Hensel, H. B~issler, Adv. Materials for Optics and Electronics 2 (1992) 179. [36] H. Scher, M.F. Shlesinger, J.T. Bendler, Physics Today 44 (1991) 26. [37] Yu.A. Berlin, Mol. Cryst. Liq. Cryst. 228 (1993) 93. [38] D. Kleinfeld, M.Y. Okamura, G. Feher, Biochemistry 23 (1984) 5780. [39] A.B. Rubin, K.V. Shaitan, A.A. Kononenko, S.K. Chamorovsky, Photosynthesis Research 22 (1989) 219. [40] H. Sumi, in: N. Mataga, T. Okada, H. Masuhara (Eds.), Dynamics and Mechanisms of Photoinduced Transfer and Related Phenomena, Elsevier, Amsterdam, 1992, p. 177. [41] K. Cosstick, T. Asano, N. Ohno, High Pressure Research 11 (1992) 37. [42] N. Agmon, J.J. Hopfield, J. Chem. Phys. 78 (1983) 6947. [43] H. Sumi, R.A. Marcus, J. Chem. Phys. 84 (1986) 4894. [44] W. Nadler, R.A. Marcus, Isr. J. Chem. 30 (1990) 69. [45] H. Heitele, Angew. Chem. Int. Ed. Engl. 32 (1993) 359. [46] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, J. Phys. A: Math. Gen. 26 (1993) 791. [47] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, J. Chem. Phys. 100 (1994) 3163. [48] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, Synth. Metals 64 (1994) 171. [49] Yu.A. Berlin, Chem. Phys. 212 (1996) 29. [50] Yu.A. Berlin, D.O. Drobnitsky, V.V. Kuzmin, Russian Chem. Phys. 12 (1993) 791. [51] N. Agmon, J.J. Hopfield, J. Chem. Phys. 79 (1983) 2042. [52] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, Chem. Phys. Lett. 197 (1992) 81. [53] E. Gudowska-Novak, J. Phys. Chem. 98 (1994) 5257. [54] Yu.A. Berlin, N.I. Chekunaev, V.I. Goldanskii, S.F. Fischer, Chem. Phys. 200 (1995) 369. [55] I. Schlichting, J. Berendzen, G.N. Phillips, R.M. Sweet, Nature 371 (1994) 808. [56] M. Lim, T.A. Jackson, P.A. Anfinrud, Science 269 (1995) 962. [57] S.J. Hagen, J. Hofrichter, W.A. Eaton, Science 269 (1995) 959. [58] D.A. Case, M. Karplus, J. Mol. Biol. 132 (1979) 343. [59] R.F. Tilton Jr., U.C. Sihgh, S.J. Weiner, M.L. Connolly, I.D. Kuntz Jr., P.A. Kollman, N. Max, D.A. Case, J. Mol. Biol. 192 (1986) 443. [60] R.F. Tilton Jr., U.C. Sihgh, I.D. Kuntz Jr., P.A. Kollman, J. Mol. Biol. 199 (1988) 195. [61] R. Elber, M. Karplus, J. Am. Chem. Soc. 112 (1990) 9161. [62] R.G. Palmer, D.L. Stein, E. Abrahams, P.W. Anderson, Phys. Rev. Lett. 53 (1984) 958. [63] M.F. Shlesinger and J. Klafter, in: L. Pietronero, E. Tosatti (Eds.), Fractals in Physics, Elsevier, Amsterdam, 1986, p. 393. [64] N. Eizenberg, J. Klafter, J. Chem. Phys. 104 (1996) 6796. [65] D.K. Lottis, R.M. White, E. Dan Dahlberg, Phys. Rev. Lett. 67 (1991) 362. [66] A.A. Zharikov, S.I. Temkin, A.I. Burshtein, Chem. Phys. 103 (1986) 1. [67] A.A. Zharikov, S.F. Fischer, Chem. Phys. Lett. 249 (1996) 459. [68] A. Erd61yi (Ed.), Tables of Integral Transforms, Vol. 1, McGraw-Hill, New York, 1954.
Y.A. Berlin et al. / Chemical Physics 220 (1997) 25-41
41
[69] G.A. Korn, T.M. Korn, Mathematical Handbook for scientists and engineers, McGraw-Hill, New York, 1968, p. 594. [70] A. Ansari, J. Berendzen, S.F. Bowne, H. Frauenfelder, I.E.T. Iben, T.B. Sauke, R. Shyamsunder, R.D. Young, Proc. Natl. Acad. Sci. USA 82 (1985) 5000. [71] Yu.A. Berlin, S.F. Fischer, Phys. Rev. E, to be submitted. [72] Yu.A. Berlin, A.L. Burin, Chem. Phys. Lett. 267 (1997) 234. [73] I.S. Gradshtein, I.M. Ryzhik, Tables of Integrals, Series and Products, Academic Press, New York, 1980. [74] A. Erd61yi (Ed.), Higher transcendental functions, Vol.3, McGraw-Hill, New York, 1955, p. 206. [75] H. B~issler, Phys. Stat. Sol. (b) 107 (1981) 9. [76] R. Zwanzig, Proc. Natl. Acad. Sci. USA 85 (1988) 2029. [77] Yu.A. Berlin, A.L. Burin, Chem. Phys. Lett. 257 (1996) 665. [78] F. Post, W. Doster, G. Karvounis, M. Settles, Biophys. J. 64 (1993) 1833. [79] N. Agmon, W. Doster, F. Post, Biophys. J. 66 (1994) 1612. [80] W.D. Tian, J.T. Sage, V. Srajer, P.M. Champion, Phys. Rev. Lett. 68 (1992) 408. [81] J.B. Johnson, D.C. Lamb, H. Frauenfelder, J.D. Miiller, B. McMahon, G.U. Nienhaus, R.D. Young, Biophys. J. 71 (1996) 1563.