J. Phys. Chem. 1981, 85, 1549-1555
1549
Numerical Solution of the Smoluchowski Equation Applied to the Radiolysis of Aliphatic Amines and Hydrazine J. A. Delaire, Laboratoire de Physico-Chimie des Interm6diaires Rgactionnels, L.A 75, Universit6 de Paris-Sud, 9 1405 Orsay Cedex, France
E. Croc, Equipe Toulonnaise de Math6matiques Appliqdes, Universlt6 de Toulon, 83 130 La Garde, France
and P. Cordier” Laboratoire de Chimie Physique des Mat6riaux Amorphes, L.A 75, Universit6 de Paris-Sud, 9 1405 Orsay Cedex, France (Received: July 28, 1980; In Final Form: November 17, 1980)
The time-dependent Smoluchowski equation is solved numerically by variational methods, for general boundary conditions and an initial Gaussian distribution function. The time and space conditions of the stability of the algorithm are discussed. The survival probability and scavenging probability of the solvated electron are computed. The numerical results are then compared with recent experiments on the solvated electron and biphenyl anion yields in amines and hydrazine. For liquids in which the rate constant for recombination between the solvated electron and the cation is high, the model applies and gives reasonable values for the mean initial separation distance t of the electron-cation pairs (t = 40 f 10 a in ethylamine and 50 f 10 A in n-propylamine). For other liquids such as ethylenediamine and hydrazine, it is shown that a competitive reaction of the solvated electron with a radical has to be taken into account to fit the experimental data concerning the nonhomogeneous decay of the solvated electron.
I. Introduction In the condensed phase, ionic species produced either by photoionization or by radiolysis are distributed nonhomogeneously. The fate of such ions has been intensively studied and considerable experimental data are available concerning electron bulk rate constants of ion neutralization,’I2 scavenger s t u d i e ~ , and ~ - ~ electric field e f f e c t ~ . ~ In ~ ~the - ~ case of isolated pairs of oppositely charged ions, theoretical paperss14 deal with the solution of the Smoluchowski equation which adequately describes the potential-dependent brownian motion of ions. The solutions of the Smoluchowski equation published so far differ from each other either by general conditions or by integration methods. Among the general conditions, one usually considers the steady-state and time-dependent problems, the different acting electric potentials, an eventual dissipative term due to the presence of an ionic ~~~~~
(1) J. Belloni, F. Billiau, P. Cordier, J. A. Delaire, M. 0. Delcourt, and 63, 58 (1977), and references M. Magat, Faraday Discuss. Chem. SOC., therein. (2) J. A. Delaire and J. R. Bazouin, Can. J . Chem., 57, 2013 (1979). (3) For a review on free-ion yields, scavenger studies, and electric field effects, see A. 0. Allen, Natl. Stand. Ref. Data Sec., Natl. Bur. Stand., No.57 (1976). (4) For a review on electron yields and scavenger studies in alcohols and water, see J. W. Hunt, Adu. Radiat. Chem., 5, 185 (1976). (5) J. W. Fletcher and W. A. Seddon, Faraday Discuss. Chem. Soc., 6% 18 -- 11977). \--. (6) J. Casanovas, R. Grob, D. Blanc, G. Brunet, and J. Mathieu, J. Chem. Phys., 63, 323 (1975). (7) W. F. Schmidt, Can. J. Chem., 55, 2197 (1977). (8)J. Bullot, P. Cordier, and M. Gauthier, J. Chem. Phys., 69, 1374, 4908 (1978). (9) For a review on models in nonpolar liquids, see A. Hummel, Adu. Radiat. Chem., 4, 1 (1974); (10) (a) G. G. Abell, A. Mozumder, and J. Magee, J. Chem. Phys., 56, 5422 (1972); (b) K. M. Hong and J. Noolandi, J. Chem. Phys., 68,5163 (1978). (11)A. Mozumder, J. Chem. Phys., 69, 1384 (1978). (12) S. A. Rice and J. K. Baird, J. Chem. Phys., 69, 1989 (1978). (13) H. Sano and M. Tachiya, J. Chem. Phys., 71, 1276 (1979). (14)Y. A. Berlin, P. Cordier, and J. A. Delaire, J. Chem. Phys., 73, 4619 (1980). --I
. I .
scavenger, the initial ion pair distribution, and the boundary conditions. In fact, all of these general conditions are not always strictly independent. Thus, the Magee-Tayler theorem15 employs the Laplace transform to correlate the steady-state solution in the presence of a scavenger with the time-dependent solution in the absence of a scavenger. Similarly, from two different sets of both initial and boundary conditions, the solutions can be connected. Thus, the Rice-Baird theorem12 shows that there exists a relationship between the escape probability and the bulk rate constant for ionic recombination. As far as the mathematical aspect of the integration is concerned, it is convenient to separate analytical methods from numerical techniques. Due to the peculiar nature of the Smoluchowski transient equation, the former generally serve to obtain approximate results16 and the latter lead to accurate solution^.^'-^^ Recently, this point has been clearly demonstrated by the comparison between computed and closed form results obtained from the calculation of the time-dependent bulk rate constant of ion neutralization.18 An alternative consists in transforming the original equation by an analytical processlo and then to proceed by a numerical technique. A typical example of such a procedure is given in the Hong and Noolandi paper.1° The distribution function is obtained by a computational method from an expanded solution of the Laplace-transformed transient equation. It should be argued that the analytical procedures provide a useful insight to predict the role of variables and of the associated parameters. However, as recently shown,20the solutions of the (15) J. L. Magee and A. B. Tayler, J. Chem. Phys., 56, 3061 (1972). (16) This is the case of the “prescribed diffusion” method (A. Mozumder, J. Chem. Phys., 47, 1859 (1967); 50, 3153 (1969)). (17) A. Hummel and P. P. Infelta, Chem. Phys. Lett. 24, 559 (1974). (18) S. A. Rice, P. R. Butler, M. J. Pilling, and J. K. Baird, J. Chem. Phys., 70, 4001 (1979). (19) 2. Schulten and K. Schulten, J . Chen. Phys.,66, 4616 (1977). (20) E. Croc, M. C. Pelissier, P. Penel, and P. Cordier, submitted to Mathematical Modelling.
0022-3654/81/2085-1549$01.25/0@ 1981 American Chemical Society
1550
The Journal of Physical Chemistry, Vol. 85, No. 11, 1981
Smoluchowski equation belong to a functional vector space in which most of the elements cannot be expressed analytically. This point is in favor of the use of direct numerical integration methods. Numerical techniques for solving the Smoluchowski transient equation have not yet been intensively reported.17-19 Most of the numerical solutions miss a basic mathematical support in elaborating the numerical procedure. Also, in this paper, the numerical integration is based upon the use of the variational method.21 By introducing a relaxation coefficient, the variational method then encompasses all singular treatments so far proposed. In addition, efforts have been made to obtain a stable algorithm with respect to time and space physical conditions. This aspect is developed in section I1 by abstracting the detailed mathematical analysis of ref 27. Our mathematical procedure has been applied and tested on a simplified physical model which supposes that the dielectric constant is not time dependent and the only reaction of the solvated electron is with a positive charge. Such a model has been largely used in recent years.”14 The numerical treatment of the Smoluchowski equation has been used to obtain the yields of solvated electrons and of scavenging products in the radiolysis of amines2 and hydrazine.22 These liquids have a peculiar interest because, as shown for some of them, the electron-cation recombination is not diffusion controlled. The efficiency of this recombination is included in the model, and then makes possible the comparison between reactions occurring in a spur and those occurring in the bulk. Moreover, the comparison between theory and experiment leads to the determination of the mean value of the initial electroncation distance in the radiolysis of these liquids. 11. Mathematical Analysis and Computational Met hod In an isolated pair, the electron experiences the Coulombic field of the cation, so that the Smoluchowski equation and associated conditions, in spherical coordinates, are written as follows:
+ Rcp(r,t)
1
=0
(3)
p(r,O) = po(r) (4) As usual p(r,t) du is the probability of finding the electron in the volume element dv at time t , D is the mutual diffusion coefficient of the electron and the cation, and Rc is the Onsager distance.23 Equation 2 expresses the first boundary condition and states that the electron flux through a sphere of radius R is consumed by a pseudo-first-order reaction with the cation. In this equation, k, represents the effectiveness of the reactant encounters, i.e., the activation-controlled rate constant, and R is the reaction radius. Usually referred to as the radiation (21)J. L. Lions, “Mathematics Applied to Physics”, Springer-Verlag, West Berlin, 1970,p 348;J. L.Lions and E. Magenes, “ProblBmes aux limites non homogenes et applications”,Dunod, Ed., Paris, 1968. (22)J. A. Delaire, P. Cordier, J. Belloni, F. Billiau, and M. 0. Delcourt, J. Phys. Chen., 80, 1687 (1976). (23)L. Onsager, Phys. Rev., 54,554 (1938).
Delaire et al.
boundary ~ondition,2~ eq 2 takes account of the chemical process and encompasses all singular cases such as in an idealized sink or the effect of unsucessful reactive encounters between species.25 The choice of the second boundary condition is sensitive and controversal. When stated, most of the theoretical work makes use of the Dirichlet condition that the probability density p(r,t) vanishes at infinity. In fact, as numerical procedures ignore infinity, it is convenient to consider a second boundary condition located at a finite but large distance R’. Under these circumstances, the Dirichlet condition does not imply that the flux entering a sphere of radius R’ is zero, so some electrons would escape through boundary, R! Then, another assumption is that the exiting flux through R’is null. Such a condition means that, at R ’, the electron backscatters without being collected. It should be noted that both Dirichlet and null flux conditions are equivalent as R ‘approaches infinitySz0Also, eq 3 has been taken as the second boundary condition. The initial condition is given by eq 4 and supposes a uniform angular distribution of the electron. Given p(r,t), the surviving electron fraction N ( t ) is
N ( t ) = 4a
& rzp(r,t) dr R’
(5)
Moreover, in the presence of an electron scavenger, the scavenger fraction F ( S ) is given by26
F(S) = k ’ S J0m N ( t )exp(-k‘St) dt
(6)
where S is the scavenger concentration and k’ is the bimolecular rate constant for the electron-scavenger reaction. The mathematical problem is to find a function p(r,t) for r E [R,+m] and t E [O,+m] satisfying the set of eq 1-4. A first mathematicalstage%consists of changing the partial differential eq 1and conditions 2-4 by a single time differential equation in an appropriate functional space. The new differential equation is then equivalent to the initial problem and is obtained as follows. Let u(r) be a numerical trial function defined in [R,+a] and smooth in this interval to make the involved integrals convergent. After multiplying eq 1 by 4ar2v(r)and integrating over [R,+m], the use of boundary conditions 2 and 3 yields d -(p,v) + a(p,v) = 0 (7) dt where ( p p ) = 4a hmr2p(r,t)v(r) dr
(8)
represents a scalar product in the chosen functional space and
is a bilinear form. Equation 7 has to be satisfied by any trial function u(r) and constitutes with initial condition 4 the basic equation of the problem. Such a formulation is referred to as the variational method in mathematics.21 Then, for smooth enough distribution po(r), it has been recently provedz7that eq 7 possesses a unique smooth (24)(a) F.C.Collins and G. E. Kimball, J.Colloid. Sci., 4,425(1949); (b) R. M. Noyes, Prog. React. Kinet., 1, 129 (1961). (25)The generality of the so-calledradiation boundary condition has already been discussed in ref 10 and 14. (26)A. Hummel, J. Chem. Phys., 49,4840 (1968).
Numerlcal Solution of the Smoluchowski Equatlon
The Journal of Physlcal Chemistty, Vol. 85, No. 11, 1981
positive solution if po(r)is positive. This result shows that the mathematical formulation of the physical model is correctly stated. In a second stage the numerical solution of the problem (7,4) is obtained by a finite difference scheme including a relaxation coefficient 8. This coefficient serves to increase the precision of the solution without obeying drastic conditions due to the space and time discretization. The procedure is as follows. (i) The space interval is restricted to [R,R’I in such a way that the distance R’is sufficiently large to simulate an infinite distance. The interval [R,R1 is divided into a set of n points ranging from r1 = R to r,, = R ! The different spacings h, = rp+l- rp,p E [l,n - 11 do not have to be constant but small when compared to R in the vicinity of R. A spatial finite difference operator 8h is classically attributed to each of the n points. Then the eq 9 discretizes as ah(p,u) = 47rD
s,
R‘
(Ahp
+ Rcp)Ghudr + k d ( R , t ) u ( R )
(10) (ii) Similarly, the time interval[O,T] is divided into a set of (m + 1) points ranging from to = 0 to t , = T with spacings k , = ts+l- t,; s E [O,m- 11. Let 8 be a given real number in [0,1] and uo(r) a numerical function on [R,R’Iwhich is constant in the interval up and equal to po(r),where up = [r - h,-J2, r, + h,/2], with ho = h, = 0. The function uofr) is said to be a step function on the set of the intervals up. The problem is now to determine m step functions us(r),by means of the m recursive equations
(y,+ IJ)
ah((1 - 8)
+ 8u,,u) = 0
(11)
where s E [l,m] and by doing this for any step function u also defined on the set cP. It should be noted that any function u,(r) is a constant u S pin any specific nterval up. As shown re~ently,~’ the above discrete eq (11) also have a unique solution which obviously depends upon time and space numerical steps. In other words, the determination of u,(r) requires the solution of successive m systems of n linear nonhomogeneous equations by means of the tridiagonal matrix T =I +k,SA (12) where I is the unity matrix, S the matrix representing the discretized scalar product, and A the matrix associated to the discretized bilinear form. When h = sup h and k = sup k, become zero, the solution ul(r), ...,u,(r! that we denote u(r,t)and the integral
LR’
47rrzu(r,t) dr
respectively converge toward p(r,t) and N ( t ) (i) without any condition between h and k if the relaxation coefficient 8 belongs to [1/2, 11, (ii) with the condition k 1 - I C(R,Rc)(13) h2 (1- 28) if the coefficient 8 belongs to [0, 1/2[, C(R,Rc)is a known function. These conditions express the mathematical stability condition of the numerical algorithm. Inequality 13 is peculiarly drastic for 0 = 0. The singluar zero value corresponds to the so-called Freed-Pedersen28 or explicit finite difference scheme, i.e. (27) E. Croc, Thesis, Universite de Provence, 1979. (28) J. H. Freed and J. B. Pedersen, Adu. Magn. Reson., 8 , 1 (1976).
1551
whereas the so-called Crank-NicolsonB method is obtained for 8 = lI2. The above stability condition then suggests choosing a 8 value larger than lI2as a convenient numerical method, providing an unconditionally stable difference finite scheme. However, in addition to the stability condition, another problem is to obtain a numerical positive solution u,(r). As can be demonstratedz0the positive solution is unconditional for 8 = l ; therefore the so-called pure implicit finite difference scheme, i.e.
has been used in the present paper. All computations have been carried out from dimensionless variables and performed on a Univac 1110 computer. In what follows, the results are expressed in terms of the numerical solution u(r,t)with the requirement that u(r,t) absolutely converges toward p(r,t). 111. Numerical and Experimental Results The present model is relevant to liquids in which the solvated electron e; disappears by neutralization with the cation. Liquid hydrocarbons probably belong to this class of solvents. Now, for relatively low polarity solvents such as aliphatic amines, the disappearance process of e; depends on the nature of the amine,2r30at least in the homogeneous stage. For example, in ethylamine and npropylamine, there is some experimental evidence t o conclude that e; disappears mainly with the cation.z In ethylenediamine, the situation is complicated by the fact that the neutralization of e; with the cation competes with its combination with a neutral r a d i ~ a l By . ~ contrast, ~ ~ ~ in liquids like a m m ~ n i aand ~ l hydrazine,22 ~~~ the radical reaction dominates the decay process at long time. Although restricted to e;-cation isolated pairs, nevertheless the present model has been applied to the above liquids, even to those for which the e,:radical reaction intervenes at long time. This first approach allows one to estimate the extent to which reactions other than e;-cation neutralization reactions have to be considered in the nonhomogeneous stage. The available experimental data from pulse radiolysis are the time-dependent radiolytic yield G,;(t) and the scavenger-concentration-dependent product yield G (S). Given the initial solvated electron yield Goe8-,the agove quantities can be normalized. The surviving electron fraction N ( t ) is then defined as G,;(t)/Goeg-and the scavenged fraction F(S)as G,(S)/GoeI-. Experimental values of G ( t ) and N ( t ) are given at times t in the last two columns of Table I. It should be noted that only selected values of G,(t) have been reported. Indeed, due to generally high dose rates delivered in pulse irradiations, experimental values of G(t) are independent of the dose rate only at short times after the end of the pulse, otherwise they must be corrected for interspur decay in order to be compared with the isolated pair model. By contrast, scavenging experiments provide a great many experimental points. (29) J. Crank and P. Nicolson, Proc. Cambridge Phil. SOC., 43, 50 (1935). (30) J. A. Delaire, these de Doctorat d’Etat, Universitg de Paris XI, 1980. (31) J. Belloni, P. Cordier, and J. A. Delaire, Chem. Phys. Lett., 27, 241 (1974). (32) J. Belloni, F. Billiau, P. Cordier, J. A. Delaire, and M. 0. Delcourt, J. Phys. Chern., 82, 532 (1978).
1552
The Journal of Physical Chemistry, Vol. 85, No. 11, 1981
Delaire et ai.
T A B L E I: Phvsical Parameters a n d ExDerimental Values dielectric visc o n - cosity, stant CP Rc,A
D, c m 2 s-' 3.5 X low4
ha, L m o l - ' s-l
G( t ) 0.9 (10 ns)" 0.5 ( - ) d n-C,H,NH2 5.3 0.46 104 5.9 3.5 X lo6 < h a < m a 3.85a 0.8 (10 ns)" 0.3 (-)" C,H,(NH,), 14 1.54 40.4 6.0 4.0 X 10.' 3 X l o 6 < h a < 1.1X l o 7 " 3.8" 2.6 (10 ns)" 2.0 (300ns)" N2H4 54 0.97 10.5 5.3 4.3 x 2 x l o 7 < k .. < 3.75 X l o 7 4.OC 3.4 110 n s ) b 2.6 (300 n ~ a F r o m ref 1, 2, and 30. F r o m ref 22. Assumed value. R e f e r e n c e 42. e R e f e r e n c e 43. f Reference 44. G(t)IGoes-. solvent
7
CA"2
81
R, A 5.7
5.5 X
GO,,4.3"
l o 6 < ha < -"
It
c 2 H5 "2
N(tF 0.21 0.12 0.21 0.08 0.68 0.53 0.85 )0.65 ~ g N ( t )=
-I,
C2H5 "2
10
a
Flgure 1. Distribution function u ( r , t ) of e,--cation pairs in ethylamine for an initial Gaussian distribution ( P = 50 A) and a partially diffusioncontrolled recombination(k, = 1.8 X 10'' L mol-' s-'): (a) t = 0; (b) t = IO-'' s; (c) t = 2 X IO-" s; (d) t = 4 X IO-'' s; (e) t = 2 X io-1o S; (f) t = 1 0 - ~ S.
Parameters used in numerical calculations are also given in Table I. The reaction radius, R, is assumed to be the sum of the cation radius R+ and the e; radius R-; R+ has been obtained from partial molar volume measurements in solutionss3 and R- has been arbitrarily assumed to be equal to 3 A.34 The mutual diffusion coefficient D is the sum of the diffusion coefficients of the solvated electron D- and of the cation D+. The former has been calculated from experimental diffusion-controlledrate constants5 and the latter deduced from the Stokes-Einstein relation. The rate constant k, has been calculated from the experimental rate constant kexptof neutralization according to the r e l a t i o n ~ h i p ~ ~ -1- 1 - exp(-Rc/R) exp(-Rc/R)
-
+
(16)
4rDRc ka Due to the large experimental uncertainty in the determination of kexpt, Table I gives the lower and upper limits of ha. Moreover, for low polarity liquids and because kerpt
(33) J. E. Desnoyers and C. Jolicoeur, "Modern Aspects of Electrochemistry", Vol. 5, J. 0. M Bockris and B. E. Conway, Ed., Plenum Press, New York, 1969. (34) This value seems realistic in liquids such as ammonia and amines. An estimation of R-comes from the measurement of the partial molar volume of e, in metal-ammonia solutions (U.Schindewolf, "Colloque Weyl 11", Butterworths, London, 1970, p 199). (35) J. A.Delaire, M. 0. Delcourt, and J. Belloni, J.Phys. Chem., 84, 1186 (1980). (36) (a) S.R.Logan, J. Chem. SOC.,Faraday Tram., 63,1712 (1967); (b) ref 9, p 10.
-
:2.-.i b
:r:-f:: 0
I
I
I
I
,
,I c
The Journal of Physical Chemistry, Vol.
Numerical Solution of the Smoluchowski Equation 1.0
I
I
I
-
C2H5 “2
10
I
I
85,No. 17, 1981 1553
I
I
I
I
C2H5 “2
t \\ I
I
I
10-5
Q)
I
10-4
10-3 10-2 [Ph2] (mol. 1-l)
t (4
Figure 4. Survival probability in ethylamine: (a) i = 50 A; (b) i = 40 A. Full curve: k, = 1.8 X 10” L mol-’s-’; dashed curve: k, = 1.8 X 1 013 L mol-‘ s-’; (.- -) asymptotic values calculated according to eq 19. The vertical bars represent experimental values.
-
O
”
I
1
10-1
Figure 6. Scavenged fraction as a function of blphenyl concentration in ethylamine: k, = 1.8 X 10” L mol-’ s-’;(a) i = 60 A; (b) i = 50 A; (c) i = 40 A. (.-e -) asymptotic values calculated according to eq 19; ( 0 )experimental points from ref 2.
C2H5 “2
t L
0
4 I
I
10-5
(0-4
I
10-3
I
I
70-2
10-1 I Ph2J (md 1-1 )
1
Flgure 5. Scavenged fraction as a function of biphenyl concentration in ethylamine: k, = 1.8 X lo* L mol-’ s-l;(a) i = 50 A; (b) i = 40 A; (c) ? = 35 A. (. - -) asymptotic values calculated according to eq 19; (0)experimental points from ref 2.
-
P corresponds to the maximum of this Gaussian distribution. In ethylamine, Figures 1 and 2 show the probability density u(r,t)vs. r at different times for two values of k, and for P = 50 A. For this compound, the time evolution of the surviving fraction N ( t ) is shown in Figures 3 and 4 for different sets of the parameters k, and P. In the presence of biphenyl, the computed scavenged fraction F ( S ) is plotted from the experimental data of ref 2 as a function of the biphenyl concentration S in Figures 5 and 6 with the same set of parameters as in Figures 3 and 4. In Figure 7, F ( S ) has been computed for n-propylamine in the presence of biphenyl and is also compared with experimental values.2 An attempt was made to fit both k, and 7 by comparing calculated curves with experimental values of N ( t ) and F(S). For ethylamine, when varying both k, and F (Figures 3-6), no value of i; leads to a good fit for k, < los L mol-l s-l. As k, is raised from lo8 to 1013,the best value of P is found to range from 35 to 45 A. Then, in the case of ethylamine, for which the bulk recombination rate constant approaches the diffusion limit,2 the best value of P is practically independent of 12,. Similarly, in n-propylamine, the bulk rate constant is close to the diffusion limit2 and the best value found for P E 55 A is also practically independent of k,. For ethylenediamine, hydrazine, or ammonia where the bulk recombination rate constant is far from the diffusion limit,192it is interesting to find out if the isolated pair model predicts the observed values of the yield. For example, in the case of ethylenediamine, realistic values of P lead to a poor correlation between computed and experimental
Flgure 7. Scavenged fraction as a function of biphenyl concentration in n-propylamine: k, = L mol-‘ s-l; (a) i = 60 A; (b) i = 50 A; (.-.-) asymptotic values calculated according to eq 19; (0) experimental points from ref 2.
10-11
I
I
10-10
10-9
I
10-8
I
10-7
, I
I,
co
t (L)
Figure 8. Survival probability in ethylenediamlne: (a) k = 2.9 X lo7 L mol-’ s-’ i = 30 A; (b) k = 2.9 X IO7 L mol-’ s-’, = 10 A; (c] k, = 2.9 X IO8 L mol-’ s-l: i = 30 A; (d) 2_.9X 10’O L mols-l , r = 30 A; (e) k, = 2.9 x 10’’L mol- s , r = 30 A; ( . - e - ) asymptotic values calculated according to eq 19. The vertical bars represent experimental values.
r, ;
values of N ( t ) ,when taking k, within the uncertainty range given in Table I (full curves of Figure 8). However, a good fit is obtained when k, is raised by a factor of 10 (upper dashed curve of Figure 8). As this last value is out of the error range, we conclude that the model fails to give experimental N(t) values because of the competitive reaction of e, with the radical associated with the cation. In hydrazine and ammonia, the neutralization rate constant is roughly lo3 times smaller than the diffusion limit.’ In the case of hydrazine, the upper curve of Figure 9 shows that, even for small values of P, there is no decay of the computed N ( t ) when considering only the neu-
1554
The Journal of Physical Chemktty, Vol, 85, No. 11, 1987 r
I
I
I
Delaire et al.
I
I
%%"2
Ii Flgure 9. Survival probability in hydrazine. Curve (a) electron-cation recombination: k, = 7 X IO6 L mol-' s-l; i = 20 A. Other curves: electron-radical reactlon. (-) k, = loi3 L mol-' s-'; (---) k, = 10" L mol-' s-'. Curves (b) and (e): P = 25 A; curves (c) and (f): i = 20 A; curves (d) and (g): i = 15 A; (.-.-) asymptotic values calculated according to eq 19. The vertical bars represent experimental values.
tralization of the e[-N2H5+ pair. However, the observed experimental decay can be accounted for by applying the isolated pair model to the combination of the electron with the radical NzH3.under the following assumptions: (i) the mutual diffusive motion of the electron and the radical is not affected by the Coulombic field of the cation; (ii) the mutual diffusion coefficients and the reaction radii have the same value for the electron and cation as for the electron and the radical; (iii) the activation-controlled rate constant k, for electron-radical combination is so large that the reaction is diffusion controlled. The lower curves of Figure 9 show that, for a mean initial egradical distance of 20 A,the calculated curves N ( t ) satisfactorily overlap the experimental data. The radical initially forms close to the cation. When the Coulombic attraction of the cation is taken into account, both solvated electron and radical are attracted closer to each other. Under these conditions the rate of disappearance of the electron is enhanced, and the mean separation distance between the electron and the radical becomes larger than 20 A.
IV. Discussion The evolution of the probability density in Figures 1and 2 for ethylamine shows the great influence of k,. Beyond a distance roughly equal to twice the reaction radius, the evolution of the density profile is not really sensitive to the k, values and is dominated by the spreading of the initial Gaussian distribution and by the shift of the maximum toward the attractive cation. Such an evolution has also been obtained by Hong and Noolandi'O for an idealized sink, i.e., R = 0. By numerical integration, Hummel and Infelta" found a shift of the maximum away from the origin because they plotted the product, 4rr2u(r,t). At short distances, the profile drastically changes with k,. As expected, the probability density increases near the reaction radius for low values of k, whereas the density rapidly falls for high values of k,. For low k, values, the accumulation of opposite charges at small distances is due to the limitation of the reaction rate at encounter. This fact has already been pointed out by H ~ m m e l , ~ and ' raises the question of the nature of such a transient electron-cation assembly. The electron survival probability dependence on k, is shown in the case of ethylamine (Figures 3 and 4). For large k, with fast electron decay, ii must be large in order (37) Reference 9,p 28.
Flgure 10. Ultimate escape probability N(m)vs. the rate constant k, for different values of the posltion r o of the initial Dirac distribution function in the case of ethylamine.
to fit the experimental data. A careful inspection of the relative effect of k, and P on N ( t ) and then on F(S)shows that (i) for large values of k,,.N(t) no longer depends on k,, but only on ii (Figure 4); (11) for low values of k,, N ( t ) is poorly sensitive to ii (curves d-f of Figure 3). This behavior is expected from theoretical predictions on the long time behavior, N ( m ) , of the survival probability.13J4 For an initial Dirac distribution centered at ro, the escape probability N ( m ) vs. k, is: N(m) = 47rDRc exp(-RdR) + k,[exp(-Rc/ro) - exp(-Rc/R)I 4rDRc exp(-Rc/R) + k,[l - exp(-R,/R)] (18) This function is plotted vs. k, in Figure 10 for different values of ro for ethylamine. The asymptotic value of N(m) is rapidly reached even for values of k, two or three orders of magnitude smaller than the critical value 4rDRc which corresponds to the experimental bulk diffusion-controlled rate constant. In the asymptotic range, the escape probability is then predominantly dependent on ro whereas, for low values of k,, this probability is weakly affected by r p It should be noted that the probability N ( a ) exactly coincides with the Onsager escape probability, Le., exp(-Rc/ro) when the rate constant k, reaches the critical value 47rDRc or when the reaction radius R equals zero and 12, approaches infinity. Another interesting point emerges from the comparison of the long time behavior of the computed N ( t )with N ( m ) as given by eq 18. However, to make the comparison significant, the escape probability N(a)has been averaged over the initial Gaussian distribution g(ro,r) to yield ~~
~
IT(-)= 4 r L E ro2N(~,ro)g(ro;t) dro
(19)
For the sake of comparison, computed values of m(m) are reported on the N(t) plots in Figures 3,4, 8, and 9 and on the F(S) plots in Figures 5-7. After a prompt decay, the time evolution of N ( t )in the 10-'-104-s range slows down and then very smoothly goes to zero, with a rate slightly dependent on the value of R'numerically chosen to simulate infinity. It is noteworthy that N ( t ) at lo4 s is practically equal to the asymptotic value, N ( m ) . In this time range, the surviving charges are very nearly homogeneously distributed and undergo homogeneous slow decay. Experimentally, the slow decay generally obeys second-order kinetics. Such kinetics is not observed from computed curves because the model is valid for isolated pairs only and takes no account of bulk recombinations.
J. Phys. Chem. 1981, 85, 1555-1558
In ethylamine and n-propylamine for which the isolated pair spur model is valid, the t values of the initial Gaussian distance distribution have been obtained by both fittings of N ( t ) and F ( S ) and are 40 f 10 A and 50 f 10 A, respectively. From a Dirac initial distance distribution function centered at ro and from the experimental value , ~ Rice and Baird relation14 yields 27 and 40 of N ( c o ) the A, respectively. Both sets of values may be compared, keeping in mind that the discrepancy is due to the methods of their calculation. In these amines of low polarity, 7 can be also compared to the thermalization range of the solvated electron found to be about 40 A in ether^.^^,^^ By contrast, in water t is either 2340or 27 A41depending upon (38) J. P. Dodelet and G. R. Freeman, Can. J. Chem., 53, 1263 (1975). (39) F. Billiau, J. Belloni, J. A. Delaire, and M. 0. Delcourt, J.Chim. Phys., 76, 1059 (1979). (40) H. A. Schwarz, J. Chem. Phys., 55,3647 (1971). (41) A. Mozumder, J. Chem. Phys., 50, 3153 (1969). (42) W. A. Seddon, J. W. Fletcher, and F. C. Sopchyshyn, Can. J. Chem., 56, 839 (1978).
1555
the model chosen. Thus the r found for ethylamine and n-propylamine conform with the general idea that, the lower the polarity of the medium, the higher the mean distance between initial ion pairs. To the extent that the experimental conditions give isolated pairs, the numerical solution of the one-pair Smoluchowski equation fits well with the experimental data as shown for ethylamine and n-propylamine. For liquids for which the electron undergoes a competitive reaction with the radical, a new model has to be chosen. This work is in progress. Acknowledgment. This work has been partially supported by the financial specific program A.T.P. 2988 under the auspices of CNRS. The authors are pleased to thank Dr. P. Penel and Dr. A. M. Pelissier of the State University of Toulon for helpful discussions of the mathematical part. (43) W. A. Seddon, private communication. (44) W. A. Seddon, J. W. Fletcher, and F. C. Sopchyshyn, Can. J. Chem.,54, 2807 (1976).
Sensitivity Analysis of Oscillating Reactions. 1. The Period of the Oregonator David Edelson” and Valerie M. Thomast Bell Laboratories, Murray Hill, New Jersey 07974 (Received: December 8, 1980; In Final Form: February 20, 1981)
A method is developed for the computation of the dependence of the period of an oscillating reaction mechanism upon the input parameters of the kinetic model. By expressing the solution of the mass-action differential equations as a delay equation, one may relate the sensitivity coefficients of the period to those of the chemical species concentrations; the latter can be computed by any of a number of established techniques. Application is made to the five-step Oregonator model proposed by Field and Noyes to represent the essential characteristics of the complex Belousov-Zhabotinsky reaction. The method is validated by comparing the results with those obtainable by direct parameter variation. The calculated period sensitivities confirm quantitatively the relative significance of the various components of the model which previous workers could only suggest on the basis of approximate analyses.
The mechanisms which have thus far been conceived to account for the behavior of real oscillating chemical systems,l such as the Bray-Liebhafsky,2 Belousov-Zhabot i n ~ k y , ~and - ~ other bromate systems6J are composed of reaction steps for which the level of independent supporting information has the widest possible range. Some reactions, primarily the inorganic ones, have been measured independently so that their rate parameters are assigned with a fair degree of certainty. Others are estimated from thermodynamic data. The basis of the organic reactions, however, ranges from reasonable confidence to pure speculation, and even the participation of some of the proposed radical intermediates is not assured. Although the mechanisms manage to simulate the gross features of the reacting systems, significant discrepancies exist in detail; furthermore, the uniqueness of the assignment always remains in doubt. The interactions themselves are so complex that the kineticist’s traditional intuition is unsuccessful in placing responsibility upon particular components of the mechanism for an observed feature of the experiment. f
Swarthmore College, Swarthmore, PA 19081. 0022-3654/81/2085-1555$01.25/0
In the past decade techniques for the study of complex chemical reactions by mathematical modeling have advanced rapidly.* These methods are mainly based on the solution of the mass-action differential equations derivable from the chemical model and also provide a route for the study of the effect of the parameters of the model upon its behavior. The use of this “sensitivity analysis” was proposed many years agogand has been well developed for equations capable of analytic solution.1° Numerical extensions with particular application to chemistry have recently been advanced.l1-l3 These could be used to gain (1) Noyes, R. M. Ber. Bunsenges. Phys. Chem. 1980,84, 295. (2) Edelson, D.; Noyes, R. M. J. Phys. Chem. 1979, 83, 212. (3) Field, R. J.; Koros, E.; Noyes, R. M. J.Am. Chem. SOC.1972,94, 8649. (4) Edelson, D.; Field, R. J.; Noyes, R. M. Int. J. Chem. Kinet. 1975, 7, 417. (5) Edelson, D.; Noyes, R. M.; Field, R. J. Znt. J. Chem. Kinet. 1979, 11, 155. (6) Noyes, R. M. J.Am. Chem. SOC. 1980,102, 4644. (7) Orbin, M.; Koros, E. J.Phys. Chem. 1978,82, 1672. (8) Edelson, D. J. Chem. Educ. 1975,52, 642. (9) Poincar6, H. “New Methods of Celestial Mechanics”, 1892-1899; NASA Technical Translation TT F-450-2, 1967. (10) Hille, E. ”Lectures on Ordinary Differential Equations”; Addison-Wesley: Menlo Park,CA, 1969; Chapter 3.
0 1981 American Chemical Society
